6

I'm looking for a Python module that can provide some / all of the following services using latitude/longitude to define points:

Distance between two points
Point in polygon 
Distance along a multi-point path 
Identifier of closest point from a dict of points with identifiers 
Distance from a point to each point in a dict containing points with identifiers
etc

You get the idea... Any suggestions?

2 Answers 2

12

You can use the GDAL Python bindings. Examples on how to use it can you find here.

For example you create points with lat/lon like this

from osgeo import ogr                   # first import the library
point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(13.381348,52.536273)    # Berlin

point2 = ogr.Geometry(ogr.wkbPoint)
point2.AddPoint(11.557617,48.136767)    # Munich

enter image description here

Create a transformation from EPSG:4326 (lat/long) to EPSG:3035 (Projected coordinate system for Europe in meters)

inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4326)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(3035)
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

Transform our points

point1.Transform(coordTransform)
>> POINT (4550348.379724434576929 3275002.535206703934819 0) 

point2.Transform(coordTransform)
>>POINT (4436977.705661337822676 2781579.793173507787287 0)

And get the distance like this

point1.Distance(point2)             # Distance in meter from Munich to Berlin
>>> 506279.480221                   # roughly 506 km

Or you can create a polygon

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-23.378906,68.974164)     # North West Corner of Europe
ring.AddPoint(-23.378906,34.307144)     # South West Corner of Europe
ring.AddPoint(31.464844,34.307144)      # South East Corner of Europe
ring.AddPoint(31.464844,68.974164)      # North East Corner of Europe
ring.AddPoint(-23.378906,68.974164)     # North West Corner of Europe (to close to polygon)

polygon = ogr.Geometry(ogr.wkbPolygon)  
polygon.AddGeometry(ring)               # rough bounding box around Europe
polygon.Transform(coordTransform)       # transform it

enter image description here

And check if the polygon contains the point

polygon.Contains(point1)                    # Does Europe contain Berlin?
>>> True                                    # It does ;)

enter image description here

And for the rest you can write functions to do what you want.

4
  • Thanks. What are the units that are used in your examples? Do I need to convert my lat/lon points to these units? Or do I need to inform ogr that I'll be providing lat/lon data and it will take care of it?
    – RoyHB
    Commented Apr 2, 2014 at 16:03
  • Ah - I've discovered that the points are expressed as wkb - now I just need to figure out how to convert to from okb/lat lon
    – RoyHB
    Commented Apr 2, 2014 at 16:23
  • OGR doesn't do geometry on the globe, so it's a little more complicated than shown above. Coordinates first need to be projected into an appropriate 2D CRS.
    – sgillies
    Commented Apr 2, 2014 at 17:04
  • I edited my answer a little bit regarding your comments. I hope it is more clear now.
    – ustroetz
    Commented Apr 3, 2014 at 7:48
1

GDAL/OGR should be enough for you.

http://www.gdal.org/ogr/classOGRGeometry.html

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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