Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

In an answer to this question ustroetz kindly provided some sample code for use of ogr to do lat/lon distance and point in polygon calculations.

I've created a Python class based on his example to do distance calcs but the results aren't correct. I realize that I could simply use haversine but the class is meant to grow to include other ogr functions, this is just the first step.

Allowing for the small difference due to the earth not being a perfect sphere, I'd expect the distance between degrees of latitude to be approximately 60 nautical miles anywhere from the equator to the poles, but the code in my class is returning approx 60 nm at the equator and 89 nm at 47 degrees latitude.

I'm geo-ignorant so don't know how to go about fixing this. I assume that the problem is in the class init where the coordinate transform is set up. Perhaps EPSG 3857 is not the correct projection to convert to from EPSG 4326?

Since distance is not calculating correctly I think that a point in polygon method (to be added to the class) will also be inaccurate unless I perform a repair.

#!/usr/bin/python

from __future__ import print_function
from osgeo import ogr
import osr
import sys

class Geo(object):
    def __init__(self):
        """
        Set up coordinate transform
        """
        inSpatialRef = osr.SpatialReference()
        inSpatialRef.ImportFromEPSG(4326)           # WGS 84
        outSpatialRef = osr.SpatialReference()
        outSpatialRef.ImportFromEPSG(3857)          #Spherical mercator
        self.coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

    def distance(self, position_1, position_2, units="km"):
        """
        Returns distance in 'units' (default km) between two Lat Lon point sets
                position_1 and position_2 are dicts containing "lon" and "lat"

                Units (optional) can be "km", "sm", or "nm"
        """
        point1 = ogr.Geometry(ogr.wkbPoint)
        point1.AddPoint(position_1["lon"], position_1["lat"])
        point1.Transform(self.coordTransform)

        point2 = ogr.Geometry(ogr.wkbPoint)
        point2.AddPoint(position_1["lon"], position_2["lat"])
        point2.Transform(self.coordTransform)

        raw_dist = point2.Distance(point1)

        print("Raw distance from ({0:6.2f}, {1:6.2f}) to ({2:6.2f}, {3:6.2f}) is {4:14.6f} km, "
            .format(position_1["lon"], position_1["lat"],
                    position_2["lon"], position_2["lat"],
                    raw_dist), end="")

        if units.lower() == "km":
            return raw_dist * 0.001
        elif units.lower() == "sm":
            #return raw_dist * 0.00062137
            return (raw_dist / 1000) * 0.621371192
        elif units.lower() == "nm":
            #return raw_dist * 0.000539956803
            return (raw_dist / 1000) * 0.539956803

if __name__ == "__main__":
        g = Geo()
        print("{0} {1}".format(g.distance({"lon":147.0, "lat":0}, {"lon":147.0, "lat":1}, "nm"), "nm"))
        print("{0} {1}".format(g.distance({"lon":147.0, "lat":47}, {"lon":147.0, "lat":48}, "nm"), "nm"))

The output produced by the code below is:

Raw distance from (147.00,   0.00) to (147.00,   1.00) is  111325.142866 M, 60.1107682357 nm
Raw distance from (147.00,  47.00) to (147.00,  48.00) is  164780.762454 M, 88.9744936905 nm
share|improve this question
    
Suggest showing the actual input and output. –  BradHards Apr 18 at 22:53
    
Some possible info here: osgeo-org.1560.x6.nabble.com/… –  mkennedy Apr 18 at 23:48
    
output added to question per Brad's suggestion –  RoyHB Apr 19 at 0:48

1 Answer 1

up vote 2 down vote accepted

If you are trying to calculate correct distance by projecting geographic coordinates (lon-lat) to map coordinates (E-N) and then using plane coordinate geometry, then you

  • either need to project using an appropriate equidistant map projection -- probably not a good idea if you have points at arbitrary locations and distances in arbitrary directions
  • or need to use a conformal projection -- which you are doing (Mercator) -- and apply the appropriate scale factor.

For long lines and accurate results, you may need to apply two or three scale factors. See Calculating distance scale factor by latitude for Mercator. At latitude 47°, the SF is 1 / cosine (47°) = approx 1 / 0.68 which, applied to your distorted-by-Mercator distance of 89NM, is now approx 60NM, i.e., a correct distance!

share|improve this answer
    
I'd suggest testing our combined procedure at a couple more latitudes to see if we're consistent. –  martin f Apr 19 at 4:15
    
Martin - I've looked at the page that you linked and can't see how to calculate the scale factor to arrive at the 0.68 you mention. Is there a lookup table somewhere or a simple means of calculating the SF? –  RoyHB Apr 19 at 5:49
    
@RoyHB - Sorry, i skipped a step or two! SF = 1 / cos (lat). Globe dist = map dist / SF = 89 * 0.68 –  martin f Apr 19 at 16:02

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.