Input data

Weather forcing

The first example uses data from ECMWF ERA5 obtained from CDS

A common format for multidimensional data is NetCDF.

An excellent python package for working with NetCDF in Python is Xarray.

To work with NetCDF data first we need to install the netcdf4 and xarray packages

Opening a NetCDF file with Xarray returns an xarray.Dataset

import xarray as xr

ds = xr.open_dataset("data/era5_small.nc")
ds

<xarray.Dataset> Size: 10MB
Dimensions:    (time: 12, latitude: 121, longitude: 121)
Coordinates:
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
Data variables:
    u10        (time, latitude, longitude) float64 1MB ...
    v10        (time, latitude, longitude) float64 1MB ...
    d2m        (time, latitude, longitude) float64 1MB ...
    t2m        (time, latitude, longitude) float64 1MB ...
    msl        (time, latitude, longitude) float64 1MB ...
    sp         (time, latitude, longitude) float64 1MB ...
    tcc        (time, latitude, longitude) float64 1MB ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • time: 12
  • latitude: 121
  • longitude: 121
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    1977-12-01 ... 1977-12-03T18:00:00

    long_name :

    time

    array(['1977-12-01T00:00:00.000000000', '1977-12-01T06:00:00.000000000',
           '1977-12-01T12:00:00.000000000', '1977-12-01T18:00:00.000000000',
           '1977-12-02T00:00:00.000000000', '1977-12-02T06:00:00.000000000',
           '1977-12-02T12:00:00.000000000', '1977-12-02T18:00:00.000000000',
           '1977-12-03T00:00:00.000000000', '1977-12-03T06:00:00.000000000',
           '1977-12-03T12:00:00.000000000', '1977-12-03T18:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    70.0 69.75 69.5 ... 40.5 40.25 40.0

    units :

    degrees_north

    long_name :

    latitude

    array([70.  , 69.75, 69.5 , 69.25, 69.  , 68.75, 68.5 , 68.25, 68.  , 67.75,
           67.5 , 67.25, 67.  , 66.75, 66.5 , 66.25, 66.  , 65.75, 65.5 , 65.25,
           65.  , 64.75, 64.5 , 64.25, 64.  , 63.75, 63.5 , 63.25, 63.  , 62.75,
           62.5 , 62.25, 62.  , 61.75, 61.5 , 61.25, 61.  , 60.75, 60.5 , 60.25,
           60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,
           57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,
           55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,
           52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,
           50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,
           47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,
           45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,
           42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,
           40.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    0.0 0.25 0.5 ... 29.5 29.75 30.0

    units :

    degrees_east

    long_name :

    longitude

    array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ,  2.25,
            2.5 ,  2.75,  3.  ,  3.25,  3.5 ,  3.75,  4.  ,  4.25,  4.5 ,  4.75,
            5.  ,  5.25,  5.5 ,  5.75,  6.  ,  6.25,  6.5 ,  6.75,  7.  ,  7.25,
            7.5 ,  7.75,  8.  ,  8.25,  8.5 ,  8.75,  9.  ,  9.25,  9.5 ,  9.75,
           10.  , 10.25, 10.5 , 10.75, 11.  , 11.25, 11.5 , 11.75, 12.  , 12.25,
           12.5 , 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, 24.  , 24.25, 24.5 , 24.75,
           25.  , 25.25, 25.5 , 25.75, 26.  , 26.25, 26.5 , 26.75, 27.  , 27.25,
           27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  ], dtype=float32)

• Data variables: (7)
  • u10
    (time, latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre U wind component

    [175692 values with dtype=float64]

  • v10
    (time, latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre V wind component

    [175692 values with dtype=float64]

  • d2m
    (time, latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre dewpoint temperature

    [175692 values with dtype=float64]

  • t2m
    (time, latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre temperature

    [175692 values with dtype=float64]

  • msl
    (time, latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Mean sea level pressure

    standard_name :

    air_pressure_at_mean_sea_level

    [175692 values with dtype=float64]

  • sp
    (time, latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Surface pressure

    standard_name :

    surface_air_pressure

    [175692 values with dtype=float64]

  • tcc
    (time, latitude, longitude)
    float64
    ...

    units :

    (0 - 1)

    long_name :

    Total cloud cover

    standard_name :

    cloud_area_fraction

    [175692 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmwf/mars-client/bin/grib_to_netcdf -S param -o /cache/data9/adaptor.mars.internal-1632204209.9956095-11632-13-5e63f48f-40a7-426e-b47b-db1b4f9d43e4.nc /cache/tmp/5e63f48f-40a7-426e-b47b-db1b4f9d43e4-adaptor.mars.internal-1632204158.805679-11632-16-tmp.grib

One of the datam variables is the t2m (Air temperature at 2m above the ground) which we can access in two different ways:

1. Using the string “t2m”

2. Or as a property of the Dataset object .t2m

ds["t2m"]

<xarray.DataArray 't2m' (time: 12, latitude: 121, longitude: 121)> Size: 1MB
[175692 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
Attributes:
    units:      K
    long_name:  2 metre temperature

xarray.DataArray

't2m'

• time: 12
• latitude: 121
• longitude: 121

• ...

  [175692 values with dtype=float64]

• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    1977-12-01 ... 1977-12-03T18:00:00

    long_name :

    time

    array(['1977-12-01T00:00:00.000000000', '1977-12-01T06:00:00.000000000',
           '1977-12-01T12:00:00.000000000', '1977-12-01T18:00:00.000000000',
           '1977-12-02T00:00:00.000000000', '1977-12-02T06:00:00.000000000',
           '1977-12-02T12:00:00.000000000', '1977-12-02T18:00:00.000000000',
           '1977-12-03T00:00:00.000000000', '1977-12-03T06:00:00.000000000',
           '1977-12-03T12:00:00.000000000', '1977-12-03T18:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    70.0 69.75 69.5 ... 40.5 40.25 40.0

    units :

    degrees_north

    long_name :

    latitude

    array([70.  , 69.75, 69.5 , 69.25, 69.  , 68.75, 68.5 , 68.25, 68.  , 67.75,
           67.5 , 67.25, 67.  , 66.75, 66.5 , 66.25, 66.  , 65.75, 65.5 , 65.25,
           65.  , 64.75, 64.5 , 64.25, 64.  , 63.75, 63.5 , 63.25, 63.  , 62.75,
           62.5 , 62.25, 62.  , 61.75, 61.5 , 61.25, 61.  , 60.75, 60.5 , 60.25,
           60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,
           57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,
           55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,
           52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,
           50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,
           47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,
           45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,
           42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,
           40.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    0.0 0.25 0.5 ... 29.5 29.75 30.0

    units :

    degrees_east

    long_name :

    longitude

    array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ,  2.25,
            2.5 ,  2.75,  3.  ,  3.25,  3.5 ,  3.75,  4.  ,  4.25,  4.5 ,  4.75,
            5.  ,  5.25,  5.5 ,  5.75,  6.  ,  6.25,  6.5 ,  6.75,  7.  ,  7.25,
            7.5 ,  7.75,  8.  ,  8.25,  8.5 ,  8.75,  9.  ,  9.25,  9.5 ,  9.75,
           10.  , 10.25, 10.5 , 10.75, 11.  , 11.25, 11.5 , 11.75, 12.  , 12.25,
           12.5 , 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, 24.  , 24.25, 24.5 , 24.75,
           25.  , 25.25, 25.5 , 25.75, 26.  , 26.25, 26.5 , 26.75, 27.  , 27.25,
           27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  ], dtype=float32)

• Attributes: (2)

  units :

  K

  long_name :

  2 metre temperature

ds.t2m

<xarray.DataArray 't2m' (time: 12, latitude: 121, longitude: 121)> Size: 1MB
[175692 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
Attributes:
    units:      K
    long_name:  2 metre temperature

xarray.DataArray

't2m'

• time: 12
• latitude: 121
• longitude: 121

• ...

  [175692 values with dtype=float64]

• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    1977-12-01 ... 1977-12-03T18:00:00

    long_name :

    time

    array(['1977-12-01T00:00:00.000000000', '1977-12-01T06:00:00.000000000',
           '1977-12-01T12:00:00.000000000', '1977-12-01T18:00:00.000000000',
           '1977-12-02T00:00:00.000000000', '1977-12-02T06:00:00.000000000',
           '1977-12-02T12:00:00.000000000', '1977-12-02T18:00:00.000000000',
           '1977-12-03T00:00:00.000000000', '1977-12-03T06:00:00.000000000',
           '1977-12-03T12:00:00.000000000', '1977-12-03T18:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    70.0 69.75 69.5 ... 40.5 40.25 40.0

    units :

    degrees_north

    long_name :

    latitude

    array([70.  , 69.75, 69.5 , 69.25, 69.  , 68.75, 68.5 , 68.25, 68.  , 67.75,
           67.5 , 67.25, 67.  , 66.75, 66.5 , 66.25, 66.  , 65.75, 65.5 , 65.25,
           65.  , 64.75, 64.5 , 64.25, 64.  , 63.75, 63.5 , 63.25, 63.  , 62.75,
           62.5 , 62.25, 62.  , 61.75, 61.5 , 61.25, 61.  , 60.75, 60.5 , 60.25,
           60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,
           57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,
           55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,
           52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,
           50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,
           47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,
           45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,
           42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,
           40.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    0.0 0.25 0.5 ... 29.5 29.75 30.0

    units :

    degrees_east

    long_name :

    longitude

    array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ,  2.25,
            2.5 ,  2.75,  3.  ,  3.25,  3.5 ,  3.75,  4.  ,  4.25,  4.5 ,  4.75,
            5.  ,  5.25,  5.5 ,  5.75,  6.  ,  6.25,  6.5 ,  6.75,  7.  ,  7.25,
            7.5 ,  7.75,  8.  ,  8.25,  8.5 ,  8.75,  9.  ,  9.25,  9.5 ,  9.75,
           10.  , 10.25, 10.5 , 10.75, 11.  , 11.25, 11.5 , 11.75, 12.  , 12.25,
           12.5 , 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, 24.  , 24.25, 24.5 , 24.75,
           25.  , 25.25, 25.5 , 25.75, 26.  , 26.25, 26.5 , 26.75, 27.  , 27.25,
           27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  ], dtype=float32)

• Attributes: (2)

  units :

  K

  long_name :

  2 metre temperature

To select the first timestep, we can use the xarray.Dataset.isel method.

ds0 = ds.isel(time=0)
ds0

<xarray.Dataset> Size: 821kB
Dimensions:    (latitude: 121, longitude: 121)
Coordinates:
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
    time       datetime64[ns] 8B 1977-12-01
Data variables:
    u10        (latitude, longitude) float64 117kB ...
    v10        (latitude, longitude) float64 117kB ...
    d2m        (latitude, longitude) float64 117kB ...
    t2m        (latitude, longitude) float64 117kB ...
    msl        (latitude, longitude) float64 117kB ...
    sp         (latitude, longitude) float64 117kB ...
    tcc        (latitude, longitude) float64 117kB ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • latitude: 121
  • longitude: 121
• Coordinates: (3)
  • latitude
    (latitude)
    float32
    70.0 69.75 69.5 ... 40.5 40.25 40.0

    units :

    degrees_north

    long_name :

    latitude

    array([70.  , 69.75, 69.5 , 69.25, 69.  , 68.75, 68.5 , 68.25, 68.  , 67.75,
           67.5 , 67.25, 67.  , 66.75, 66.5 , 66.25, 66.  , 65.75, 65.5 , 65.25,
           65.  , 64.75, 64.5 , 64.25, 64.  , 63.75, 63.5 , 63.25, 63.  , 62.75,
           62.5 , 62.25, 62.  , 61.75, 61.5 , 61.25, 61.  , 60.75, 60.5 , 60.25,
           60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,
           57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,
           55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,
           52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,
           50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,
           47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,
           45.  , 44.75, 44.5 , 44.25, 44.  , 43.75, 43.5 , 43.25, 43.  , 42.75,
           42.5 , 42.25, 42.  , 41.75, 41.5 , 41.25, 41.  , 40.75, 40.5 , 40.25,
           40.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    0.0 0.25 0.5 ... 29.5 29.75 30.0

    units :

    degrees_east

    long_name :

    longitude

    array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ,  2.25,
            2.5 ,  2.75,  3.  ,  3.25,  3.5 ,  3.75,  4.  ,  4.25,  4.5 ,  4.75,
            5.  ,  5.25,  5.5 ,  5.75,  6.  ,  6.25,  6.5 ,  6.75,  7.  ,  7.25,
            7.5 ,  7.75,  8.  ,  8.25,  8.5 ,  8.75,  9.  ,  9.25,  9.5 ,  9.75,
           10.  , 10.25, 10.5 , 10.75, 11.  , 11.25, 11.5 , 11.75, 12.  , 12.25,
           12.5 , 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, 24.  , 24.25, 24.5 , 24.75,
           25.  , 25.25, 25.5 , 25.75, 26.  , 26.25, 26.5 , 26.75, 27.  , 27.25,
           27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  ], dtype=float32)

  • time
    ()
    datetime64[ns]
    1977-12-01

    long_name :

    time

    array('1977-12-01T00:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (7)
  • u10
    (latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre U wind component

    [14641 values with dtype=float64]

  • v10
    (latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre V wind component

    [14641 values with dtype=float64]

  • d2m
    (latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre dewpoint temperature

    [14641 values with dtype=float64]

  • t2m
    (latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre temperature

    [14641 values with dtype=float64]

  • msl
    (latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Mean sea level pressure

    standard_name :

    air_pressure_at_mean_sea_level

    [14641 values with dtype=float64]

  • sp
    (latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Surface pressure

    standard_name :

    surface_air_pressure

    [14641 values with dtype=float64]

  • tcc
    (latitude, longitude)
    float64
    ...

    units :

    (0 - 1)

    long_name :

    Total cloud cover

    standard_name :

    cloud_area_fraction

    [14641 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmwf/mars-client/bin/grib_to_netcdf -S param -o /cache/data9/adaptor.mars.internal-1632204209.9956095-11632-13-5e63f48f-40a7-426e-b47b-db1b4f9d43e4.nc /cache/tmp/5e63f48f-40a7-426e-b47b-db1b4f9d43e4-adaptor.mars.internal-1632204158.805679-11632-16-tmp.grib

The dataset from the first timestep is also a Dataset, except it no longer has a time dimension.

To plot a map of the air temperature we can use select the t2m property and call the plot method.

ds0.t2m.plot()

<matplotlib.collections.QuadMesh at 0x7f3bf51e6e40>

Pressure

Atmospheric pressure comes in different flavors.

Surface pressure (sp) is the pressure at ground level, i.e. following the terrain.

In the following plot you can clearly see the decrease in surface pressure with altitude in areas with mountains.

This is NOT the variable used by MIKE 21 FM HD.

ds0.sp.plot()

<matplotlib.collections.QuadMesh at 0x7f3bf5091a90>

MIKE 21 FM HD expects the atmospheric pressure reduced to Mean Sea Level. (msl)

In the following plot the pressure field has a much smaller range, since the variation with altitude have been removed.

ds0.msl.plot()

<matplotlib.collections.QuadMesh at 0x7f3bf1f53b10>

Wind

The wind velocity components U,V at 10m above the ground, are found in the two data variables, u10 and v10

ds0.u10.long_name

'10 metre U wind component'

ds0.v10.long_name

'10 metre V wind component'

Convert to dfs2

Sub-region

In this case the NetCDF file covers a much larger area, than is required for our modelling work, so first step is to create a spatial subset to match our Area of Interest.

Since the latitude dimension is ordered in decreasing order in this dataset, the slice has to go from high to low.

ds_aoi = ds.sel(longitude=slice(10,30), 
                latitude=slice(60,45)) # N -> S
ds_aoi

<xarray.Dataset> Size: 3MB
Dimensions:    (time: 12, latitude: 61, longitude: 81)
Coordinates:
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
  * latitude   (latitude) float32 244B 60.0 59.75 59.5 59.25 ... 45.5 45.25 45.0
  * longitude  (longitude) float32 324B 10.0 10.25 10.5 ... 29.5 29.75 30.0
Data variables:
    u10        (time, latitude, longitude) float64 474kB ...
    v10        (time, latitude, longitude) float64 474kB ...
    d2m        (time, latitude, longitude) float64 474kB ...
    t2m        (time, latitude, longitude) float64 474kB ...
    msl        (time, latitude, longitude) float64 474kB ...
    sp         (time, latitude, longitude) float64 474kB ...
    tcc        (time, latitude, longitude) float64 474kB ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • time: 12
  • latitude: 61
  • longitude: 81
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    1977-12-01 ... 1977-12-03T18:00:00

    long_name :

    time

    array(['1977-12-01T00:00:00.000000000', '1977-12-01T06:00:00.000000000',
           '1977-12-01T12:00:00.000000000', '1977-12-01T18:00:00.000000000',
           '1977-12-02T00:00:00.000000000', '1977-12-02T06:00:00.000000000',
           '1977-12-02T12:00:00.000000000', '1977-12-02T18:00:00.000000000',
           '1977-12-03T00:00:00.000000000', '1977-12-03T06:00:00.000000000',
           '1977-12-03T12:00:00.000000000', '1977-12-03T18:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    60.0 59.75 59.5 ... 45.5 45.25 45.0

    units :

    degrees_north

    long_name :

    latitude

    array([60.  , 59.75, 59.5 , 59.25, 59.  , 58.75, 58.5 , 58.25, 58.  , 57.75,
           57.5 , 57.25, 57.  , 56.75, 56.5 , 56.25, 56.  , 55.75, 55.5 , 55.25,
           55.  , 54.75, 54.5 , 54.25, 54.  , 53.75, 53.5 , 53.25, 53.  , 52.75,
           52.5 , 52.25, 52.  , 51.75, 51.5 , 51.25, 51.  , 50.75, 50.5 , 50.25,
           50.  , 49.75, 49.5 , 49.25, 49.  , 48.75, 48.5 , 48.25, 48.  , 47.75,
           47.5 , 47.25, 47.  , 46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25,
           45.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    10.0 10.25 10.5 ... 29.5 29.75 30.0

    units :

    degrees_east

    long_name :

    longitude

    array([10.  , 10.25, 10.5 , 10.75, 11.  , 11.25, 11.5 , 11.75, 12.  , 12.25,
           12.5 , 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, 24.  , 24.25, 24.5 , 24.75,
           25.  , 25.25, 25.5 , 25.75, 26.  , 26.25, 26.5 , 26.75, 27.  , 27.25,
           27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  ], dtype=float32)

• Data variables: (7)
  • u10
    (time, latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre U wind component

    [59292 values with dtype=float64]

  • v10
    (time, latitude, longitude)
    float64
    ...

    units :

    m s**-1

    long_name :

    10 metre V wind component

    [59292 values with dtype=float64]

  • d2m
    (time, latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre dewpoint temperature

    [59292 values with dtype=float64]

  • t2m
    (time, latitude, longitude)
    float64
    ...

    units :

    K

    long_name :

    2 metre temperature

    [59292 values with dtype=float64]

  • msl
    (time, latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Mean sea level pressure

    standard_name :

    air_pressure_at_mean_sea_level

    [59292 values with dtype=float64]

  • sp
    (time, latitude, longitude)
    float64
    ...

    units :

    Pa

    long_name :

    Surface pressure

    standard_name :

    surface_air_pressure

    [59292 values with dtype=float64]

  • tcc
    (time, latitude, longitude)
    float64
    ...

    units :

    (0 - 1)

    long_name :

    Total cloud cover

    standard_name :

    cloud_area_fraction

    [59292 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2021-09-21 06:03:31 GMT by grib_to_netcdf-2.20.0: /opt/ecmwf/mars-client/bin/grib_to_netcdf -S param -o /cache/data9/adaptor.mars.internal-1632204209.9956095-11632-13-5e63f48f-40a7-426e-b47b-db1b4f9d43e4.nc /cache/tmp/5e63f48f-40a7-426e-b47b-db1b4f9d43e4-adaptor.mars.internal-1632204158.805679-11632-16-tmp.grib

Make a plot to check that the area is the one we expect.

ds_aoi.isel(time=0).u10.plot(cmap='jet')

<matplotlib.collections.QuadMesh at 0x7f3bf1e5ae90>

Time

Now we are ready to convert this data to dfs2.

First step is to create a time axis, which can be understood by MIKE IO.

import pandas as pd

time = pd.DatetimeIndex(ds.time)
time

DatetimeIndex(['1977-12-01 00:00:00', '1977-12-01 06:00:00',
               '1977-12-01 12:00:00', '1977-12-01 18:00:00',
               '1977-12-02 00:00:00', '1977-12-02 06:00:00',
               '1977-12-02 12:00:00', '1977-12-02 18:00:00',
               '1977-12-03 00:00:00', '1977-12-03 06:00:00',
               '1977-12-03 12:00:00', '1977-12-03 18:00:00'],
              dtype='datetime64[ns]', freq=None)

The second step is to collect the variables that we need, and at the same time do unit conversion where it is necessary. In this case for pressure and temperature.

data = [ds_aoi.msl.values / 100.0, # Pa -> hPa
        ds_aoi.u10.values,
        ds_aoi.v10.values,
        ds_aoi.t2m.values - 273.15] # K -> °C

The third step is to create the proper EUM item information for each variable.

from mikeio import EUMType, EUMUnit, ItemInfo

items = [
        ItemInfo("MSLP", EUMType.Pressure, EUMUnit.hectopascal),
        ItemInfo("U 10m", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),
        ItemInfo("V 10m", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),
        ItemInfo("Temperature", EUMType.Temperature, EUMUnit.degree_Celsius)
]

The fourth step is to create the mikeio.Grid2D geometry.

import mikeio

g = mikeio.Grid2D(x=ds_aoi.longitude, y=ds_aoi.latitude[::-1], projection="LONG/LAT")
g

<mikeio.Grid2D>
x: [10, 10.25, ..., 30] (nx=81, dx=0.25)
y: [45, 45.25, ..., 60] (ny=61, dy=0.25)
projection: LONG/LAT

The fifth step is to gather the data, time, item and geometry information in a mikeio.Dataset. And remember to flip the dataset upside-down since the y-axis in MIKE is south-to-north.

my_ds = mikeio.Dataset(data=data, time=time, items=items, geometry=g).flipud()
my_ds

/opt/hostedtoolcache/Python/3.13.9/x64/lib/python3.13/site-packages/mikeio/dataset/_dataset.py:153: FutureWarning: Supplying data as a list of numpy arrays is deprecated. Use Dataset.from_numpy
  warnings.warn(

<mikeio.Dataset>
dims: (time:12, y:61, x:81)
time: 1977-12-01 00:00:00 - 1977-12-03 18:00:00 (12 records)
geometry: Grid2D (ny=61, nx=81)
items:
  0:  MSLP <Pressure> (hectopascal)
  1:  U 10m <Wind Velocity> (meter per sec)
  2:  V 10m <Wind Velocity> (meter per sec)
  3:  Temperature <Temperature> (degree Celsius)

Let’s plot to confirm that

my_ds.U_10m.plot();

The dimensions for at DFS2 are expected to be (t, y, x), which matches the ones used by this dataset.

This convention is recommended by the CF convention

ds.Conventions

'CF-1.6'

ds.u10.dims

('time', 'latitude', 'longitude')

my_ds.dims

('time', 'y', 'x')

The final step is to write the dataset to a dfs2 file. The Dataset already contains all the information needed.

my_ds.to_dfs("era5_aoi.dfs2")

Screenshot of U 10m from MIKE Zero

Unstructured grid forcing

E.g. roughness map

dfs = mikeio.open("data/FakeLake.dfsu")
dfs.geometry.plot();

import numpy as np
import matplotlib.pyplot as plt

z = np.linspace(0,-30)
manning = 32 + (40 - 32)*(1 - np.exp(z/10))

plt.plot(z, manning)
plt.xlim(0,-30)
plt.ylabel("Manning")
plt.xlabel("Depth (m)")
plt.title("This is only an example...");

Extract the Z values in each element.

ze = dfs.geometry.element_coordinates[:,2]

Create a new array of manning numbers based on the depth. This is only a hypothetical example of a data transformation based on the data in the file, not a best practise for bed friction.

me = 32 + (40 - 32)*(1 - np.exp(ze/10))

my_da = mikeio.DataArray(me,
        item=mikeio.ItemInfo("Roughness", mikeio.EUMType.Mannings_M),
        geometry=dfs.geometry)
my_da

<mikeio.DataArray>
name: Roughness
dims: (element:1011)
time: 2018-01-01 00:00:00 (time-invariant)
geometry: Dfsu2D (1011 elements, 798 nodes)
values: [34.21, 33.59, ..., 32.32]

my_da.plot();

my_da.to_dfs("manning.dfsu")

dfs = mikeio.open("manning.dfsu")
dfs

<mikeio.Dfsu2DH>
number of elements: 1011
number of nodes: 798
projection: PROJCS["UTM-17",GEOGCS["Unused",DATUM["UTM Projections",SPHEROID["WGS 1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-81],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1]]
items:
  0:  Roughness <Mannings M> (meter pow  1 per 3  per sec)
time: time-invariant file (1 step) at 2018-01-01 00:00:00

ds = dfs.read()
ds["Roughness"].plot(title="Manning");