使用Atlite创建巴伐利亚多面行政区Cutout遇索引错误求助

问题根源 #

你遇到的IndexError是因为Cutout未生成有效网格，核心原因是shapefile坐标系不匹配：巴伐利亚官方shapefile使用投影坐标系（如EPSG:32632/32633），而Atlite要求bounds必须是WGS84（EPSG:4326）的经纬度范围（西、南、东、北）。此外，旧版GeoPandas的cascaded_union已被弃用，可能导致边界计算异常。

修复步骤 #

1. 转换坐标系并修正边界计算 #

替换边界计算代码，先将shapefile转为WGS84，再用稳定的方法计算整体边界：

counties = ["Bezirksverwaltung Oberbayern", "Bezirksverwaltung Niederbayern", "Bezirksverwaltung Oberpfalz",
            "Bezirksverwaltung Oberfranken", "Bezirksverwaltung Mittelfranken", "Bezirksverwaltung Unterfranken",
            "Bezirksverwaltung Schwaben"]

# 读取shapefile并转换到WGS84经纬度坐标系
bavaria = gpd.read_file(input_path + '\Verwaltungseinheit.shp').to_crs(epsg=4326)
bavaria_shapes = bavaria[bavaria.name.isin(counties)].explode(index_parts=True).set_index("name")

# 用unary_union替代已弃用的cascaded_union计算整体区域
bavaria_union = bavaria_shapes.unary_union
# 缓冲1度（WGS84下单位为经纬度，按需调整）
bavaria_bounds = bavaria_union.buffer(1).bounds

2. 验证边界格式 #

确保bavaria_bounds输出为(west, south, east, north)的经纬度数值，示例：

print("Cutout Bounds:", bavaria_bounds)
# 正确输出类似：(8.01, 47.27, 13.99, 50.52)

3. 重新创建并验证Cutout #

# 初始化Cutout
cutout = atlite.Cutout(
    path="bavaria_cutout_2017.nc",
    module="era5",
    bounds=bavaria_bounds,
    time="2017",
)

# 先执行prepare()完成数据加载
cutout.prepare()

# 绘制验证网格
plt.rc("figure", figsize=[10, 7])
fig, ax = plt.subplots()
bavaria_shapes.plot(ax=ax)
cutout.grid.plot(ax=ax, edgecolor="grey", color="None")
plt.show()

关键注意事项 #

• 坐标系转换是核心：Atlite依赖ERA5的WGS84经纬度网格，必须确保输入边界为经纬度格式。
• 按公里缓冲的处理：若需要按实际距离缓冲，需先转投影坐标系缓冲，再转回WGS84：

  # 示例：缓冲10公里
  bavaria_proj = bavaria_shapes.to_crs(epsg=32632)  # UTM 32N适配巴伐利亚
  union_proj = bavaria_proj.unary_union.buffer(10000)  # 缓冲10000米
  union_wgs84 = gpd.GeoSeries([union_proj], crs=epsg=32632).to_crs(epsg=4326).iloc[0]
  bavaria_bounds = union_wgs84.bounds
  

内容的提问来源于stack exchange，提问作者m325