Alpha coefficients - Wind scaling power law

In the following Notebook, we will have a look at how to compute wind speed at a desired height (e.g. 84 meters and 100 meters) for climate projection models, starting from wind speed at a reference height (e.g. 10 meters) making use of wind speed scaling coeffiecient (alpha) through a power law. The alpha coefficients file contains data for each latitude-longitude grid point and for each month and hour of the day. It has been produced starting from ERA5 data of wind speed at 10 meters and at 100 meters, both bias adjusted with the delta method using data coming from the Global Wind Atlas version 2 (GWA).

For illustration purposes, this notebook uses a short temporal subsample of ws10 from the C3S Energy Lot1 dataset. This approach is intended to demonstrate the workflow and methodology.

It is important to emphasize that for rigorous climate impact studies, the use of such limited data is not sufficient. Reliable analyses should be based on longer time series (typically at least 30 years).
We will consider just one year (2020) of data from a CMIP6 climate model (ECE3) and a specific climate change scenario (SSP126). The analysis could be extended to consider more years and also others climate models and/or scenarios.

After reading the needed data from the CDS, we will apply the wind speed scaling power law to compute wind speed data at 84 meters and 100 meters, then we will visualize it by generating figures that will help us understanding the data generated, namely a map of mean wind speed for may at 4:00 PM and a map of mean wind speed for december at 9:00 PM. All the code will be in Python.

Learning objectives 🎯

Throughout this notebook, you will:

• understand the Wind Shear Power Law

• learn how to read C3S Energy Lot1 NetCDF files of wind speed and alpha scaling coefficients

• how to define a function in order to apply the power_law to compute wind speed at any desired height

• how to plot a map of average wind speed at a specified month and hour of the day.

Target Audience 🎯

Anyone interested in learning what scaling cofficients are and how to use them to compute wind speed at a desired height given a dataset of wind speed at a reference height, making use of wind speed data from Lot1 climate projection models.

The physical meaning of the Wind Speed Power Law

The wind speed power law is an empirical relationship used to estimate how wind speed changes with height above ground level within the atmospheric boundary layer (ABL). It reflects how friction and turbulence near the surface slow down the wind and how this effect decreases with height. The physical idea is that near the ground, roughness elements (such as trees, buildings, terrain) create friction, which dissipates momentum. As you climb, the effects of friction weaken, so the wind accelerates. The rate of this increase (wind shear) depends on the roughness of the surface and the stability of the atmosphere (which varies depending on the time of day and season).

Mathematical formula

The power law expresses the wind speed \(w(z)\) at height \(h\) as: $\(w(z)=w_{ref}\left(\frac{h(z)}{h_{ref}}\right)^{\alpha}\)$ where:

• \(w_{ref}\) : known wind speed at reference height \(h_{ref}\), the typical refernce height is 10 meters

• \(\alpha\) : wind shear exponent (dimensionless)

Prepare your environment

Import required packages

We start by importing some useful libraries. These include: os, which provides a way to interact with the operating system and it is used here to manipulate paths; pandas, one of the most common and easy to use tools for data analysis and manipulation; numpy, fundamental for numerical computing; matplotlib, a standard library for creating visualizations in Python; xarray, a library to work with multidimensional arrays in NetCDF format.

# importing libraries
import os
import glob
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
import dask

Load the data

Now we load alpha coefficients file into a Xarray datarray. For each grid point, it contains data for each month and hour of the day

# define alpha file path
alpha_file = "./data/alph/ANCI_ALPH-coeff_C3S2LOT1_025d_v1.00.nc"
# load into alpha datarray
alpha = xr.open_dataset(alpha_file).alpha
# show result datarray
alpha

<xarray.DataArray 'alpha' (time: 288, latitude: 721, longitude: 1440)> Size: 1GB
[299013120 values with dtype=float32]
Coordinates:
  * longitude  (longitude) float32 6kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * time       (time) datetime64[ns] 2kB 2011-01-01 ... 2011-12-01T23:00:00
Attributes:
    units:          -
    long_name:      alpha power law coefficient
    short_name:     alpha
    standard_name:  alpha_power_law_coefficient

xarray.DataArray

'alpha'

• time: 288
• latitude: 721
• longitude: 1440

• ...

  [299013120 values with dtype=float32]

• Coordinates: (3)
  • longitude
    (longitude)
    float32
    0.0 0.25 0.5 ... 359.2 359.5 359.8

    units :

    degrees_east

    long_name :

    longitude

    short_name :

    longitude

    axis :

    X

    array([0.0000e+00, 2.5000e-01, 5.0000e-01, ..., 3.5925e+02, 3.5950e+02,
           3.5975e+02], dtype=float32)

  • latitude
    (latitude)
    float32
    90.0 89.75 89.5 ... -89.75 -90.0

    units :

    degrees_north

    long_name :

    latitude

    short_name :

    latitude

    standard_name :

    latitude

    axis :

    Y

    array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)

  • time
    (time)
    datetime64[ns]
    2011-01-01 ... 2011-12-01T23:00:00

    long_name :

    time

    short_name :

    time

    standard_name :

    time

    axis :

    T

    array(['2011-01-01T00:00:00.000000000', '2011-01-01T01:00:00.000000000',
           '2011-01-01T02:00:00.000000000', ..., '2011-12-01T21:00:00.000000000',
           '2011-12-01T22:00:00.000000000', '2011-12-01T23:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (3)
  • longitude
    PandasIndex

    PandasIndex(Index([   0.0,   0.25,    0.5,   0.75,    1.0,   1.25,    1.5,   1.75,    2.0,
             2.25,
           ...
            357.5, 357.75,  358.0, 358.25,  358.5, 358.75,  359.0, 359.25,  359.5,
           359.75],
          dtype='float32', name='longitude', length=1440))

  • latitude
    PandasIndex

    PandasIndex(Index([  90.0,  89.75,   89.5,  89.25,   89.0,  88.75,   88.5,  88.25,   88.0,
            87.75,
           ...
           -87.75,  -88.0, -88.25,  -88.5, -88.75,  -89.0, -89.25,  -89.5, -89.75,
            -90.0],
          dtype='float32', name='latitude', length=721))

  • time
    PandasIndex

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

• Attributes: (4)

  units :

  -

  long_name :

  alpha power law coefficient

  short_name :

  alpha

  standard_name :

  alpha_power_law_coefficient

We can notice that longitude values of alpha spans [0,359.75]. We want to convert it to [-180, 180]; to do this, we will define a function to handle the calculation:

# function to convert longitude values from 0, 360 to -180,180
def shift_longitude(ds):
    ds = ds.assign_coords(longitude=(((ds.longitude + 180) % 360) - 180))
    return ds.sortby('longitude')
# convert longitude
alpha = shift_longitude(alpha)
alpha

<xarray.DataArray 'alpha' (time: 288, latitude: 721, longitude: 1440)> Size: 1GB
[299013120 values with dtype=float32]
Coordinates:
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * time       (time) datetime64[ns] 2kB 2011-01-01 ... 2011-12-01T23:00:00
  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8
Attributes:
    units:          -
    long_name:      alpha power law coefficient
    short_name:     alpha
    standard_name:  alpha_power_law_coefficient

xarray.DataArray

'alpha'

• time: 288
• latitude: 721
• longitude: 1440

• ...

  [299013120 values with dtype=float32]

• Coordinates: (3)
  • latitude
    (latitude)
    float32
    90.0 89.75 89.5 ... -89.75 -90.0

    units :

    degrees_north

    long_name :

    latitude

    short_name :

    latitude

    standard_name :

    latitude

    axis :

    Y

    array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)

  • time
    (time)
    datetime64[ns]
    2011-01-01 ... 2011-12-01T23:00:00

    long_name :

    time

    short_name :

    time

    standard_name :

    time

    axis :

    T

    array(['2011-01-01T00:00:00.000000000', '2011-01-01T01:00:00.000000000',
           '2011-01-01T02:00:00.000000000', ..., '2011-12-01T21:00:00.000000000',
           '2011-12-01T22:00:00.000000000', '2011-12-01T23:00:00.000000000'],
          dtype='datetime64[ns]')

  • longitude
    (longitude)
    float32
    -180.0 -179.8 ... 179.5 179.8

    array([-180.  , -179.75, -179.5 , ...,  179.25,  179.5 ,  179.75],
          dtype=float32)

• Indexes: (3)
  • latitude
    PandasIndex

    PandasIndex(Index([  90.0,  89.75,   89.5,  89.25,   89.0,  88.75,   88.5,  88.25,   88.0,
            87.75,
           ...
           -87.75,  -88.0, -88.25,  -88.5, -88.75,  -89.0, -89.25,  -89.5, -89.75,
            -90.0],
          dtype='float32', name='latitude', length=721))

  • time
    PandasIndex

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

  • longitude
    PandasIndex

    PandasIndex(Index([ -180.0, -179.75,  -179.5, -179.25,  -179.0, -178.75,  -178.5, -178.25,
            -178.0, -177.75,
           ...
             177.5,  177.75,   178.0,  178.25,   178.5,  178.75,   179.0,  179.25,
             179.5,  179.75],
          dtype='float32', name='longitude', length=1440))

• Attributes: (4)

  units :

  -

  long_name :

  alpha power law coefficient

  short_name :

  alpha

  standard_name :

  alpha_power_law_coefficient

Relation to surface roughness

Surface roughness determines the amount of kinetic energy that is subtracted from the flow: smooth surfaces (e.g., open water, ice, short grass) have low friction, which translates into low alpha values, while rough surfaces (e.g., forests, cities) have high friction, which means high alpha values. This is because rougher surfaces create more turbulence, extending the layer of slowed air higher into the atmosphere.

Let’s plot the scaling coeffients for a given month and hour

plt.figure(figsize=(12,6))
alpha.mean(dim="time").plot()
title=plt.title(f"alpha monthly average")

Relation to Atmospheric Stability, Daily & Seasonal Cycles

The exponent \(\alpha\) also varies with atmospheric stability, which changes over day and night and across seasons.

Let’s select one location over land and one over sea. We will then we plot the seasonal cycle by taking the average over time after applying the groupby method over month. The exponent should change across seasons:

• In winter, long nights and stable stratification → higher \(\alpha\)

• In summer, stronger solar heating → more convective mixing → lower \(\alpha\)

It varies strongly also between land and sea because of the surface roughness:

• Over the sea: water surfaces are very smooth → low surface roughness → wind changes slowly with height → lower \(\alpha\)

• over land: Llnd has higher and more variable roughness (vegetation, buildings, terrain) → higher \(\alpha\).

land_lat, land_lon = 53, 11 # near Hamburg
sea_lat, sea_lon = 55, 6 # North Sea

alpha.sel(latitude=land_lat, longitude=land_lon, method='nearest').groupby('time.month').mean().plot(label='land')
alpha.sel(latitude=sea_lat, longitude=sea_lon, method='nearest').groupby('time.month').mean().plot(label='sea')
plt.title("alpha average seasonal cycle")
plt.legend()

<matplotlib.legend.Legend at 0x14a0efaa0>

One can see how the average (over hours of the day) scaling coefficient varies as a function of the month for the given location. Indeed values are lower in the summer and higher in the winter season for the land location; for the sea location, there are some months of lag in the seasonal cycle.

Next we will plot the average over all months for each hour of the day. The exponent should vary in the following way:

• Daytime (unstable) Surface heating causes convective turbulence, enhancing mixing. The vertical temperature gradient is small (or negative), so momentum is well mixed. Wind profile is more uniform → small \(\alpha\) (≈ 0.1)

• Nighttime (stable) The ground cools, suppressing turbulence. Wind speed near the surface drops sharply, but remains strong aloft. Wind shear is large → large \(\alpha\) (≈ 0.4)

alpha.sel(latitude=land_lat, longitude=land_lon, method='nearest').groupby('time.hour').mean().plot(label='land')
alpha.sel(latitude=sea_lat, longitude=sea_lon, method='nearest').groupby('time.hour').mean().plot(label='sea')
plt.title("alpha average diurnal cycle")
plt.legend()

<matplotlib.legend.Legend at 0x14a184110>

Indeed we see that the coefficient has the maximum values in the night and lowest during the sunny ones. This is less evident for the sea location.

Now we are going to use this alpha data in order to compute wind speed at 100 meters.

First of all let’s load the wind speed at 10 meters files for the selected year into a dataset using open_mfdataset method. We will use the year 2030 (considering only the months of May and december for illustrative purpouses), model “ECE3” and scenario SP126. The data is already available in the fiolder “data/nc/ws10”

We can see that the resulting dataset covers one year of data at hourly resolution.

# year, model, scenario
year = 2030
model = "ECE3"
scenario = "SP126"

# WS10 folder
ws10_dir = f"./data/ws10/"
print(f"{ws10_dir}/*S{year}*.nc")
# ws10 files
ws10_files = sorted(glob.glob(f"{ws10_dir}/*S{year}*.nc"))
ws10 = xr.open_mfdataset(ws10_files)
ws10 = shift_longitude(ws10)
ws10

./data/ws10//*S2030*.nc

<xarray.Dataset> Size: 12GB
Dimensions:    (time: 1488, latitude: 721, longitude: 1440)
Coordinates:
  * time       (time) datetime64[ns] 12kB 2030-05-01 ... 2030-12-31T23:00:00
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8
Data variables:
    ws10       (time, latitude, longitude) float64 12GB dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
Attributes:
    title:        bias adjusted 10 metre wind speed
    institution:  Copernicus Climate Change Service, Sectoral Information Sys...
    source:       Institution=AEMET, Spain; BSC, Spain; CNR-ISAC, Italy; DMI,...
    history:      2025-01-08 18:25:20 GMT by write_hist_xr_to_nc-0.1.dev51+g1...
    references:   
    Conventions:  CF-1.12-draft

xarray.Dataset

• Dimensions:
  • time: 1488
  • latitude: 721
  • longitude: 1440
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2030-05-01 ... 2030-12-31T23:00:00

    long_name :

    time

    short_name :

    time

    standard_name :

    time

    axis :

    T

    array(['2030-05-01T00:00:00.000000000', '2030-05-01T01:00:00.000000000',
           '2030-05-01T02:00:00.000000000', ..., '2030-12-31T21:00:00.000000000',
           '2030-12-31T22:00:00.000000000', '2030-12-31T23:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    90.0 89.75 89.5 ... -89.75 -90.0

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    short_name :

    latitude

    array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    -180.0 -179.8 ... 179.5 179.8

    array([-180.  , -179.75, -179.5 , ...,  179.25,  179.5 ,  179.75],
          dtype=float32)

• Data variables: (1)
  • ws10
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>

    units :

    m s**-1

    long_name :

    10 metre wind speed

    short_name :

    ws10

    standard_name :

    10_metre_wind_speed

    Array Chunk Bytes 11.51 GiB 7.92 MiB Shape (1488, 721, 1440) (1, 721, 1440) Dask graph 1488 chunks in 6 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • time
    PandasIndex

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

  • latitude
    PandasIndex

    PandasIndex(Index([  90.0,  89.75,   89.5,  89.25,   89.0,  88.75,   88.5,  88.25,   88.0,
            87.75,
           ...
           -87.75,  -88.0, -88.25,  -88.5, -88.75,  -89.0, -89.25,  -89.5, -89.75,
            -90.0],
          dtype='float32', name='latitude', length=721))

  • longitude
    PandasIndex

    PandasIndex(Index([ -180.0, -179.75,  -179.5, -179.25,  -179.0, -178.75,  -178.5, -178.25,
            -178.0, -177.75,
           ...
             177.5,  177.75,   178.0,  178.25,   178.5,  178.75,   179.0,  179.25,
             179.5,  179.75],
          dtype='float32', name='longitude', length=1440))

• Attributes: (6)

  title :

  bias adjusted 10 metre wind speed

  institution :

  Copernicus Climate Change Service, Sectoral Information System for the Energy Sector

  source :

  Institution=AEMET, Spain; BSC, Spain; CNR-ISAC, Italy; DMI, Denmark; ENEA, Italy; FMI, Finland; Geomar, Germany; ICHEC, Ireland; ICTP, Italy; IDL, Portugal; IMAU, The Netherlands; IPMA, Portugal; KIT, Karlsruhe, Germany; KNMI, The Netherlands; Lund University, Sweden; Met Eireann, Ireland; NLeSC, The Netherlands; NTNU, Norway; Oxford University, UK; surfSARA, The Netherlands; SMHI, Sweden; Stockholm University, Sweden; Unite ASTR, Belgium; University College Dublin, Ireland; University of Bergen, Norway; University of Copenhagen, Denmark; University of Helsinki, Finland; University of Santiago de Compostela, Spain; Uppsala University, Sweden; Utrecht University, The Netherlands; Vrije Universiteit Amsterdam, the Netherlands; Wageningen University, The Netherlands. Mailing address: EC-Earth consortium, Rossby Center, Swedish Meteorological and Hydrological Institute/SMHI, SE-601 76 Norrkoping, Sweden; project=CMIP6; source=EC-Earth3; variant=r1i1p1f1; experiment=ssp126; spatial_resolution=100 km; frequency=3hrPt; variable=uas;

  history :

  2025-01-08 18:25:20 GMT by write_hist_xr_to_nc-0.1.dev51+g18f5059: Saved file with scale factor and add_offset;2025-01-08 18:25:20 GMT by cdft-0.1.dev5+gee082ea: Bias adjustment with cdft method;2024-08-30 06:03:41 GMT by write_hist_xr_to_nc-0.1.dev51+g18f5059: Saved file with scale factor and add_offset;2024-08-30 06:02:22 GMT by calc_ws-0.1.dev51+g18f5059: WIND:calc_ws: calculation of wind speed from components;Original Conventions=CF-1.12-draft; 2024-08-30 06:02:15 GMT by spline-0.1.dev51+g18f5059: Temporal interpolation with cubic spline and moving windows;2024-08-29 22:20:50 GMT by write_hist_xr_to_nc-0.1.dev51+g18f5059: Saved file with scale factor and add_offset;2024-08-29 22:18:39 GMT by standardize_cmip6_variables-0.1.dev5+gee082ea: Standardize internal variables according to c3s standard;Thu Aug 29 20:39:23 2024: cdo remapbil,lot1_grid.txt uas_3hr_EC-Earth3_ssp126_r1i1p1f1_gr_203001010000-203012312100.nc P_CMI6_ECEC_ECE3_U--_0010m_Glob_025d_S203001010000_E203012312100_INS_MAP_03h_NA-_noc_org_NA_SP126_NA---_NA---_v1.00.nc 2019-09-26T14:59:10Z ; CMOR rewrote data to be consistent with CMIP6, CF-1.7 CMIP-6.2 and CF standards.; processed by ece2cmor v1.2.1, git rev. 37fbd03c29f5012299f660ded8d41abfdfda8dda The cmor-fixer version v2.0 script has been applied.The cmor-fixer version v2.1 script has been applied.;2019-09-26T15:12:31Z altered by CMOR: Treated scalar dimension: 'height'. 2019-09-26T15:12:31Z altered by CMOR: Reordered dimensions, original order: lat lon time.

  references :

  Conventions :

  CF-1.12-draft

Define power law function

Let’s define a function called power_law, which will allow us to calculate the wind speed at a desired height, given a dataset of wind speeds at a reference height and an array of scaling coefficients. The function accepts the two data sets as input, along with optional parameters representing the input and output heights and the names of the input and output variables. The function then applies the power law by looping over each month-hour pair in the input wind speed data set. The function will implement the following formula:

def power_law(
    ds_in: xr.Dataset,
    alpha: xr.DataArray,
    h_in: float = 10.0,
    h_out: float = 100.0,
    var_in: str = "ws10",
    var_out: str | None = None,
) -> xr.Dataset:
    """
    Scale wind speed from a reference height h_in to h_out using a power-law profile.

    Example:
        ws100 = power_law(ws10, alpha)              # 10 -> 100 m
        ws84  = power_law(ws10, alpha, 10, 84)      # 10 -> 84 m
    """

    # If no output name is given, build one from the output height
    if var_out is None:
        # Simple convention: ws<h_out>
        var_out = f"ws{int(h_out)}"

    # list of months in ds_in.time datarray
    months = ds_in["time"].dt.month
    # list of hours in ds_in.time datarray
    hours = ds_in["time"].dt.hour
    # list of unique month-hour pairs
    month_hour_pairs = list(zip(months.values, hours.values))
    unique_month_hour_pairs = list(set(month_hour_pairs))
    unique_month_hour_pairs.sort()

    list_ds = []
    for month, hour in unique_month_hour_pairs:
        alphai = alpha.sel(
            time=(
                (alpha["time"].dt.month == month)
                & (alpha["time"].dt.hour == hour)
            )
        ).values
        dsi = ds_in.sel(
            time=(
                (ds_in["time"].dt.month == month)
                & (ds_in["time"].dt.hour == hour)
            )
        )

        # rename input variable to output name
        dsi = dsi.rename({var_in: var_out})

        # apply power-law scaling
        dsi = dsi * ((h_out / h_in) ** alphai)
        list_ds.append(dsi)

    ds_out = xr.concat(list_ds, dim="time")
    ds_out = ds_out.sortby("time")

    return ds_out

Apply power law to compute wind speed at 84 and 100 meters

Now that we have loaded wind speed at 10 meters data in a xarray dataset, we can call the function to compute

• wind speed at 100 meters (first case) where the function uses the defaults h_in=10, h_out=100 → var_out=”ws100”

• wind speed at 100 meters (second case) where h_in=10, h_out=84 → the output variable is called var_out=”ws84”

ws100 = power_law(ws10, alpha)            
ws84  = power_law(ws10, alpha, 10, 84) 

print(ws100)
print(ws84)

<xarray.Dataset> Size: 12GB
Dimensions:    (time: 1488, latitude: 721, longitude: 1440)
Coordinates:
  * time       (time) datetime64[ns] 12kB 2030-05-01 ... 2030-12-31T23:00:00
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8
Data variables:
    ws100      (time, latitude, longitude) float64 12GB dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
Attributes:
    title:        bias adjusted 10 metre wind speed
    institution:  Copernicus Climate Change Service, Sectoral Information Sys...
    source:       Institution=AEMET, Spain; BSC, Spain; CNR-ISAC, Italy; DMI,...
    history:      2025-01-08 18:25:20 GMT by write_hist_xr_to_nc-0.1.dev51+g1...
    references:   
    Conventions:  CF-1.12-draft
<xarray.Dataset> Size: 12GB
Dimensions:    (time: 1488, latitude: 721, longitude: 1440)
Coordinates:
  * time       (time) datetime64[ns] 12kB 2030-05-01 ... 2030-12-31T23:00:00
  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8
Data variables:
    ws84       (time, latitude, longitude) float64 12GB dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
Attributes:
    title:        bias adjusted 10 metre wind speed
    institution:  Copernicus Climate Change Service, Sectoral Information Sys...
    source:       Institution=AEMET, Spain; BSC, Spain; CNR-ISAC, Italy; DMI,...
    history:      2025-01-08 18:25:20 GMT by write_hist_xr_to_nc-0.1.dev51+g1...
    references:   
    Conventions:  CF-1.12-draft

Plot the data

Now that we have the wind speed at 84 and 100 meters, we can inspect the values using matplotlib plotting library.

Plot a map of average wind speed for selected month/hour

In the following we are going to specify a month and hour, and we will calculate the average wind speed at 100 for those values. In order to do so we will make use of .sel method of xarray together with .dt accessor, then we will use .mean method over time dimension in order to have a map of longitude latitude average values. We are going to plot data for May and hour 4:00 PM

# Compute average of ws100 for month 5 and hour 16
month, hour = 5, 16
plt.figure(figsize=(12,6))
ds1 = ws100.sel(time=ws100.time.dt.month==5)
ds1 = ds1.sel(time=ds1.time.dt.hour==16)
ds1.mean(dim="time").ws100.plot()
title=plt.title(f"ws100 average year {year} month {month} hour {hour}") 

In the same way, we will calculate the average wind speed at 84 meters for those values.

# Compute average of ws100 for month 5 and hour 16
month, hour = 5, 16
plt.figure(figsize=(12,6))
ds1 = ws84.sel(time=ws84.time.dt.month==5)
ds1 = ds1.sel(time=ds1.time.dt.hour==16)
ds1.mean(dim="time").ws84.plot()
title=plt.title(f"ws84 average year {year} month {month} hour {hour}") 

Now we are going to plot wind speeds at 10 meters, at 84 meters and at 100 meters to compare them, both for a land location and sea location. We will select a specific month and group values based on hour of the day, then we will take the average of those values:

ws10.ws10.sel(time=ws10.time.dt.month==5).sel(latitude=land_lat, longitude=land_lon, method='nearest').groupby('time.hour').mean().plot(label='10 m')
ws84.ws84.sel(time=ws84.time.dt.month==5).sel(latitude=land_lat, longitude=land_lon, method='nearest').groupby('time.hour').mean().plot(label='84 m')
ws100.ws100.sel(time=ws100.time.dt.month==5).sel(latitude=land_lat, longitude=land_lon, method='nearest').groupby('time.hour').mean().plot(label='100 m')
plt.title(f"average wind speed month 5 lat {land_lat} lon {land_lon} land")
plt.legend()

<matplotlib.legend.Legend at 0x14b7dd640>

We can see that wind speed at 100 meters is higher than at 10 meters especially at night, when alpha coefficients are higher due to boundary layer stability

ws10.ws10.sel(time=ws10.time.dt.month==5).sel(latitude=sea_lat, longitude=sea_lon, method='nearest').groupby('time.hour').mean().plot(label='10 m')
ws84.ws84.sel(time=ws84.time.dt.month==5).sel(latitude=sea_lat, longitude=sea_lon, method='nearest').groupby('time.hour').mean().plot(label='84 m')
ws100.ws100.sel(time=ws100.time.dt.month==5).sel(latitude=sea_lat, longitude=sea_lon, method='nearest').groupby('time.hour').mean().plot(label='100 m')
plt.title(f"average wind speed month 5 lat {sea_lat} lon {sea_lon}")
plt.legend()

<matplotlib.legend.Legend at 0x14c0702c0>

Now we do the same plot, selecting December and hour 9:00 PM; here there is the plot for ws100

# Compute average of ws100 for month 12 and hour 21
month, hour = 12, 21
plt.figure(figsize=(12,6))
ds1 = ws100.sel(time=ws100.time.dt.month==12)
ds1 = ds1.sel(time=ds1.time.dt.hour==21)
ds1.mean(dim="time").ws100.plot()
title=plt.title(f"ws100 average year {year} month {month} hour {hour}") 

# Compute average of ws100 for month 12 and hour 21
month, hour = 12, 21
ds1 = ws84.sel(time=ws84.time.dt.month==12)
ds1 = ds1.sel(time=ds1.time.dt.hour==21)
plt.figure(figsize=(12,6))
ds1.mean(dim="time").ws84.plot()
title=plt.title(f"ws84 average year {year} month {month} hour {hour}") 

Take home messages 📌

• Wind shear coefficient can be used to estimate the vertical wind profile and derive wind speed information at different heights above the ground;

• Alpha coefficients are strongly dependent on the hour of the day, terrain type and season;

• Using alpha coefficients stratified by hour of the day and month is useful for more realistic wind power modeling;

• Visualization of average wind speed map for given month and hour of the day using matplotlib