Announcement

Collapse
No announcement yet.

distance between 2 points using lat/lon

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • distance between 2 points using lat/lon

    I have a JavaScript function which will tell you the distance between 2 points using latitude/longitude

    Code:
    // figure distance between 2 lat/lon points
    function distance(lat1, lon1, lat2, lon2, unit) {
        if ((lat1 == lat2) && (lon1 == lon2)) {
            return 0;
        }
        else {
            var radlat1 = Math.PI * lat1/180;
            var radlat2 = Math.PI * lat2/180;
            var theta = lon1-lon2;
            var radtheta = Math.PI * theta/180;
            var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
            if (dist > 1) {
                dist = 1;
            }
            dist = Math.acos(dist);
            dist = dist * 180/Math.PI;
            dist = dist * 60 * 1.1515;
            // 'M' - miles is default
            if (unit=="K") { dist = dist * 1.609344 }
            if (unit=="N") { dist = dist * 0.8684 }
            if (unit=="F") { dist = dist * 5280 }
            return dist;
        }
    }
    Lat/lon are passed as numbers, like 41.39733579 and -90.33258726.

    I'm trying to convert to PB, but it doesn't work.

    I have:

    Code:
    pi = 3.141592653589793
    
    ' returns distance in rounded off feet
    Function distance(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Long
        Local radlat1 As Double
        Local radlat2 As Double
        Local theta As Double
        Local radtheta As Double
        Local dist As Double        
    
        If ((lat1 = lat2) And (lon1 = lon2)) Then
            Function = 0
            Exit Function    
        Else 
            radlat1 = pi * lat1/180
            radlat2 = pi * lat2/180
            theta = lon1-lon2
            radtheta = pi * theta/180
            dist = Sin(radlat1) * Sin(radlat2) + Cos(radlat1) * Cos(radlat2) * Cos(radtheta)       
    
            If (dist > 1) Then
                dist = 1
            End If
    
            dist = arccos(dist)           
            dist = dist * 180/pi
            dist = dist * 60 * 1.1515
    
            Function=dist*5280 ' return in feet
    
        End If
    End Function              
    
    Function arccos(v As Double) As Double
        Function = pi / 2 - Atn(v / Sqr(1 - v * v))          
    End Function

    what am I doing wrong?

    For context, you can do an online test here

  • #2
    Shawn,
    These are digital degrees.
    41.39733579 and -90.33258726. To stay digital you merely take the difference between the Lats (Delta X) and the difference between the Lons (Delta Y). That will give you the (Delta X/ Delta Y) slope of the arc. And the distance is determined by ABS(square root of [(Delta X) sq + (Delta Y) sq]) = Distance.

    Comment


    • #3
      Never mind... this code is good.
      The problem was elsewhere in my code.

      Comment


      • #4
        Though I'm interested in your idea... can you give me sample code? Your way sounds much simpler.

        Comment


        • #5
          Originally posted by Jim Fritts View Post
          Shawn,
          These are digital degrees.
          41.39733579 and -90.33258726. To stay digital you merely take the difference between the Lats (Delta X) and the difference between the Lons (Delta Y). That will give you the (Delta X/ Delta Y) slope of the arc. And the distance is determined by ABS(square root of [(Delta X) sqr + (Delta Y) sqr]) = Distance.
          True if you live on a flat plane, not if you live on the surface of a sphere (to a first approximation)

          Comment


          • #6
            Originally posted by Shawn Anderson View Post
            Though I'm interested in your idea... can you give me sample code? Your way sounds much simpler.
            There is no simple way. Even your code fails in many situations (try the link with a combination of E and W (+ and -) longitudes

            I've got a full set of VBA functions for Great Circle and Rhumb Line distance and bearing calculations in an Access database which contain Lat and Lon of over 12,000 locations worldwide. . I'll look at converting the functions to PB and posting them if you're interested.

            To give you a taste, the actual formula used for Great Circle Distance (in nautical miles) is:

            Code:
            GCDist = RadToNm(2 * arcsin(Sqr((Sin((lat1 - lat2) / 2)) ^ 2 + Cos(lat1) * Cos(lat2) * (Sin((lon1 - lon2) / 2)) ^ 2)))

            Comment


            • #7
              I'd definitely be interested in seeing your functions.
              We do a lot of stuff with the Google maps API which I'd love to not use when possible. Today I just needed compare distances between 2 sets of GPS coords.
              I had a list of about 600 coords that I need to compare to see which ones where within 1000 feet of another list of about 1800.
              PBCC read the data from CSV and brute-force matched them (over 1 million combinations) in about 30 seconds. It would have been faster with no screen output.


              Comment


              • #8
                I threw this together... Like Stuart points out it is a simple planar approximation.
                Code:
                #COMPILE EXE
                #DIM ALL
                
                FUNCTION PBMAIN () AS LONG
                    LOCAL SILAT AS EXT
                    LOCAL SILON AS EXT
                
                    LOCAL DILAT AS EXT
                    LOCAL DILON AS EXT
                
                    LOCAL DIFX  AS EXT
                    LOCAL DIFY  AS EXT
                
                    LOCAL DIST  AS EXT
                
                    'Springfield IL LAT 39.799358 LON -89.643623
                    'Decatur IL     LAT 39.841450 LON -88.955850
                    'Distance via the web apps is 41.6 to 38 miles(walking on roads) or AVERAGE 39.8 miles
                    SILAT = 39.799358
                    SILON = -89.643623
                    DILAT = 39.841450
                    DILON = -88.955850
                
                    DIFX = SILAT - DILAT
                    DIFY = SILON - DILON
                    '? "DELTA_X= " & USING$("###.##############", DIFX) & "  DELTA_Y= " & USING$("###.##############", DIFY)
                
                    DIST = ABS(SQR(DIFX * DIFX + DIFY * DIFY))
                    DIST = DIST * 100
                    DIST = DIST /2    '= 34.45
                    DIST = DIST * 1.155  'Variance
                    ? "The distance between Springfield and Decatur is" & USING$("###.##", DIST) & " miles."  'Result 39.79
                END FUNCTION
                Disclaimer:
                The code was not tested for other locations.

                NOTE:
                Web distances are usually aligned with road surfaces even while walking. To do it righteously you would need to incorporate elevation into the equation for every foot step.

                Ok break is over...

                Comment


                • #9
                  Here's some distance code that I've been using for awhile. I can't remember where I found it but it looks suspiciously like the Java code in the original post. I seem to recall that the output matched test vectors from a few web sites. No guarantees.

                  Code:
                  MACRO rad2deg (rad) = CDBL(rad * 180 / PI) ' Convert radians to degrees
                  MACRO deg2rad (deg) = CDBL(deg * PI / 180) ' Convert degrees to radians
                  
                  '//////////////////////////////////////////////////////////////////////////////////////////////////
                  ' zMapDistance - Returns distance between two locations.
                  '
                  ' Calulate distance based on longitude and latitude of both locations. String arg unit determines
                  ' return type (statute miles, kilometers or nautical miles).
                  '//////////////////////////////////////////////////////////////////////////////////////////////////
                  FUNCTION zMapDistance( _
                      BYREF lat1 AS DOUBLE, _ 'IN:  Location 1 latitude
                      BYREF lon1 AS DOUBLE, _ 'IN:    and longitude
                      BYREF lat2 AS DOUBLE, _ 'IN:  Same for location 2
                      BYREF lon2 AS DOUBLE, _ 'IN:    Ditto
                      BYREF unit AS STRING _  'IN:  M=statute miles K=kilometers N=nautical miles
                      ) AS LONG               'RTN: Integer distance between locations
                  
                    LOCAL theta, dist AS DOUBLE
                  
                    IF (lat1 <> lat2) OR (lon1 <> lon2) THEN
                      theta = lon1 - lon2
                      dist =  SIN(deg2rad(lat1)) * SIN(deg2rad(lat2)) + _
                              COS(deg2rad(lat1)) * COS(deg2rad(lat2)) * _
                              COS(deg2rad(theta))
                      dist = ACOS(dist)
                      dist = rad2deg(dist) * 60 * 1.1515
                      SELECT CASE UCASE$(unit)
                        CASE "M" : FUNCTION = dist
                        CASE "K" : FUNCTION = dist * 1.609344
                        CASE "N" : FUNCTION = dist * 0.8684
                      END SELECT
                    END IF
                  END FUNCTION  ' zMapDistance

                  Comment


                  • #10
                    Here you go. A quick conversion of the relevant VBA functions to PB with a simple demo including conversion of "Degree and decimal minute" strings to decimal degrees. (The database I have came with this common format)

                    Results are in Nautical Miles which is standard for sailors and pilots. Conversion to statute miles, yards, feet, kilometers, meters, furlongs etc. is a trivial exercise which is left to you

                    Code:
                    #COMPILE EXE
                    #DIM ALL
                    MACRO Pi = 3.141592653589793##
                    MACRO TwoPi = 6.283185307179586##
                    MACRO HalfPi = 1.5707963267948966##
                    MACRO QtrPi = 0.7853981633974483##
                    MACRO DegToRad(dpDegrees) = (dpDegrees*0.0174532925199433##)
                    MACRO RadToDeg(dpRadians) = (dpRadians*57.29577951308232##)
                    MACRO arcsin(X) = ATN(X / SQR(-X * X + 1))
                    MACRO acos(angle) = ATN(-angle / SQR(-angle * angle + 1)) + 1.5707963267949##
                    MACRO NmToRad(distance_nm)  =  .000290888208665722##  * distance_nm
                    MACRO RadToNm(distance_radians) = 3437.746770784939252## * distance_radians
                    
                    FUNCTION atan2(y AS DOUBLE, X AS DOUBLE) AS DOUBLE
                        'modification of ATN to deal with different quadrants. Used in Rhumb line calcuations.
                        SELECT CASE X
                        CASE IS > 0
                            FUNCTION = ATN(y / X) '     x>0
                        CASE IS = 0
                            SELECT CASE y
                                CASE IS > 0
                                    FUNCTION = HalfPi ' x=0, y>0
                                CASE IS < 0
                                    FUNCTION = - HalfPi 'x=0, y<0
                                CASE 0
                                    FUNCTION = 0      'x = 0, y = 0
                            END SELECT
                        CASE IS < 0
                             SELECT CASE y
                                CASE IS >= 0
                                    FUNCTION = ATN(y / X) + pi ' x<0, y>=0
                                CASE IS < 0
                                    FUNCTION = ATN(y / X) - pi 'x<0, y<0
                              END SELECT
                        END SELECT
                    END FUNCTION
                    
                    FUNCTION ModDbl(X AS DOUBLE, y AS DOUBLE) AS DOUBLE
                      'Modulus of floating point values - used in Rhumb Line calculations
                      IF y >= 0 THEN
                        FUNCTION = y - X * INT(y / X)
                      ELSE
                        FUNCTION = y + X * (INT(-y / X) + 1)
                      END IF
                    END FUNCTION
                    
                    
                    FUNCTION Latitude(sLat AS STRING) AS DOUBLE
                        'convert latitude string in form 9°27.9' S to latitude in decimal degrees
                        LOCAL d, m AS DOUBLE
                        LOCAL ns AS LONG
                        d = VAL(sLat) 'left numeric part
                        m = VAL(MID$(sLat, INSTR(sLat, "°") + 1))
                        ns = IIF(RIGHT$(sLat, 1) = "N", 1, -1)
                        FUNCTION = (d + m / 60) *ns
                    END FUNCTION
                    
                    FUNCTION Longitude(sLon AS STRING) AS DOUBLE
                        'convert  longitude string in form 147°11.5' E to longitude  in decimal degrees
                        LOCAL d, m AS DOUBLE
                        LOCAL ew AS LONG
                        d = VAL(sLon) 'left numeric part
                        m = VAL(MID$(sLon, INSTR(sLon, "°") + 1))
                        ew = IIF(RIGHT$(sLon, 1) = "E", 1, -1)
                        FUNCTION = (d + m / 60) * ew
                    END FUNCTION
                    
                    FUNCTION GCDist(la1 AS DOUBLE, lo1 AS DOUBLE, la2 AS DOUBLE, lo2 AS DOUBLE) AS DOUBLE
                        'Great Circle distance in Nautical miles
                        LOCAL lat1, lat2,lon1, lon2 AS DOUBLE
                        lat1 = DegToRad(la1)
                        lon1 = DegToRad(lo1)
                        lat2 = DegToRad(la2)
                        lon2 = DegToRad(lo2)
                        FUNCTION = RadToNm(2 * arcsin(SQR((SIN((lat1 - lat2) / 2)) ^ 2 + COS(lat1) * COS(lat2) * (SIN((lon1 - lon2) / 2)) ^ 2)))
                    END FUNCTION
                    
                    FUNCTION GCDir(la1 AS DOUBLE, lo1 AS DOUBLE, la2 AS DOUBLE, lo2 AS DOUBLE) AS LONG
                        'Great Circle Initial Bearing in degrees
                        LOCAL lat1, lat2 , lon1, lon2, Dist AS DOUBLE
                        lat1 = DegToRad(la1)
                        lon1 = DegToRad(lo1)
                        lat2 = DegToRad(la2)
                        lon2 = DegToRad(lo2)
                    
                        IF (COS(lat1) < 0.0000001) THEN '(a small number ~ machine precision) to deal with polar circles
                          IF (lat1 > 0) THEN
                             FUNCTION   = 180       '  starting from N pole
                          ELSE
                             FUNCTION   = 0        '  starting from S pole
                          END IF
                         EXIT FUNCTION
                        END IF
                        Dist = NmToRad(GCDist(la1,lo1,la2,lo2))
                        IF SIN(lon2 - lon1) < 0 THEN
                           FUNCTION = RadToDeg(acos((SIN(lat2) - SIN(lat1) * COS(Dist)) / (SIN(Dist) * COS(lat1))))
                        ELSE
                           FUNCTION = RadToDeg(TwoPi - acos((SIN(lat2) - SIN(lat1) * COS(Dist)) / (SIN(Dist) * COS(lat1))))
                        END IF
                    END FUNCTION
                    
                    FUNCTION RhumbDir(la1 AS DOUBLE, lo1 AS DOUBLE, la2 AS DOUBLE, lo2 AS DOUBLE) AS DOUBLE
                        'Rhumb Line direction in degrees
                        LOCAL lat1, lat2 , lon1, lon2, Dist AS DOUBLE
                        LOCAL dlon_W, dlon_E, dphi, q AS DOUBLE
                        lat1 = DegToRad(la1)
                        lon1 = DegToRad(lo1)
                        lat2 = DegToRad(la2)
                        lon2 = DegToRad(lo2)
                        dlon_W = ModDbl((TwoPi), (lon2 - lon1))
                        dlon_E = ModDbl((TwoPi), (lon1 - lon2))
                        dphi = LOG(TAN(lat2 / 2 + QtrPi) / TAN(lat1 / 2 + QtrPi))
                        IF (ABS(lat2 - lat1) < 0.000000000001) THEN
                          q = COS(lat1)
                        ELSE
                          q = (lat2 - lat1) / dphi
                       END IF
                       IF (dlon_W < dlon_E) THEN ' Westerly rhumb line is the shortest
                            FUNCTION = RadToDeg(ModDbl((TwoPi), (atan2(-dlon_W, dphi))))
                          ELSE
                            FUNCTION = RadToDeg(ModDbl((TwoPi), (atan2(dlon_E, dphi))))
                      END IF
                    END FUNCTION
                    
                    FUNCTION RhumbDist(la1 AS DOUBLE, lo1 AS DOUBLE, la2 AS DOUBLE, lo2 AS DOUBLE) AS DOUBLE
                        'Rhumb line distance in nautical miles
                        LOCAL lat1, lat2 , lon1, lon2, Dist AS DOUBLE
                        LOCAL dlon_W, dlon_E, dphi, q AS DOUBLE
                        lat1 = DegToRad(la1)
                        lon1 = DegToRad(lo1)
                        lat2 = DegToRad(la2)
                        lon2 = DegToRad(lo2)
                    
                        dlon_W = ModDbl(TwoPi, (lon2 - lon1))
                        dlon_E = ModDbl(TwoPi, (lon1 - lon2))
                        dphi = LOG(TAN(lat2 / 2 + QtrPi) / TAN(lat1 / 2 + QtrPi))
                    
                        IF (ABS(lat2 - lat1) < 0.00000001) THEN
                             q = COS(lat1)
                        ELSE
                             q = (lat2 - lat1) / dphi
                        END IF
                        IF (dlon_W < dlon_E) THEN ' Westerly rhumb line is the shortest
                              FUNCTION = RadToNm(SQR(q ^ 2 * dlon_W ^ 2 + (lat2 - lat1) ^ 2))
                        ELSE
                              FUNCTION = RadToNm(SQR(q ^ 2 * dlon_E ^ 2 + (lat2 - lat1) ^ 2))
                        END IF
                    END FUNCTION
                    
                    
                    FUNCTION PBMAIN() AS LONG
                        LOCAL la1,lo1,la2,lo2 AS DOUBLE
                        LOCAL sla1,slo1,sla2,slo2 AS STRING
                        LOCAL strResult AS STRING
                        sla1 = "9°27.9'S" :slo1 = "147°11.5'E" 'Port Moresby
                        sla2 = "51°28' N" :slo2 = "0°12' E" ' Greenwich
                        la1= latitude(sla1):lo1= longitude(slo1)
                        la2= latitude(sla2):lo2= longitude(slo2)
                        strResult = "From: Lat " & FORMAT$(la1,"0.000") & " Lon " & FORMAT$(lo1,"0.000")& $CRLF & _
                                    "  To: Lat " & FORMAT$(la2,"0.000") & " Lon " & FORMAT$(lo2,"0.000") &  $CRLF & $CRLF & _
                                    "Great Circle Distance: " & FORMAT$(gcdist(la1,lo1,la2,lo2),"0.0") & " NM" &  $CRLF & _
                                    "Great Circle Bearing: " & FORMAT$(gcdir(la1,lo1,la2,lo2),"000") & "°" &  $CRLF & _
                                    "Rhmb Line Dstance: " & FORMAT$(rhumbdist(la1,lo1,la2,lo2),"0.0") & " NM" & $CRLF & _
                                    "Rhumb Line Bearing: " & FORMAT$(rhumbdir(la1,lo1,la2,lo2),"000")  & "°"
                        ? strResult
                    END FUNCTION
                    Last edited by Stuart McLachlan; 9 Nov 2019, 05:55 AM.

                    Comment


                    • #11
                      Jerry,
                      Your function gets 37 miles.
                      'Distance via the web apps is 41.6 to 38 miles(walking on roads) or AVERAGE 39.8 miles

                      Comment


                      • #12
                        When did Jerry say anything about roads?
                        Dale

                        Comment


                        • #13
                          See Post #8. I actually got 34.45 but had to tack on a variance to get it to the road mileage.

                          Comment


                          • #14
                            Originally posted by Jim Fritts View Post
                            I threw this together... Like Stuart points out it is a simple planar approximation.
                            NOTE:
                            Web distances are usually aligned with road surfaces even while walking. To do it righteously you would need to incorporate elevation into the equation for every foot step.

                            Ok break is over...
                            If you are aligning with road surfaces, you need to incorporate every bend in the road too

                            How far is from the bottom to the top of this road?

                            Click image for larger version  Name:	5000311978_e15ae0008a_b.jpg Views:	0 Size:	385.2 KB ID:	786345

                            Comment


                            • #15
                              Originally posted by Jim Fritts View Post
                              See Post #8. I actually got 34.45 but had to tack on a variance to get it to the road mileage.
                              Too late! You should have said that then. And that still does not answer that you reported travel by road distance and compared it to straight line. If you want to say your code got a different straight line than Jerry's, then you should have said that. And given the subject is distance by latitude and longitude, the flat surface code is off topic anyway
                              Dale

                              Comment


                              • #16
                                Originally posted by Dale Yarker View Post
                                When did Jerry say anything about roads?
                                Well, he didn't. If you want driving/walking distance, a LOT of potential route data and some fancy logic will be required - as in Google Maps. I suspect that Google or someone else provides some type of service and an API to get that information. For me, I needed to convert IP addresses to longitude/latitude (also an available service) so that I could "connect" client applications to the nearest available server. Straight-line miles were fine and resolution was not a significant issue since a single server could support one or more states/provinces.

                                Comment


                                • #17
                                  I've just edited Post #10 to take advantage of PB's Macros. It should be much more efficient with reduced function calls and arithmetic operations.

                                  Comment


                                  • #18
                                    I have an application that requires range and bearing over short (less than 5 miles) distances. Over really short distances bearing in post #10 for rhumb and great circle point in opposite directions. Rhumb line seems to be correct for the few examples I tried, will that always be the case. (if necessary you can ignore crossing from E to W and N to S)

                                    Comment


                                    • #19
                                      Dale --> That's ok you are fine.
                                      Stuart --> Love the picture.

                                      NOTE:
                                      In my example the 2 sites were nearly the same Latitude. Also I did not consider the river drift current when crossing them or any head winds that would interfere with sure footedness. Ha Ha
                                      Last edited by Jim Fritts; 9 Nov 2019, 06:08 PM.

                                      Comment


                                      • #20
                                        Originally posted by James McNab View Post
                                        I have an application that requires range and bearing over short (less than 5 miles) distances. Over really short distances bearing in post #10 for rhumb and great circle point in opposite directions. Rhumb line seems to be correct for the few examples I tried, will that always be the case. (if necessary you can ignore crossing from E to W and N to S)
                                        What are "really short " distances? (I'll check the code above)

                                        If you are calculating over less than 5 miles, then plane geometry is simpler, faster and accurate enough for any practical purposes.

                                        Comment

                                        Working...
                                        X