Cartopy中Albers等面积投影无法使用问题求助

错误原因 #

你通过ccrs.Projection(proj_crs)创建的是Cartopy的通用投影包装类，而非针对Albers等面积投影的专用实现类（ccrs.AlbersEqualArea）。通用包装类仅提供基础投影转换能力，并未实现boundary属性——而Cartopy的GeoAxes需要该属性来绘制轴的边界，因此触发NotImplementedError。

解决方案（保留自动获取投影的逻辑） #

无需手动输入EPSG代码，只需从PyProj解析后的CRS中提取投影参数，用Cartopy对应的内置投影类实例化即可。这样既保留了从GeoTIFF自动读取投影的流程，又能获得完整的投影功能（包括边界支持）。

完整代码示例 #

from osgeo import gdal
from pyproj.crs import CRS
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# 1. 读取GeoTIFF数据及元数据（原代码不变）
data_object = gdal.Open(path_geotiff)
data_array = data_object.ReadAsArray()
data_transform = data_object.GetGeoTransform()
proj_wkt = data_object.GetProjection()

# 2. 从WKT加载CRS并提取参数
proj_crs = CRS.from_wkt(proj_wkt)
proj_params = proj_crs.to_dict()

# 3. 构建Cartopy兼容的Albers投影实例
# 配置基准面（NAD83 + GRS80椭球）
globe = ccrs.Globe(datum='NAD83', ellipse='GRS80')
# 提取Albers关键参数
proj_cartopy = ccrs.AlbersEqualArea(
    central_longitude=proj_params['lon_0'],
    central_latitude=proj_params['lat_0'],
    standard_parallels=(proj_params['lat_1'], proj_params['lat_2']),
    false_easting=proj_params['x_0'],
    false_northing=proj_params['y_0'],
    globe=globe
)

# 4. 创建GeoAxes并绘制栅格
subplot_kw = dict(projection=proj_cartopy)
fig, ax = plt.subplots(subplot_kw=subplot_kw)

# 计算栅格的地理范围
x_min = data_transform[0]
x_max = data_transform[0] + data_transform[1] * data_object.RasterXSize
y_min = data_transform[3] + data_transform[5] * data_object.RasterYSize
y_max = data_transform[3]
extent = [x_min, x_max, y_min, y_max]

# 绘制地形栅格
ax.imshow(data_array, extent=extent, origin='upper', transform=proj_cartopy)
plt.show()

关键说明 #

• 核心思路：将通用CRS解析为Cartopy内置的投影类（如AlbersEqualArea），这类内置类预先实现了boundary等GeoAxes所需的方法。
• 参数提取：通过proj_crs.to_dict()可以将PyProj的CRS转换为字典格式，直接提取投影的关键参数（标准纬线、中央经线等）。
• 基准面配置：手动指定globe参数确保与原数据的NAD83基准面匹配，避免投影偏差。

内容的提问来源于stack exchange，提问作者Fab von Bellingshausen