Docker>> basemaps>> 返回
项目作者： charlie1347

项目描述 ：
Playing with basemaps for geographic data visualisations
高级语言： Jupyter Notebook
项目地址: git://github.com/charlie1347/basemaps.git
创建时间： 2017-03-22T18:28:56Z
项目社区：https://github.com/charlie1347/basemaps
开源协议：
下载


Basemaps in matplotlib/bokeh

Matplotlib’s “basemap” ability allows the option of using shapefiles to add basemaps to your plots. When experimenting with this function, I ran into the following error:

  1. ValueError: shapefile must have lat/lon vertices  - it looks like this one has vertices
  2. in map projection coordinates. You can convert the shapefile to geographic
  3. coordinates using the shpproj utility from the shapelib tools
  4. (http://shapelib.maptools.org/shapelib-tools.html)

Some googling turned up a post stating that GDAL’s ogr2ogr could handle these conversions, as opposed to figuring out how to use the shapelib tools.

The first requirement is to set up a working version of GDAL on your computer - I compiled this from source in Linux. Once I had GDAL working, I used ogr2ogr to convert my shapefile. ogr2ogr contains documentation on how to run it from the command line here. Alternatively, I experimented with the ogr2ogr Python module, taken from here. The Python module takes roughly the same syntax as the command line ogr2ogr, with individual arguments passed as a list.

The first shapefile I was working with was a shapefile of all the buildings in London, TQTL_LON.shp, downloaded from here. Looking at the original TQTL_LON.prj file, the shapefile is using the British National Grid projection - ogr2ogr is smart enough to figure this out on its own. ogr2ogr can also clip your plots to just a specific area by feeding it a bounding box of coordinates. 

  1. ogr2ogr.main(["", "-f", "ESRI Shapefile", 
  2.                 "-clipdst", "-0.237947", "51.449814", "0.007965", "51.553407",
  3.                 "Converted basemaps/London_buildings.shp", 
  4.                 "Raw basemaps/TQTL_LON.shp", 
  5.                 "-t_srs", "EPSG:4326"])

This code takes the original “TQTL_LON.shp” shapefile, clips it to a bounding box, and reprojects it into a new ESRI Shapefile called “London_buildings.shp”, using the coordinate system EPSG:4326 (also known as WGS84).

I took a second shapefile, TQ_Roadlink.shp, from the Ordnance Survey website here, this time of all of the UK’s roads. Running the above code on this shapefile however produced the following error:

  1. ValueError: readshapefile can only handle 2D shape types

A simple fix - ogr2ogr can deal with this by forcing the coordinate dimensions (from XYZ to just XY). 

  1. ogr2ogr.main(["", "-f", "ESRI Shapefile",
  2.                 "-clipdst", "-0.237947", "51.449814", "0.007965", "51.553407",
  3.                 "Converted basemaps/London_roads.shp",
  4.                 "Raw basemaps/TQ_RoadLink.shp",
  5.                 "-t_srs", "EPSG:4326",
  6.                 "-dim", "2"])

Once I had converted both shapefiles, matplotlib could then easily plot both as basemaps (see ipynb for full code):

Buildings

Matplotlib buildings

Roads

Matplotlib roads

Bokeh

To add basemaps onto bokeh plots, I needed to convert the shapefiles into geoJSON format. Again, ogr2ogr can handle this:

  1. ogr2ogr.main(["", "-f", "GeoJSON", 
  2.                 "-clipdst", "-0.237947", "51.449814", "0.007965", "51.553407",
  3.                 "Converted basemaps/London_buildings.geojson", 
  4.                 "Raw basemaps/TQTL_LON.shp", 
  5.                 "-t_srs", "EPSG:4326"])

Once the shapefiles are in a geoJSON format, I used bokeh’s GeoJSONDataSource to add them both as “patches” to my bokeh plots (again, see ipynb file for full code):

Bokeh