NetCDF to dfs2

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 numpy as np
import pandas as pd
import xarray as xr
import mikeio

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

<xarray.Dataset> Size: 10MB
Dimensions:    (longitude: 121, latitude: 121, time: 12)
Coordinates:
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
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:
  • longitude: 121
  • latitude: 121
  • time: 12
• Coordinates: (3)
  • 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)

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

  • 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]')

• 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]

• 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,
           ...
           27.75,  28.0, 28.25,  28.5, 28.75,  29.0, 29.25,  29.5, 29.75,  30.0],
          dtype='float32', name='longitude', length=121))

  • latitude
    PandasIndex

    PandasIndex(Index([ 70.0, 69.75,  69.5, 69.25,  69.0, 68.75,  68.5, 68.25,  68.0, 67.75,
           ...
           42.25,  42.0, 41.75,  41.5, 41.25,  41.0, 40.75,  40.5, 40.25,  40.0],
          dtype='float32', name='latitude', length=121))

  • time
    PandasIndex

    PandasIndex(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]', name='time', freq=None))

• 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:
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
Attributes:
    units:      K
    long_name:  2 metre temperature

xarray.DataArray

't2m'

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

• ...

  [175692 values with dtype=float64]

• Coordinates: (3)
  • 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)

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

  • 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]')

• 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,
           ...
           27.75,  28.0, 28.25,  28.5, 28.75,  29.0, 29.25,  29.5, 29.75,  30.0],
          dtype='float32', name='longitude', length=121))

  • latitude
    PandasIndex

    PandasIndex(Index([ 70.0, 69.75,  69.5, 69.25,  69.0, 68.75,  68.5, 68.25,  68.0, 67.75,
           ...
           42.25,  42.0, 41.75,  41.5, 41.25,  41.0, 40.75,  40.5, 40.25,  40.0],
          dtype='float32', name='latitude', length=121))

  • time
    PandasIndex

    PandasIndex(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]', name='time', freq=None))

• 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:
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.0
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
Attributes:
    units:      K
    long_name:  2 metre temperature

xarray.DataArray

't2m'

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

• ...

  [175692 values with dtype=float64]

• Coordinates: (3)
  • 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)

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

  • 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]')

• 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,
           ...
           27.75,  28.0, 28.25,  28.5, 28.75,  29.0, 29.25,  29.5, 29.75,  30.0],
          dtype='float32', name='longitude', length=121))

  • latitude
    PandasIndex

    PandasIndex(Index([ 70.0, 69.75,  69.5, 69.25,  69.0, 68.75,  68.5, 68.25,  68.0, 67.75,
           ...
           42.25,  42.0, 41.75,  41.5, 41.25,  41.0, 40.75,  40.5, 40.25,  40.0],
          dtype='float32', name='latitude', length=121))

  • time
    PandasIndex

    PandasIndex(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]', name='time', freq=None))

• 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:    (longitude: 121, latitude: 121)
Coordinates:
  * longitude  (longitude) float32 484B 0.0 0.25 0.5 0.75 ... 29.5 29.75 30.0
  * latitude   (latitude) float32 484B 70.0 69.75 69.5 69.25 ... 40.5 40.25 40.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:
  • longitude: 121
  • latitude: 121
• Coordinates: (3)
  • 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)

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

  • 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]

• Indexes: (2)
  • longitude
    PandasIndex

    PandasIndex(Index([  0.0,  0.25,   0.5,  0.75,   1.0,  1.25,   1.5,  1.75,   2.0,  2.25,
           ...
           27.75,  28.0, 28.25,  28.5, 28.75,  29.0, 29.25,  29.5, 29.75,  30.0],
          dtype='float32', name='longitude', length=121))

  • latitude
    PandasIndex

    PandasIndex(Index([ 70.0, 69.75,  69.5, 69.25,  69.0, 68.75,  68.5, 68.25,  68.0, 67.75,
           ...
           42.25,  42.0, 41.75,  41.5, 41.25,  41.0, 40.75,  40.5, 40.25,  40.0],
          dtype='float32', name='latitude', length=121))

• 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 0x7fd221546550>

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 0x7fd2215be690>

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 0x7fd21c310b90>

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:    (longitude: 81, latitude: 61, time: 12)
Coordinates:
  * longitude  (longitude) float32 324B 10.0 10.25 10.5 ... 29.5 29.75 30.0
  * latitude   (latitude) float32 244B 60.0 59.75 59.5 59.25 ... 45.5 45.25 45.0
  * time       (time) datetime64[ns] 96B 1977-12-01 ... 1977-12-03T18:00:00
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:
  • longitude: 81
  • latitude: 61
  • time: 12
• Coordinates: (3)
  • 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)

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

  • 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]')

• 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]

• Indexes: (3)
  • longitude
    PandasIndex

    PandasIndex(Index([ 10.0, 10.25,  10.5, 10.75,  11.0, 11.25,  11.5, 11.75,  12.0, 12.25,
            12.5, 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,  24.0, 24.25,  24.5, 24.75,
            25.0, 25.25,  25.5, 25.75,  26.0, 26.25,  26.5, 26.75,  27.0, 27.25,
            27.5, 27.75,  28.0, 28.25,  28.5, 28.75,  29.0, 29.25,  29.5, 29.75,
            30.0],
          dtype='float32', name='longitude'))

  • latitude
    PandasIndex

    PandasIndex(Index([ 60.0, 59.75,  59.5, 59.25,  59.0, 58.75,  58.5, 58.25,  58.0, 57.75,
            57.5, 57.25,  57.0, 56.75,  56.5, 56.25,  56.0, 55.75,  55.5, 55.25,
            55.0, 54.75,  54.5, 54.25,  54.0, 53.75,  53.5, 53.25,  53.0, 52.75,
            52.5, 52.25,  52.0, 51.75,  51.5, 51.25,  51.0, 50.75,  50.5, 50.25,
            50.0, 49.75,  49.5, 49.25,  49.0, 48.75,  48.5, 48.25,  48.0, 47.75,
            47.5, 47.25,  47.0, 46.75,  46.5, 46.25,  46.0, 45.75,  45.5, 45.25,
            45.0],
          dtype='float32', name='latitude'))

  • time
    PandasIndex

    PandasIndex(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]', name='time', freq=None))

• 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()

<matplotlib.collections.QuadMesh at 0x7fd21c20f190>

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.

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)

Then create the spatial axis

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 next 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.

das = [
    mikeio.DataArray(np.flip(ds_aoi.msl.values, axis=1) / 100.0,
                     time=time, geometry=g,
                     item=mikeio.ItemInfo("MSLP", mikeio.EUMType.Pressure, mikeio.EUMUnit.hectopascal)),
    
    mikeio.DataArray(np.flip(ds_aoi.u10.values, axis=1),
                     time=time, geometry=g,
                     item=mikeio.ItemInfo("U 10m", mikeio.EUMType.Wind_Velocity, mikeio.EUMUnit.meter_per_sec)),
    
    mikeio.DataArray(np.flip(ds_aoi.v10.values, axis=1), 
                     time=time, geometry=g, 
                     item=mikeio.ItemInfo("V 10m", mikeio.EUMType.Wind_Velocity, mikeio.EUMUnit.meter_per_sec)),
    
    mikeio.DataArray(np.flip(ds_aoi.t2m.values, axis=1) - 273.15,
                     time=time, geometry=g,
                     item=mikeio.ItemInfo("Temperature", mikeio.EUMType.Temperature, mikeio.EUMUnit.degree_Celsius))
]

my_ds = mikeio.Dataset(das)
my_ds

<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 it looks alright

my_ds.U_10m.plot(cmap="coolwarm",vmin=-8,vmax=8);

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