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()
Best Answer
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:
With shapely and the wkt format:
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)
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 ?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)
Plot the feature
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 formatResult:
It is easy to adapt these solutions with:
New: solution with Fiona and shapely
Polygon
Intersection:
Plot:
Result:
And you can save the resulting shapefile: