为具有多边形的地理形状创建 atlite 切口

Question

我想使用 Atlite 计算巴伐利亚的光伏和风能可用性，得出所有 7 个县的结果，但我无法创建切口。对于数据基础，我使用此处的 shapefile https://geodaten.bayern.de/opengeodata/OpenDataDetail.html?pn=atkis_basis_dlm，其中包括巴伐利亚的所有行政区域。从那里我选择了想要的区域（=县），并按照 atlite 示例了解土地利用可用性（https://atlite.readthedocs.io/en/latest/examples/landuse-availability.html#Calculate-Landuse-Availability）因为它还包含多个切口区域。

counties = ["Bezirksverwaltung Oberbayern", "Bezirksverwaltung Niederbayern", "Bezirksverwaltung Oberpfalz",
                "Bezirksverwaltung Oberfranken", "Bezirksverwaltung Mittelfranken", "Bezirksverwaltung Unterfranken",
                "Bezirksverwaltung Schwaben"]
    bavaria = gpd.read_file(input_path+'\Verwaltungseinheit.shp')
    bavaria_shapes = bavaria[bavaria.name.isin(counties)].explode(index_parts=True).set_index("name")

绘图效果很好

bavaria_shapes.plot(figsize=(10, 7))
    plt.show()

然后要创建切口，我继续这样

bavaria_bounds = bavaria_shapes.cascaded_union.buffer(1).bounds
cutout = atlite.Cutout(
        path="bavaria_cutout_2017.nc",  # path to where the cutout will be stored
        module="era5",
        bounds=bavaria_bounds,
        time="2017",
    )

进一步按照 Atlite 文档的示例，我想用覆盖网格绘制切口。

plt.rc("figure", figsize=[10, 7])
fig, ax = plt.subplots()
bavaria_shapes.plot(ax=ax)
cutout.grid.plot(ax=ax, edgecolor="grey", color="None")

这里我得到以下错误：

INFO:atlite.cutout:Building new cutout bavaria_cutout_2017.nc
Traceback (most recent call last):
  File "C:\Users\xxx\main.py", line 219, in <module>
    cutout.grid.plot(ax=ax, edgecolor="grey", color="None")
    ^^^^^^^^^^^
  File "C:\Users\xxx\Lib\site-packages\atlite\utils.py", line 171, in __get__
    result = self.method(inst)
             ^^^^^^^^^^^^^^^^^
  File "C:\Users\xxx\Lib\site-packages\atlite\cutout.py", line 402, in grid
    span = (coords[self.shape[1] + 1] - coords[0]) / 2
            ~~~~~~^^^^^^^^^^^^^^^^^^^
IndexError: index 1 is out of bounds for axis 0 with size 0

我很高兴了解输入数据与示例中的数据有何不同，或者如何更改我的形状文件以创建有效的剪切图。

为了找出巴伐利亚形状和示例形状的不同之处，我研究了数据类型，发现巴伐利亚形状不仅由 7 个多边形组成，而且其中一些实际上是多边形。我解决这个问题的方法是在 bavaria_shapes 上使用 geopandas.GeoDataFrame.explode() 来分割多边形。错误与之前相同。 我也尝试过，如果这只发生在 cutout.grid 函数中，但它发生在任何尝试访问剪切的函数中（接下来我尝试的是 cutout.prepare() ）。

Answer 1

答案是改变坐标参考系。 ERA5 使用 EPSG 4326，因此通过更改输入 CRS 解决了该问题。 为了重新投影，我使用了 .to_crs