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 atlite
import xarray as xr
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.plot import show
from atlite.gis import shape_availability, 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")
)

/tmp/ipykernel_36247/2540265064.py:1: FutureWarning: The 'cascaded_union' attribute is deprecated, use 'unary_union' instead
  bounds = shapes.cascaded_union.buffer(1).bounds
/home/fabian/vres/py/atlite/atlite/cutout.py:185: UserWarning: Arguments module, bounds, time are ignored, since cutout is already built.
  warn(

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

<AxesSubplot: >

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

The eligibility share is: 78.99%

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)

<AxesSubplot: title={'center': 'Eligible area (green) 78.99%'}>

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

Compute availability matrix: 100%|██████████| 4/4 [00:01<00:00,  3.42 gridcells/s]

<xarray.DataArray (name: 4, y: 31, x: 45)>
array([[[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]])
Coordinates:
  * name     (name) object 'Croatia' 'Bosnia and Herz.' 'Serbia' 'Montenegro'
  * y        (y) float64 40.0 40.25 40.5 40.75 41.0 ... 46.75 47.0 47.25 47.5
  * x        (x) float64 12.75 13.0 13.25 13.5 13.75 ... 23.0 23.25 23.5 23.75

xarray.DataArray

• name: 4
• y: 31
• x: 45

• 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

  array([[[0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          ...,
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.]],
  
         [[0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          ...,
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.]],
  
         [[0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          ...,
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.]],
  
         [[0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          ...,
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.],
          [0., 0., 0., ..., 0., 0., 0.]]])

• Coordinates: (3)
  • name
    (name)
    object
    'Croatia' ... 'Montenegro'

    array(['Croatia', 'Bosnia and Herz.', 'Serbia', 'Montenegro'], dtype=object)

  • y
    (y)
    float64
    40.0 40.25 40.5 ... 47.0 47.25 47.5

    array([40.  , 40.25, 40.5 , 40.75, 41.  , 41.25, 41.5 , 41.75, 42.  , 42.25,
           42.5 , 42.75, 43.  , 43.25, 43.5 , 43.75, 44.  , 44.25, 44.5 , 44.75,
           45.  , 45.25, 45.5 , 45.75, 46.  , 46.25, 46.5 , 46.75, 47.  , 47.25,
           47.5 ])

  • x
    (x)
    float64
    12.75 13.0 13.25 ... 23.5 23.75

    array([12.75, 13.  , 13.25, 13.5 , 13.75, 14.  , 14.25, 14.5 , 14.75, 15.  ,
           15.25, 15.5 , 15.75, 16.  , 16.25, 16.5 , 16.75, 17.  , 17.25, 17.5 ,
           17.75, 18.  , 18.25, 18.5 , 18.75, 19.  , 19.25, 19.5 , 19.75, 20.  ,
           20.25, 20.5 , 20.75, 21.  , 21.25, 21.5 , 21.75, 22.  , 22.25, 22.5 ,
           22.75, 23.  , 23.25, 23.5 , 23.75])

• Indexes: (3)
  • name
    PandasIndex

    PandasIndex(Index(['Croatia', 'Bosnia and Herz.', 'Serbia', 'Montenegro'], dtype='object', name='name'))

  • y
    PandasIndex

    PandasIndex(Float64Index([ 40.0, 40.25,  40.5, 40.75,  41.0, 41.25,  41.5, 41.75,  42.0,
                  42.25,  42.5, 42.75,  43.0, 43.25,  43.5, 43.75,  44.0, 44.25,
                   44.5, 44.75,  45.0, 45.25,  45.5, 45.75,  46.0, 46.25,  46.5,
                  46.75,  47.0, 47.25,  47.5],
                 dtype='float64', name='y'))

  • x
    PandasIndex

    PandasIndex(Float64Index([12.75,  13.0, 13.25,  13.5, 13.75,  14.0, 14.25,  14.5, 14.75,
                   15.0, 15.25,  15.5, 15.75,  16.0, 16.25,  16.5, 16.75,  17.0,
                  17.25,  17.5, 17.75,  18.0, 18.25,  18.5, 18.75,  19.0, 19.25,
                   19.5, 19.75,  20.0, 20.25,  20.5, 20.75,  21.0, 21.25,  21.5,
                  21.75,  22.0, 22.25,  22.5, 22.75,  23.0, 23.25,  23.5, 23.75],
                 dtype='float64', name='x'))

• Attributes: (0)

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

<AxesSubplot: title={'center': 'name = Croatia'}, xlabel='x', ylabel='y'>

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