python - 使用 ccrs.epsg() 绘制带有 EPSG 4326 坐标系的邮政编码周长 shapefile

Question

我从这里获得了一个邮政编码周长的 shapefile，并希望将它们绘制在 Cartopy 地图的顶部，就像我在这个例子中所做的那样。

根据消息来源，该数据位于 EPSG 4326 坐标系中。当我尝试绘制数据时

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt 
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 = (mapsize,mapsize))
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, 8, zorder = 0)

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.epsg(4326), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

我看到以下错误：

ValueError: EPSG code does not define a projection

可能相关——当我查看ccrs.epsg()函数时，它说这个 EPSG 代码不受支持

help(ccrs.epsg)
Help on function epsg in module cartopy.crs:

epsg(code)
    Return the projection which corresponds to the given EPSG code.

    The EPSG code must correspond to a "projected coordinate system",
    so EPSG codes such as 4326 (WGS-84) which define a "geodetic coordinate
    system" will not work.

    Note
    ----
        The conversion is performed by querying https://epsg.io/ so a
        live internet connection is required.

鉴于此结果，我还尝试使用ccrs.Geodetic()：

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.Geodetic(), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

但这也无法显示邮政编码周长，并显示以下警告消息：

UserWarning: Approximating coordinate system <cartopy._crs.Geodetic object at 0x1a2d2375c8> with the PlateCarree projection.
  'PlateCarree projection.'.format(crs))

我也试过ccrs.PlateCarree()，但没有运气。请帮忙！

标签： pythonmatplotlibshapefileshapelycartopy

Answer 1

要将来自不同来源的数据绘制在一起，必须coordinate reference system为每个数据集声明正确。对于 shapefile，您可以在其随附xxx.prj文件中找到它。

这是工作代码：

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

shapefile_name= "./data/ZIPCODE.shp"
mapwidth, mapheight = 8, 8
pad = 0.25

stamen_terrain = cimgt.Stamen('terrain-background')
stm_crs = stamen_terrain.crs

fig = plt.figure(figsize = (mapwidth, mapheight))
ax = fig.add_subplot(1, 1, 1, projection=stm_crs)  #world mercator

# Set extent of map
ax.set_extent([-123.3-pad, -121.5+pad, 37.05-pad, 38.75+pad], crs=ccrs.Geodetic())
# Plot base map
ax.add_image(stamen_terrain, 8, zorder=0)

# Add polygons from shapefile
# Note: the use of `ccrs.epsg(26910)`
shape_feature = ShapelyFeature(Reader(shapefile_name).geometries(), ccrs.epsg(26910))

# You can choose one of the 2 possible methods to plot
# ... the geometries from shapefile
# Styling is done here.
method = 1
if method==1:
    # iteration is hidden
    ax.add_feature(shape_feature, facecolor='b', edgecolor='red', alpha=0.4, zorder = 15)
if method==2:
    # iterate and use `.add_geometries()`
    # more flexible to manipulate particular items
    for geom in shape_feature.geometries():
        ax.add_geometries([geom], crs=shape_feature.crs, facecolor='b', edgecolor='red', alpha=0.4)

plt.show()

输出图：