Solutions rasterio & atlite

Exercise 1: Land Exclusion Sensitivities

Tweak some parameters of the land exclusion calculation for wind deployment and observe the changes in available area and resulting wind power potential. For instance:

1. Vary the CORINE land classes,

2. Vary the buffer distance, for instance, to 0m or 2000m,

3. Add the GEBCO topography data loaded at the beginning of the tutorial to the exclusion calculator and exclude areas above 1000m elevation.

Solution:

1. Change numbers in codes_to_include or codes_to_exclude.

2. Change buffer keyword argument in excluder.add_raster(..., buffer=2000).

3. Add line excluder.add_raster(fn_to_gebco, crs=4326). The crs=4326 is essential as data is given in another different projection.

Exercise 2: Solar PV Orientation Sensitivities

The orientation of solar panels plays a decisive role in determining their potential power generation time series (as capacity factors using the per_unit=True keyword argument). Explore this aspect by performing the following sensitivity analysis on the solar PV setup for Portugal in January:

Task 1: Plot the total energy yield in January as a function of the slope of the panels oriented towards the South (0° means horizontal orientation, 90° means vertical orientation) in steps of 10°. What is the optimal tilt angle? Why might the value be higher than typically recommended for PV panels?

Task 2: Plot the power production profile for January 1 for the azimuth angles 0°, 90°, 180°, 270° of panels tilted at a 30° angle (i.e. varying north, east, south west orientations). What orientation does each angle represent? Describe the differences you observe.

Task 3: The cutout.pv() function can also consider multiple options for tracking features of solar panels as listed below. Tracking in solar panel systems refers to the mechanism that allows panels to move and orient themselves towards the sun, maximizing their exposure to sunlight and thereby increasing their electricity output. The default settings in atlite assume no tracking. Explore and describe the power output profile for different tracking options for January 31. How much more expensive in percent can the equipment for tracking be to achieve the same cost per unit of electricity produced as a fixed panel (if it only operated in January)?

Available tracking options in atlite are:

• tracking='horizontal' for 1-axis horizontal tracking: This setting positions solar panels to rotate around a horizontal axis, allowing them to follow the sun’s path from east to west.

• tracking='vertical' for 1-axis vertical tracking: In this mode, the solar panels rotate around a vertical axis, which can be advantageous in high latitudes for better sun capture throughout the day.

• tracking='dual' for 2-axis tracking: This advanced system allows panels to track the sun both horizontally and vertically, maximizing solar exposure and energy capture throughout the day and year.

For all calculations, disregard land eligibility criteria and assume a uniform distribution of solar panels across Portugal.

import atlite
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from urllib.request import urlretrieve
from os.path import basename, exists

url = (
    "https://tubcloud.tu-berlin.de/s/2oogpgBfM5n4ssZ/download/PORTUGAL-2013-01-era5.nc"
)
if not exists(basename(url)):
    urlretrieve(url, basename(url))

cutout = atlite.Cutout("PORTUGAL-2013-01-era5.nc")

url = "https://tubcloud.tu-berlin.de/s/PisSQdRzwbWmzPt/download/country_shapes_reduced.geojson"
countries = gpd.read_file(url).set_index("name")
shapes = countries.loc[["PT"]].geometry

cutout.pv(
    panel="CSi",
    orientation={"slope": 30.0, "azimuth": 180.0},
    shapes=shapes,
    per_unit=True,
)

<xarray.DataArray (time: 768, name: 1)> Size: 6kB
array([[0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [4.23259202e-02],
       [1.54659473e-01],
       [2.96705202e-01],
       [3.54283778e-01],
       [4.23966612e-01],
       [4.00369283e-01],
       [3.08207520e-01],
       [1.92481709e-01],
       [8.77487256e-02],
       [0.00000000e+00],
       [0.00000000e+00],
...
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [4.97159643e-02],
       [1.44374187e-01],
       [2.51331753e-01],
       [3.70007494e-01],
       [4.49929139e-01],
       [4.70099805e-01],
       [4.42353101e-01],
       [3.43244557e-01],
       [1.89051756e-01],
       [5.64463039e-02],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00],
       [0.00000000e+00]])
Coordinates:
  * time     (time) datetime64[ns] 6kB 2013-01-01 ... 2013-02-01T23:00:00
  * name     (name) object 8B 'PT'
Attributes:
    units:    p.u.

xarray.DataArray

• time: 768
• name: 1

• 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.1891 0.05645 0.0 0.0 0.0 0.0 0.0

  array([[0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [4.23259202e-02],
         [1.54659473e-01],
         [2.96705202e-01],
         [3.54283778e-01],
         [4.23966612e-01],
         [4.00369283e-01],
         [3.08207520e-01],
         [1.92481709e-01],
         [8.77487256e-02],
         [0.00000000e+00],
         [0.00000000e+00],
  ...
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [4.97159643e-02],
         [1.44374187e-01],
         [2.51331753e-01],
         [3.70007494e-01],
         [4.49929139e-01],
         [4.70099805e-01],
         [4.42353101e-01],
         [3.43244557e-01],
         [1.89051756e-01],
         [5.64463039e-02],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00],
         [0.00000000e+00]])

• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    2013-01-01 ... 2013-02-01T23:00:00

    array(['2013-01-01T00:00:00.000000000', '2013-01-01T01:00:00.000000000',
           '2013-01-01T02:00:00.000000000', ..., '2013-02-01T21:00:00.000000000',
           '2013-02-01T22:00:00.000000000', '2013-02-01T23:00:00.000000000'],
          shape=(768,), dtype='datetime64[ns]')

  • name
    (name)
    object
    'PT'

    array(['PT'], dtype=object)

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2013-01-01 00:00:00', '2013-01-01 01:00:00',
                   '2013-01-01 02:00:00', '2013-01-01 03:00:00',
                   '2013-01-01 04:00:00', '2013-01-01 05:00:00',
                   '2013-01-01 06:00:00', '2013-01-01 07:00:00',
                   '2013-01-01 08:00:00', '2013-01-01 09:00:00',
                   ...
                   '2013-02-01 14:00:00', '2013-02-01 15:00:00',
                   '2013-02-01 16:00:00', '2013-02-01 17:00:00',
                   '2013-02-01 18:00:00', '2013-02-01 19:00:00',
                   '2013-02-01 20:00:00', '2013-02-01 21:00:00',
                   '2013-02-01 22:00:00', '2013-02-01 23:00:00'],
                  dtype='datetime64[ns]', name='time', length=768, freq=None))

  • name
    PandasIndex

    PandasIndex(Index(['PT'], dtype='object', name='name'))

• Attributes: (1)

  units :

  p.u.

Task 1 Solution:

Show code cell content

Hide code cell content

kwargs = dict(
    panel="CSi",
    shapes=shapes,
    per_unit=True,
)

pv_orientation = {}
for slope in range(0, 100, 10):
    print(f"slope: {slope}°")
    pv_orientation[slope] = (
        cutout.pv(orientation=dict(slope=slope, azimuth=180.0), **kwargs)
        .to_pandas()
        .squeeze()
    )

df = pd.DataFrame(pv_orientation)

ax = df.mean().plot(
    ylabel="time-averaged capacity factor [per-unit]", xlabel="slope [°]"
)

ax.set_ylim(0, 0.13)

best = df.mean().idxmax()
print("The optimal slope is {:.0f}°".format(best))

slope: 0°
slope: 10°
slope: 20°
slope: 30°
slope: 40°
slope: 50°
slope: 60°
slope: 70°
slope: 80°
slope: 90°
The optimal slope is 60°

Task 2 Solution:

Show code cell content

Hide code cell content

kwargs = dict(
    panel="CSi",
    shapes=shapes,
    per_unit=True,
)

pv_azimuth = {}
for angle in [0, 90, 180, 270]:  # North, East, South, West
    print(f"azimuth: {angle}°")
    pv_azimuth[angle] = (
        cutout.pv(orientation=dict(slope=30, azimuth=angle), **kwargs)
        .to_pandas()
        .squeeze()
    )

df = pd.DataFrame(pv_azimuth)

df.loc["2013-01-01"].plot(ylabel="capacity factor [per-unit]");

azimuth: 0°
azimuth: 90°
azimuth: 180°
azimuth: 270°

Task 3 Solution:

Show code cell content

Hide code cell content

kwargs = dict(
    panel="CSi",
    orientation={"slope": 30.0, "azimuth": 180.0},
    shapes=shapes,
    per_unit=True,
)

pv = cutout.pv(**kwargs)

pv_h = cutout.pv(tracking="horizontal", **kwargs)

pv_v = cutout.pv(tracking="vertical", **kwargs)

pv_d = cutout.pv(tracking="dual", **kwargs)

DAY = "2013-01-31"
fig, ax = plt.subplots()
pv.loc[DAY].plot(ax=ax, label="fixed")
pv_v.loc[DAY].plot(ax=ax, label="vertical")
pv_h.loc[DAY].plot(ax=ax, label="horizontal")
pv_d.loc[DAY].plot(ax=ax, label="dual")
plt.legend()

<matplotlib.legend.Legend at 0x7de4effcc4a0>