在Cartopy地图上绘制shapefile城市边界

问题描述

我正在尝试使用获得的shapefile在Cartopy地形图上绘制湾区城市/城镇边界的轮廓线，

解决方案

正如 @ImportanceOfBeingErnest 和 @swatchai 所建议的，

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

解决方案

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

这篇关于在Cartopy地图上绘制shapefile城市边界的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！