0

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()
4

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:
enter image description here

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:

enter image description here

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:

enter image description here

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)
  • Thank you very much for your nice answer! Perhaps I have not expressed myself right, but I have NO special interest in GeoJson format. Actually, I am fond of Numpy, and usually I convert things to numpy arrays. Also, I use OGR because I could not use Fiona to read GPX files (perhaps I did something wrong). Your suggestion of using Fiona + Shapely is very interesting to me, would you mind posting some example code? – heltonbiker Oct 3 '14 at 13:22
  • look above for a solution with Fiona-shapely – gene Oct 4 '14 at 17:55
  • Thank you a lot. I haven't tested the answer yet, but I learned many new ways of doing things here, and it helped me very much. – heltonbiker Oct 4 '14 at 18:31
1

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()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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