Improving this script - GPX reading and processing with Python OGR

Question

I am trying to plot every GPX trace I have in a folder, by reading it with Python-OGR, clipping it by a bounding box, and plotting

Curently I have a working script, but it seems very uggly, and I suspect there should be a more straightforward way to achieve the same results.

Specific questions in comments:

#coding:utf-8

import glob
import ogr
import json
import numpy
import matplotlib.pyplot as plt


# 1) Is there a more direct or idiomatic way to create this rectangle?
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-51.22, -30.027)
ring.AddPoint(-51.22, -30.047)
ring.AddPoint(-51.18, -30.047)
ring.AddPoint(-51.18, -30.027)
ring.AddPoint(-51.22, -30.027)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)



for filename in sorted(glob.glob('../../../04 GPS/REGISTROS/*/*.gpx'), reverse=True)[:1]:

    # 2) Are all those "if not X: continue" statements the right way to check?
    datasource = ogr.Open(filename)
    if not datasource:
        continue

    layer = datasource.GetLayerByName('tracks')
    if not layer:
        continue

    feature = layer.GetNextFeature()  # Shouldn't I get *every* feature?
    if not feature:
        continue

    geometry = feature.GetGeometryRef().Intersection(poly)
    if not geometry:
        continue


    # 3) Exporting was the only way I found
    # to get numeric values. Is there any other?
    j = json.loads(geometry.ExportToJson())


    # 4) one of returned objects have only "geometry" key...
    if not 'coordinates' in j:
        continue

    for segment in j['coordinates']:
        seg = numpy.array(segment, ndmin=2)
        plt.plot(seg[:,0], seg[:,1], lw=0.5, alpha=0.1)


plt.axis('equal')
plt.show()

Answer 1

You have many solutions but you choose to use ogr and the GeoJSON format (you could have chosen Fiona and shapely or one of the many modules to parse gpx files: Pypi:gpx or to parse XML files).

1) Is there a more direct or idiomatic way to create this rectangle?

ogr is is verbose, look at Python GDAL/OGR Cookbook

With the GeoJSON format:

import json
polygs = {"type": "Polygon", "coordinates": [[[-51.22, -30.027, 0.0], [-51.22, -30.047, 0.0], [-51.18, -30.047, 0.0], [-51.18, -30.027, 0.0], [-51.22, -30.027, 0.0]]]}
poly = ogr.CreateGeometryFromJson(json.dumps(polygs))

With shapely and the wkt format:

from shapely.geometry import Polygon
polys = Polygon([(-51.22,-30.027, 0),(-51.22, -30.047, 0),(-51.18, -30.047, 0),(-51.18, -30.027, 0),(-51.22, -30.027, 0)])
poly = ogr.CreateGeometryFromWkt(polys.wkt)

2) Are all those "if not X: continue" statements the right way to check?

Use a single for loop to find the layers in a gpx file (replace your script for loops)

from osgeo import ogr
datasource = ogr.Open('my.gpx')
for i in range(datasource.GetLayerCount()):
     print datasource.GetLayer(i).GetName()," :", datasource.GetLayer(i).GetFeatureCount()     
 waypoints  : 0
 routes  : 0
 tracks  : 1
 route_points  : 0
 track_points  : 547

3) Exporting was the only way I found to get numeric values. Is there any other?

Examine the first layer with one feature: tracks: what is the geometry ?

 # geometry type
 layer = datasource.GetLayerByName('tracks')
 first = layer.GetFeature(0)
 geom = first.GetGeometryRef()
 print geom.GetGeometryName()
 MULTILINESTRING

a) Processing with the GeoJSON format: define two simple functions to plot the lines and the multilines formats from GeoJSON coordinates with matplotlib (you don't need the segments computation)

 def plot_line(ax, ob, col):
      x,y = zip(*[(j[0],j[1]) for j in ob])
      ax.plot(x, y, color=col)

 def plot_multiline(ax, ob,col):
      for line in ob:
         plot_line(ax,line,col)

Plot the feature

 # transform the geometry to GeoJSON format
 import json
 geomjson = json.loads(geom.ExportToJson())
 # plot the feature
 import matplotlib.pyplot as plt
 fig = plt.figure() 
 ax = fig.gca() 
 geomcoor = geomjson['coordinates']
 if geomjson['type'] == 'LineString':
     plot_line(ax, a,'g')
 if  geomjson['type'] == 'MultiLineString':
     plot_multiline(ax, a,'g')

Result:

If there are many features, use a simple for loop

b) Examine the second layer with features: tracks_points (Points) and simple processing without the GeoJSON format

layer2 = datasource.GetLayerByName('track_points')
fig = plt.figure() 
ax = fig.gca() 
for i in range(layer2.GetFeatureCount()):
     point = layer2.GetFeature(i)
     geom_point = point.GetGeometryRef()
     x,y,z = geom_point.GetPoint()
     ax.scatter(x,y)

Result:

It is easy to adapt these solutions with:

geometryfeature.GetGeometryRef().Intersection(poly) 

New: solution with Fiona and shapely

import fiona
fiona.listlayers('my.gpx')
[u'waypoints', u'routes', u'tracks', u'route_points', u'track_points']
# open layer tracks
layer = fiona.open('my.gpx', layer='tracks')
geom = layer[0] # only one feature here)
geom['geometry']['type']
'MultiLineString'
# convert to shapely geometry
from shapely.geometry import shape, mapping
tracks = shape(geom['geometry'])

Polygon

polygs = {'type': 'Polygon', 'coordinates': [[(-79.159610379971326, 35.615399458474776), (-79.153546268732796, 35.615411586697256), (-79.153109652723629, 35.611215221720187), (-79.159392071966749, 35.611118195940371), (-79.159610379971326, 35.615399458474776)]]}
poly = shape(polygs)

Intersection:

geometry = poly.intersection(tracks)
mapping(geometry)['type']
'MultiLineString'

Plot:

import matplotlib.pyplot as plt
from descartes.patch import PolygonPatch
def plot_poly(ax,ob):
   patch = PolygonPatch(ob, facecolor='#6699cc', edgecolor='#6699cc')
   ax.add_patch(patch)
def plot_line(ax, ob, col):
   x, y = ob.xy
   ax.plot(x, y, color=col)
def plot_multiline(ax, ob,col):
   for line in ob:
      plot_line(ax,line,col)

fig = plt.figure() 
ax = fig.gca()
plot_multiline(ax, tracks,'b')
plot_poly(ax,poly)
plot_multiline(ax, geometry,'r')

Result:

And you can save the resulting shapefile:

schema = {'geometry': 'MultiLineString','properties': {'test': 'int'}}
with fiona.open('intersect.shp', 'w', 'ESRI Shapefile', schema) as output:
     feature = {'geometry': {'coordinates': mapping(geometry)['coordinates'] ,'type': 'MultiLineString'},'properties': {'test': 1}}
     output.write(feature)

Answer 2

here's a slightly modified version of your code. my changes are in the main loop and I made it such that you iterate over all features. I did not execute it so I apologize if there are any syntax errors!

#coding:utf-8

import glob
import ogr
import json
import numpy
import matplotlib.pyplot as plt


# 1) Is there a more direct or idiomatic way to create this rectangle?
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-51.22, -30.027)
ring.AddPoint(-51.22, -30.047)
ring.AddPoint(-51.18, -30.047)
ring.AddPoint(-51.18, -30.027)
ring.AddPoint(-51.22, -30.027)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

for filename in sorted(glob.glob('../../../04 GPS/REGISTROS/*/*.gpx'), reverse=True)[:1]:

    # 2) Are all those "if not X: continue" statements the right way to check?
    datasource = ogr.Open(filename)
    if not datasource:
        continue

    layer = datasource.GetLayerByName('tracks')
    if not layer:
        continue

    feature = layer.GetNextFeature()  # Shouldn't I get *every* feature?
    while feature is not None:
        geometry = feature.GetGeometryRef().Intersection(poly)
        if geometry is not None:  #haven't run this...I assume None is returned if there's no intersection
            seg = []
            for i in range(0, geometry.GetPointCount()):
                seg.append(geom.GetPoint(i))

            seg = np.array(seg)
            plt.plot(seg[:,0], seg[:,1], lw=0.5, alpha=0.1)

        feature = layer.GetNextFeature()

plt.axis('equal')
plt.show()