Interpolating in a spherical triangle

Hello gurus,

I'm not sure if this is the best place to start with these questions, because they're really more related to mathematics, but I've got a good feeling y'all will help. Right now I'm really looking more for some input on how I'm planning to approach the problem rather than the code itself. I'd like to give people the opportunity to suggest alternatives before I pour my heart into this.... I have searched far and wide for similar questions/solutions, but nothing seems to fit quite right.

Here's the scenario:

I have a set of about 200 Lat/Long points each with an associated value of some magnitude, B. These points represent "stations" where data is collected.

The desire is for an end user to enter in any Lat/Long point, Pw, and have a value B returned that is interpolated from the known values....below is how I plan to interpolate this value:

1. Calculate the distances from each data point to Pw using the haversine formula

2. Define data points P1 and P2 as the two closest points to Pw

3. Define data point Pn as the next closest point that satisfies the criteria: Pw is in triangle P1P2Pn (note: I plan to use haversine distances for proximity, but lat long values to test that Pw is enveloped, is there anything wrong with that logic?)

4. Calculate the distances P1P2, P2P3, P1P3, P1Pw and P2Pw using the haversine formula

5. "flatten" the triangle (generate cartesian coordinates) that honor the distances calculated in step 4 (I think this skews distance P3Pw, not sure how problematic that is)

6. Now that my triangle is flat, and all points are located, use barycentric coordinates and interpolate B for Pw


How convoluted is this? Is there any accuracy to this? What would be the magnitude of error if all my data points are in the same state and I just used Lat/Long coords and jumped to barycentric coordinates?

My search for answers keeps leading me to dissertations and abstracts as if this is some unsolved paradox, but it seems like the answer should be simple... (simpler than my logic)


ANY help is appreciated, and I will post code as I write it once I get some feedback.

Thanks guys!

MarvinP, I appreciate the strategy. Solving your example problem with pencil and paper helped me to understand the barycentric component of the problem, and for what it's worth, I arrive at 5/3 as the final answer. My biggest qualm with keeping it this simple are units. Using lat and long as if they are Cartesian values doesn't seem wise. Because I'm enjoying the learning experience, I'd like to keep this spherical, and have started to code it with the same basic approach I previously outlined. I made it pretty far before my first snag. Please see the attached spreadsheet.

Here's my new scenario:

1. A haversine UDF is used to calculate the distance from each data point to my target point in a helper column

2. Another UDF for testing if a point is in the triangle created by the three data points is defined

3. a Third UDF defines the station names and coordinates, using the 2 closest data points, while finding a third that satisfies the triangle test.

My problem is, in my code, if I define i as any value and delete the Do Until Loop then I get expected results. As soon as I add the Loop statement and try to call the PTINPOLY function I get #VALUES.

What's wrong with my Loop?

Function HAVERSINEDIST(proj_lat, proj_long, stat_lat, stat_long)

earth_radius = 3959 'miles

Pi = 3.14159265

deg2rad = Pi / 180

dLat = deg2rad * (stat_lat - proj_lat)

dLon = deg2rad * (stat_long - proj_long)

a = Sin(dLat / 2) * Sin(dLat / 2) + Cos(deg2rad * proj_lat) * Cos(deg2rad * stat_lat) * Sin(dLon / 2) * Sin(dLon / 2)

c = 2 * WorksheetFunction.Asin(Sqr(a))

d = earth_radius * c

HAVERSINEDIST = d

End Function

Function PTINPOLY(stat_coords As Variant, proj_lat As Double, proj_long As Double) As Boolean

  Dim x As Long, m As Double, b As Double, poly As Variant, num_sides_crossed As Long
  
  poly = stat_coords
  
  For x = 1 To UBound(poly) - 1
  
    If poly(x, 1) > proj_long Xor poly(x + 1, 1) > proj_long Then
    
      m = (poly(x + 1, 2) - poly(x, 2)) / (poly(x + 1, 1) - poly(x, 1))
      
      b = (poly(x, 2) * poly(x + 1, 1) - poly(x, 1) * poly(x + 1, 2)) / (poly(x + 1, 1) - poly(x, 1))
      
      If m * proj_long + b > proj_lat Then num_sides_crossed = num_sides_crossed + 1
      
    End If
    
  Next
  
  PTINPOLY = num_sides_crossed Mod 2
  
End Function

Function STATIONS(proj_lat As Double, proj_long As Double) As Variant
   
    Dim stat_names(2) As Variant
    
    Dim stat_coords(2, 1) As Variant
    
    i = 3
    
        With WorksheetFunction
        
            stat_names(0) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Name
            
            stat_coords(0, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Latitude
            
            stat_coords(0, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Longitude
            
            stat_names(1) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station
            
            stat_coords(1, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station Latitude
            
            stat_coords(1, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station Longitude
            
            Do Until PTINPOLY(stat_coords, proj_lat, proj_long) = True
            
            stat_names(2) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station
            
            stat_coords(2, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station Latitude
            
            stat_coords(2, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station Longitude
        
            i = i + 1
        
            Loop
        
        End With
        
STATIONS = WorksheetFunction.Transpose(stat_names)
        
End Function

BDE Testing.xlsm

Can someone have a look, please? If I can troubleshoot the Loop I think I can finish this thing!

Thanks for the link, shg. Looks like I have some studying to do! In the meantime, could you please look at the Do Until Loop... Even if my approach is flawed, this is eating me up. Each UDF works independently, and I can manually change i until I arrive at a solution, but I can't get it to iterate without #VALUE!s

I think it's a simple syntax error, but I'm guessing at this point.

Function STATIONS(proj_lat As Double, proj_long As Double) As Variant
   
    Dim stat_names(2) As Variant
    
    Dim stat_coords(2, 1) As Variant
    
    i = 3
    
        With WorksheetFunction
        
            stat_names(0) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Name
            
            stat_coords(0, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Latitude
            
            stat_coords(0, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 1), ['BDE VALUES'!$E$4:$E$178], 0)) '1st Closest Station Longitude
            
            stat_names(1) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station
            
            stat_coords(1, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station Latitude
            
            stat_coords(1, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], 2), ['BDE VALUES'!$E$4:$E$178], 0)) '2nd Closest Station Longitude
            
            Do Until PTINPOLY(stat_coords, proj_lat, proj_long) = True
            
            stat_names(2) = .Index(['BDE VALUES'!$A$4:$A$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station
            
            stat_coords(2, 0) = .Index(['BDE VALUES'!$C$4:$C$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station Latitude
            
            stat_coords(2, 1) = .Index(['BDE VALUES'!$D$4:$D$178], .Match(.Small(['BDE VALUES'!$E$4:$E$178], i), ['BDE VALUES'!$E$4:$E$178], 0)) 'ith Closest Station Longitude
        
            i = i + 1
        
            Loop
        
        End With
        
STATIONS = WorksheetFunction.Transpose(stat_names)
        
End Function

....full code in previous post