Skip to main content
added 21 characters in body; edited title
Source Link
PolyGeo
  • 65.4k
  • 29
  • 114
  • 345

My class using Using ogr to calculate distance is returning incorrect results away from the equator

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

My class using ogr to calculate distance is returning incorrect results away from the equator

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.

Using ogr to calculate distance is returning incorrect results away from equator

In an answer to Looking for Python module to do lat lon calculations @ustroetz provided some sample code for use of ogr to do lat/lon distance and point in polygon calculations.

replaced http://gis.stackexchange.com/ with https://gis.stackexchange.com/
Source Link

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

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.

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.

Tweeted twitter.com/#!/StackGIS/status/457628245681512448
added tags (including new 'scale-factor') moved code output to bottom (away from code itself)
Source Link

inIn 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'm (obviously) 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.

Help much appreciated.

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

#!/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

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'm (obviously) 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.

Help much appreciated.

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

#!/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"))

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'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
deleted 2 characters in body
Source Link
RoyHB
  • 332
  • 2
  • 9
Loading
added 247 characters in body
Source Link
RoyHB
  • 332
  • 2
  • 9
Loading
Source Link
RoyHB
  • 332
  • 2
  • 9
Loading