CODEWITHSUNDEEP
pythongeopandasshapefilepypsa
m325
m325

Reputation: 21

Create atlite cutout for geographic shapes with multipolygons

I want to use Atlite to compute pv and wind availabilities for bavaria in a resultion that yields all 7 counties, but I cannot create a cutout. For the data basis, I use the shapefile from here https://geodaten.bayern.de/opengeodata/OpenDataDetail.html?pn=atkis_basis_dlm which includes all administrative regions in bavaria. From there I selected the wanted regions (=counties) and followed the atlite example for landuse availabilities (https://atlite.readthedocs.io/en/latest/examples/landuse-availability.html#Calculate-Landuse-Availabilities) as it also conists of more than one area for the cutout.

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

The plotting works just fine with

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

then the cutout is to be created and I continued like this

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",
    )

following the example of the Atlite documentation further, I want to plot the cutout with an overlaying grid.

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

Here I get the following error:

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

I'd be happy for any idea on how the input data differs from the one in the example or on how I can change my shapefile to create a valid cutout.

To find in what way the Bavaria shape and the example shape differ I looked into the datatypes and found, that the bavaria shape doesn't only consits of 7 polygons, but that some of them are in fact multipolygons. My approach to solve this was to use geopandas.GeoDataFrame.explode() on bavaria_shapes to split the multipolygons. The error is the same as before. I also tried if this only occured with the cutout.grid function, but it occurs in any function that tries to access the cutout (next I tried was cutout.prepare() ).

Upvotes: 1

Views: 117

Answers (1)

m325
m325

Reputation: 21

The answer was to change the coordinate reference system. The ERA5 uses EPSG 4326 and therefore with changing the input CRS solved the problem. For reprojecting I used .to_crs

Upvotes: 1

Related Questions