coordinatesfionamappingpythonshapely
I have written a parser in Python that extracts the coordinates out of a GPX file and I would like to view those points on a map. I have done a little research and it seems like either Fiona or Shapely is capable of doing this, but all of the examples I find handle shape files.
My coordinates are stored as tuples in a list. I was wondering if there was an easy way to use either of these packages to plot my coordinates.
I had the same issue when packaging a fiona project with py2exe. I think py2exe simply cannot resolve all the fiona dependencies. But fortunately fiona is well behaved: just add it under packages into the py2exe script and it will probably also work for you:
opts = {
'py2exe': {
'packages':['fiona'],
'bundle_files': 3
}
}
setup(
console=[{"script" : "test.py"}],
options=opts,
zipfile=None
)
As Fiona works with dictionaries, it is easy.
With the example of Finding the common borders between polygons in the same shapefile
import fiona
from shapely.geometry import shape, mapping
# load the shapefile as an iterable
polys = fiona.open("polygons.shp")
Now, with the itertools module:
import itertools
for pol1,pol2 in itertools.combinations(polys, 2):
if shape(pol1['geometry']).touches(shape(pol2['geometry'])):
print "Polygon ",pol1['properties']['id'], " touches Polygon " , pol2['properties']['id']
Polygon 1 touches Polygon 4
Polygon 2 touches Polygon 4
Polygon 3 touches Polygon 4
And directly:
schema = {'geometry': 'LineString','properties': OrderedDict([(u'poly1', 'int:10'), (u'poly2', 'int:10')])}
with fiona.open('intersection2.shp','w','ESRI Shapefile', schema) as output:
with fiona.open("intersections.shp") as polys:
for pol1,pol2 in itertools.combinations(polys, 2):
geompol1 = shape(pol1['geometry'])
geompol2 = shape(pol2['geometry'])
if geompol1.touches(geompol2):
output.write({'geometry':mapping(geompol1.intersection(geompol2)), 'properties': {'poly1':pol1['properties']['id'],'poly2':pol2['properties']['id']}})
Control
for line in fiona.open("intersection2.shp"):
print line['properties']
OrderedDict([(u'poly1', 1), (u'poly2', 4)])
OrderedDict([(u'poly1', 2), (u'poly2', 4)])
OrderedDict([(u'poly1', 3), (u'poly2', 4)])
Best Answer
Given you're already familiar with Python I'd recommend taking a look at the
cartopylibrary for plotting spatial data with matplotlib. It's a much nicer replacement to matplotlib-basemap.The basic idea of cartopy is that you can apply a projection to both your matplotlib frame as well as everything you're plotting to the frame. It also includes a lot of helper methods for adding base data, displaying shapely geometries, etc (even in v0.12 released earlier this month, reading in mapbox tiles as a basemap!).
In your case you've got a set of points that are (presumably) in WGS84 longitude, latitude pairs. You could build a shapely geometry and plot that, but the easiest method would be to use the standard matplotlib scatter to see the data. In which case I'm going to assume that your points rather than being a list of tuples, are in a 2d numpy array called
points. For my example I grabbed a trail from the OpenStreetMap gpx trail store.First you should define some projections. See the cartopy projection list for a whole stack of different options.
Then build the axis and plot some points with a bit of basemap information.
Finally add your data (and apply the transform)