Участок шейп-файла с границами города на вершине картографической карты

Question

Я пытаюсь нанести контур границ города / района залива на вершине картографической карты местности, используя шейп-файл, полученный здесь и следуя этому примеру .По какой-то причине границы не отображаются, даже когда я указываю границы сверху через zorder.Я что-то упустил?

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

Answer 1

Как подсказывают @ImportanceOfBeingErnest и @swatchai, параметр CRS (система координат) в ShapelyFeature cartopy.feature.ShapelyFeature() был неверным.

Правильный код EPSG (European Petroleum Survey Group?) Можно найти в одном из XML-файлов, включенных в шейп-файл:

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

и передача его в качестве второго параметра в ShapelyFeature() - это все, что нужно, чтобы получить шейп-файл для правильного построения границ города:

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