name = '2015-08-03-fiona_gpx'
title = "Reading GPX files directly with fiona"
%matplotlib inline
import os
from datetime import datetime
from IPython.core.display import HTML
with open('creative_commons.txt', 'r') as f:
html = f.read()
hour = datetime.utcnow().strftime('%H:%M')
comments="true"
date = '-'.join(name.split('-')[:3])
slug = '-'.join(name.split('-')[3:])
metadata = dict(title=title,
date=date,
hour=hour,
comments=comments,
slug=slug,
name=name)
markdown = """Title: {title}
date: {date} {hour}
comments: {comments}
slug: {slug}
{{% notebook {name}.ipynb cells[2:] %}}
""".format(**metadata)
content = os.path.abspath(os.path.join(os.getcwd(),
os.pardir,
os.pardir,
'{}.md'.format(name)))
with open('{}'.format(content), 'w') as f:
f.writelines(markdown)
html = '''
<small>
<p> This post was written as an IPython notebook.
It is available for <a href='https://ocefpaf.github.com/python4oceanographers/downloads/notebooks/%s.ipynb'>download</a>
or as a static <a href='https://nbviewer.ipython.org/url/ocefpaf.github.com/python4oceanographers/downloads/notebooks/%s.ipynb'>html</a>.</p>
<p></p>
%s''' % (name, name, html)
In a previous post I said I would try to build a GPX loader for GeoPandas. If you try to load a GPX file right now this is what you get.
import geopandas as gpd
fname = './data/2014_08_05_farol.gpx'
gdf = gpd.read_file(fname)
gdf
import fiona
fiona.listlayers(fname)
[u'waypoints', u'routes', u'tracks', u'route_points', u'track_points']
Let's get that track layer into fiona.
layer = fiona.open(fname, layer='tracks')
layer
<open Collection './data/2014_08_05_farol.gpx:tracks', mode 'r' at 0x7f17e450e850>
The layer object has a lot of properties/methods. Some of the properties carries the underlying metadata. Here are some important ones:
layer.crs, layer.bounds
({'init': u'epsg:4326'}, (-38.53207166666667, -13.01158, -38.49701, -13.001495))
How many items we have in this object?
len(list(layer.items()))
1
What are these items?
geom = layer[0]
type(geom)
dict
geom.keys()
['geometry', 'type', 'id', 'properties']
geom['type'], geom['id'], geom['properties']
('Feature', '0', OrderedDict([(u'name', None), (u'cmt', None), (u'desc', None), (u'src', None), (u'link1_href', None), (u'link1_text', None), (u'link1_type', None), (u'link2_href', None), (u'link2_text', None), (u'link2_type', None), (u'number', None), (u'type', None)]))
OK. It is a dictionary with more data and other dictionaries in it.
We are interested in the data stored into geometry['coordinates']
.
Instead of dumping a bunch of numbers in the screen lets create a
Shapely
object to inspect the geometry coordinates.
from shapely.geometry import shape
data = {'type': 'MultiLineString',
'coordinates': geom['geometry']['coordinates']}
shp = shape(data)
shp
Nice! A simple and quick way to extract tracks from GPX and export to Shapely. But this image makes no sense without a map ;-)
(Note that I created a new dictionary with a structure I know relates to a GeoJson before feeding it to Shapely.)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io import shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
def make_map(projection=ccrs.PlateCarree()):
fig, ax = plt.subplots(figsize=(9, 13),
subplot_kw=dict(projection=projection))
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
import cartopy.io.img_tiles as cimgt
request = cimgt.OSM()
extent = [-38.54, -38.48,
-13.02, -12.98]
fig, ax = make_map(projection=request.crs)
ax.set_extent(extent)
img = ax.add_image(request, 14)
s = ax.add_geometries(shp, ccrs.PlateCarree(),
facecolor='none',
edgecolor='crimson',
linewidth=2)
Have fun exploring your GPX files!
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.