Plot shapefile city borders on top of cartopy map

2019-08-27 20:42发布

问题:

I'm trying to plot the outline of Bay Area city/town borders on top of a cartopy terrain map using a shapefile obtained here and following this example. For some reason, the borders don't show up, even when I specify borders are on top via zorder. Am I missing something?

# import functions
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

# Create a Stamen terrain background instance
stamen_terrain = cimgt.Stamen('terrain-background')
fig = plt.figure(figsize = (10, 10))
ax = fig.add_subplot(1, 1, 1, projection=stamen_terrain.crs)

# Set range of map, stipulate zoom level
ax.set_extent([-122.7, -121.5, 37.15, 38.15], crs=ccrs.Geodetic())
ax.add_image(stamen_terrain, 12, zorder = 0)

# Add city borders - not working
filename = r'./shapefile/ba_cities.shp' # from https://earthworks.stanford.edu/catalog/stanford-vj593xs7263
shape_feature = ShapelyFeature(Reader(filename).geometries(), ccrs.PlateCarree(), edgecolor='black')
ax.add_feature(shape_feature, zorder = 1)
plt.show()

回答1:

As @ImportanceOfBeingErnest and @swatchai suggest, the CRS (coordinate reference system) parameter in ShapelyFeature cartopy.feature.ShapelyFeature() was incorrect.

The proper EPSG (European Petroleum Survey Group?) code can be found in one of the .xml files included with the shapefile:

   <gco:CharacterString>26910</gco:CharacterString>
</code>
<codeSpace>
   <gco:CharacterString>EPSG</gco:CharacterString>

and passing this as the second parameter in ShapelyFeature() is all it takes to get the shapefile to plot the city borders properly:

# Add city borders
filename = r'./shapefile/ba_cities.shp'
shape_feature = ShapelyFeature(Reader(filename).geometries(), ccrs.epsg(26910), 
                               linewidth = 1, facecolor = (1, 1, 1, 0), 
                               edgecolor = (0.5, 0.5, 0.5, 1))
ax.add_feature(shape_feature)
plt.show()