Note

Download this example as a Jupyter notebook here: pypsa/atlite

# Calculate Landuse Availabilities#

This example shows how `atlite`

can deal with landuse restrictions. For the demonstration the effective availability per weather cell is calculated while excluding areas specified by the CORINE CLC raster.

```
[ ]:
```

```
import geopandas as gpd
import matplotlib.pyplot as plt
import xarray as xr
import atlite
from atlite.gis import ExclusionContainer
```

With `geopandas`

we retrieve the geographic shapes for 5 countries on the Balkan Peninsula, namely Serbia, Croatia, Macedonia, Bosnia & Herzegovina and Montenegro.

```
[ ]:
```

```
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
countries = ["Serbia", "Croatia", "Bosnia and Herz.", "Montenegro"]
shapes = world[world.name.isin(countries)].set_index("name")
shapes.plot(figsize=(10, 7));
```

We create an `atlite.Cutout`

which covers the whole regions and builds the backbone for our analysis. Later, it will enable to retrieve the needed weather data.

```
[ ]:
```

```
bounds = shapes.cascaded_union.buffer(1).bounds
cutout = atlite.Cutout(
"balkans", module="era5", bounds=bounds, time=slice("2013-01-01", "2013-01-02")
)
```

Let’s see how the grid cells and the regional shapes overlap.

```
[ ]:
```

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

The CORINE Land Cover (CLC) database provides a 100 m x 100 m raster which, for each raster cell, indicates the type of landuse (forest, urban, industrial). In total there are 44 classes. Download the raster (.tif file) from the download page and store the raster as `corine.tif`

.

For calculating the availability per cutout weather cells, an `ExclusionContainer`

must be defined beforehand. It serves as a container for all rasters and geometries we want to exclude (or possibly include).

In many cases, rasters and geometries have different projections and resolutions. Therefore, the `ExclusionContainer`

is initialized by a CRS and a resolution which suits as a basis for all added rasters and geometries. Per default the CRS is 3035 and the resoultion 100, which leads set a raster of 100 meter resolution. All rasters and geometries will be converted to this `(crs, res)`

config if they don’t match it.

When adding a raster to the `ExclusionContainer`

you can specify which codes (integers) to exclude. By setting `invert=True`

, you can also restrict the **inclusion** to a set of codes. Further you can buffer around codes (see the docs for detail). Here we are going to exclude the first twenty landuse codes.

```
[ ]:
```

```
CORINE = "corine.tif"
excluder = ExclusionContainer()
excluder.add_raster(CORINE, codes=range(20))
```

For the demonstration we want to see how the landuse availability behaves within one specific shape, e.g. Croatia.

Note that, since the excluder is in `crs=3035`

, we convert to geometry of Croatia to `excluder.crs`

for plotting it…

```
[ ]:
```

```
croatia = shapes.loc[["Croatia"]].geometry.to_crs(excluder.crs)
```

…and use the `shape_availability`

function of `atlite`

to calculate a mask for the `ExclusionContainer`

excluder.

```
[ ]:
```

```
masked, transform = excluder.compute_shape_availability(croatia)
```

The `masked`

object is a numpy array. Eligile raster cells have a 1 and excluded cells a 0. Note that this data still lives in the projection of `excluder`

. For calculating the eligible share we can use the following routine.

```
[ ]:
```

```
eligible_share = masked.sum() * excluder.res**2 / croatia.geometry.item().area
print(f"The eligibility share is: {eligible_share:.2%}")
```

For plotting the geometry together with the excluded areas, we can use the function `plot_shape_availability`

which uses `rasterio`

’s and `geopandas`

’ plot function in the background.

```
[ ]:
```

```
fig, ax = plt.subplots()
excluder.plot_shape_availability(croatia)
```

How does is look when we add our cutout grid to the plot? How do the weather cells intersect with the available area?

```
[ ]:
```

```
fig, ax = plt.subplots()
excluder.plot_shape_availability(croatia, ax=ax)
cutout.grid.to_crs(excluder.crs).plot(edgecolor="grey", color="None", ax=ax, ls=":")
ax.set_title(f"Eligible area (green) {eligible_share:.2%}")
ax.axis("off");
```

We see that the weather cells are much larger than the raster cells. `GDAL`

provides a fast reprojection function for averaging fine-grained to coarse-grained rasters. `Atlite`

automates this calculation for all geometries in `shapes`

when calling the `cutout.availabilitymatrix`

function. Let’s see how this function performs. (Note that the steps before are not necessary for this calculation.)

**INFO**: *For large sets of shapes set nprocesses to a number > 1 for parallelization.*

```
[ ]:
```

```
A = cutout.availabilitymatrix(shapes, excluder)
A
```

`A`

is an DataArray with 3 dimensions (`shape`

, `x`

, `y`

) and very sparse data. It indicates the relative overlap of weather cell `(x, y)`

with geometry `shape`

while excluding the area specified by the `excluder`

.

How does the availability look for our example of Croatia?

```
[ ]:
```

```
fig, ax = plt.subplots()
A.sel(name="Croatia").plot(cmap="Greens")
shapes.loc[["Croatia"]].plot(ax=ax, edgecolor="k", color="None")
cutout.grid.plot(ax=ax, color="None", edgecolor="grey", ls=":")
```

Note that now the projection is in `cutout.crs`

. In the north-west, where most of the areas were excluded, the availability is lower than 0.5. That means less than the half of these weather cells and their potentials can be exploited.

And for the other shapes…

```
[ ]:
```

```
fg = A.plot(row="name", col_wrap=3, cmap="Greens")
fg.set_titles("{value}")
for i, c in enumerate(shapes.index):
shapes.loc[[c]].plot(ax=fg.axs.flatten()[i], edgecolor="k", color="None")
```

The availibility matrix `A`

can now be used as a layoutmatrix in the conversion functions of `atlite`

, i.e. `cutout.pv`

, `cutout.wind`

. The normal approach would be to weigh the availabilities with the area per grid cell and the capacity per area.

```
[ ]:
```

```
cap_per_sqkm = 1.7
area = cutout.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6
area = xr.DataArray(area, dims=("spatial"))
capacity_matrix = A.stack(spatial=["y", "x"]) * area * cap_per_sqkm
```

After the cutout preparation, we can calculate the static and dynamic capacity factors of each region.

```
[ ]:
```

```
cutout.prepare()
pv = cutout.pv(
matrix=capacity_matrix,
panel=atlite.solarpanels.CdTe,
orientation="latitude_optimal",
index=shapes.index,
)
```

Finally let’s see how the total power potential per region look.

```
[ ]:
```

```
pv.to_pandas().div(1e3).plot(ylabel="Solar Power [GW]", ls="--", figsize=(10, 4));
```