Plot shapefile city borders on top of cartopy map

Question

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

It doesn't look like the shape is even added correctly, meaning there is no geometry which intersects with the extent of the axes. — ImportanceOfBeingErnest, Jun 30 '19 at 20:08
Could you elaborate on what it means to add the shape correctly? Not sure what change you're recommending. — Michael Boles, Jun 30 '19 at 20:41
I do not recommend anything; I just wanted to point out that the shape is not even drawn. So one would need to find the reason for that, not look into zorder. — ImportanceOfBeingErnest, Jun 30 '19 at 20:45
I used [some files](https://www.naturalearthdata.com/downloads/110m-physical-vectors/) from naturalearth.com and with those the code runs fine. So potentially the file you are using is simply corrupted or is at least incompatible with cartopy. — ImportanceOfBeingErnest, Jun 30 '19 at 21:47
Based on the metadata of the shapefile, the CRS code is EPSG:26910 or NAD83 / UTM zone 10N. Using `ccrs.PlateCarree()` in ShapelyFeature() is not appropriate. — swatchai, Jul 01 '19 at 12:23

score 4 · Accepted Answer · answered Jul 02 '19 at 22:14

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

Plot shapefile city borders on top of cartopy map

1 Answers1

Linked