3

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?

RyanKDalton
  • 23,068
  • 17
  • 110
  • 178
culebrón
  • 2,314
  • 2
  • 23
  • 32
  • Maybe elaborate on which 3D functionality you need. It's not immediately clear. – underdark Jun 30 '13 at 08: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 Jun 30 '13 at 08:40

1 Answers1

2
  • 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)
point.AddPoint(5,4,4)
point.GetZ()
4.0

from shapely.geometry import Point
point1 = Point(5,4,4)
point1.has_z
True
point1.z
4.0

  • 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.AddPoint(5,4,4)
point_shapely = loads(point.ExportToWkb())
point_shapely.has_z
True

  • 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

Community
  • 1
gene
  • 54,868
  • 3
  • 110
  • 187
  • 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 Jun 30 '13 at 17:14
  • 2
    yes but you can use ogr for that or the module shapely.affinity allowing 3D affine transformations – gene Jun 30 '13 at 17:26