I need to do the following in Python

  • open OSM & GPX files (I have packages for this)
  • transform points from GPX (like stick a track to roads)
  • calculate the results (distances, cumulative distances, etc.)

I need to do this repeatedly, so a Python script is much more preferred. Desktop software with plugins is not suitable. I'll prefer some Python & C modules than a quest of installing plugins. PostGIS may be an option too.

Shapely (Python package) seems to not be able to do this (it works only in 2D on a plane, and mentions it has no projections). I don't need 3D, but I have the input as lat&lon coordinates, and need to do geometric transformations (project a point on a polyline) and calculate distances in metres.

What modules should I use?

  • Maybe elaborate on which 3D functionality you need. It's not immediately clear.
    – underdark
    Commented Jun 30, 2013 at 8:25
  • Explained that in the last paragraph. The thing is that the input is lat&lon points, and the output should be metres. I can't see how Shapely does that.
    – culebrón
    Commented Jun 30, 2013 at 8:40

1 Answer 1

  • osgeo.ogr can read all these formats: OGR Vector Formats
  • osgeo.ogr and shapely support 3D:

    from osgeo import ogr
    point = ogr.Geometry(ogr.wkbPoint25D)
    from shapely.geometry import Point
    point1 = Point(5,4,4)
  • you can change projections with osgeo.ogr: see Projecting shapefile with transformation using OGR in python and many, many other examples

  • transform the geometries between ogr and shapely is easy:

    from shapely.wkb import loads
    point = ogr.Geometry(ogr.wkbPoint25D)
    point_shapely = loads(point.ExportToWkb())
  • inverse

    point_ogr = ogr.CreateGeometryFromWkb(point_shapely.wkb)
    print point_ogr.GetX(), point_ogr.GetY(), point_ogr.GetZ()
    5.0 4.0 0.0
  • so you can use ogr or pyproj to change the projection of a shapely geometry, (see Measuring distance in spherical Mercator vs zoned UTM for example)
  • and shapely or analytical geometry allows to project a point on a PoLyline (see How to draw perpendicular lines in QGIS?, with PyQGIS, but it is similar with ogr)

As one example of the process, here are the results of the creation of geological cross-sections from 3D points (from Python: Using vector and raster layers in a geological perspective, without GIS software, in French, but the scripts and the figures are universal).

3D representation (distance between points):

enter image description here

cumulative distance (geological cross-section)

enter image description here

  • I'll try doing thsi. But I have tested the latest version of Shapely: Point(0, 0, 0).distance(Point(0, 0, 1)) == 0.0 They explicitly say the z-value is not used in calculations.
    – culebrón
    Commented Jun 30, 2013 at 17:14
  • 2
    yes but you can use ogr for that or the module shapely.affinity allowing 3D affine transformations
    – gene
    Commented Jun 30, 2013 at 17:26

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.