Gridded NetCDF modelresults

2D modelresults stored in NetCDF or Grib can be loaded to ModelSkill using xarray. In this way, MIKE 21 modelresults in dfsu format can easily be compared to model results from third party providers often stored in NetCDF.

import xarray as xr
import modelskill as ms

Observations

o1 = ms.PointObservation('../data/SW/HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
o2 = ms.PointObservation("../data/SW/eur_Hm0.dfs0", item=0, x=3.2760, y=51.9990, name="EPL")
o3 = ms.TrackObservation("../data/SW/Alti_c2_Dutch.dfs0", item=3, name="c2")

MIKE ModelResult

mrMIKE = ms.model_result('../data/SW/HKZN_local_2017_DutchCoast.dfsu', name='MIKE21SW', item=0)

NetCDF ModelResult

fn = "../data/SW/ERA5_DutchCoast.nc"
xr.open_dataset(fn)

<xarray.Dataset> Size: 590kB
Dimensions:    (time: 67, latitude: 11, longitude: 20)
Coordinates:
  * time       (time) datetime64[ns] 536B 2017-10-27 ... 2017-10-29T18:00:00
  * latitude   (latitude) float32 44B 55.0 54.5 54.0 53.5 ... 51.0 50.5 50.0
  * longitude  (longitude) float32 80B -1.0 -0.5 0.0 0.5 1.0 ... 7.0 7.5 8.0 8.5
Data variables:
    mwd        (time, latitude, longitude) float64 118kB ...
    mwp        (time, latitude, longitude) float64 118kB ...
    mp2        (time, latitude, longitude) float64 118kB ...
    pp1d       (time, latitude, longitude) float64 118kB ...
    swh        (time, latitude, longitude) float64 118kB ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-06-07 12:25:02 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • time: 67
  • latitude: 11
  • longitude: 20
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2017-10-27 ... 2017-10-29T18:00:00

    long_name :

    time

    array(['2017-10-27T00:00:00.000000000', '2017-10-27T01:00:00.000000000',
           '2017-10-27T02:00:00.000000000', '2017-10-27T03:00:00.000000000',
           '2017-10-27T04:00:00.000000000', '2017-10-27T05:00:00.000000000',
           '2017-10-27T06:00:00.000000000', '2017-10-27T07:00:00.000000000',
           '2017-10-27T08:00:00.000000000', '2017-10-27T09:00:00.000000000',
           '2017-10-27T10:00:00.000000000', '2017-10-27T11:00:00.000000000',
           '2017-10-27T12:00:00.000000000', '2017-10-27T13:00:00.000000000',
           '2017-10-27T14:00:00.000000000', '2017-10-27T15:00:00.000000000',
           '2017-10-27T16:00:00.000000000', '2017-10-27T17:00:00.000000000',
           '2017-10-27T18:00:00.000000000', '2017-10-27T19:00:00.000000000',
           '2017-10-27T20:00:00.000000000', '2017-10-27T21:00:00.000000000',
           '2017-10-27T22:00:00.000000000', '2017-10-27T23:00:00.000000000',
           '2017-10-28T00:00:00.000000000', '2017-10-28T01:00:00.000000000',
           '2017-10-28T02:00:00.000000000', '2017-10-28T03:00:00.000000000',
           '2017-10-28T04:00:00.000000000', '2017-10-28T05:00:00.000000000',
           '2017-10-28T06:00:00.000000000', '2017-10-28T07:00:00.000000000',
           '2017-10-28T08:00:00.000000000', '2017-10-28T09:00:00.000000000',
           '2017-10-28T10:00:00.000000000', '2017-10-28T11:00:00.000000000',
           '2017-10-28T12:00:00.000000000', '2017-10-28T13:00:00.000000000',
           '2017-10-28T14:00:00.000000000', '2017-10-28T15:00:00.000000000',
           '2017-10-28T16:00:00.000000000', '2017-10-28T17:00:00.000000000',
           '2017-10-28T18:00:00.000000000', '2017-10-28T19:00:00.000000000',
           '2017-10-28T20:00:00.000000000', '2017-10-28T21:00:00.000000000',
           '2017-10-28T22:00:00.000000000', '2017-10-28T23:00:00.000000000',
           '2017-10-29T00:00:00.000000000', '2017-10-29T01:00:00.000000000',
           '2017-10-29T02:00:00.000000000', '2017-10-29T03:00:00.000000000',
           '2017-10-29T04:00:00.000000000', '2017-10-29T05:00:00.000000000',
           '2017-10-29T06:00:00.000000000', '2017-10-29T07:00:00.000000000',
           '2017-10-29T08:00:00.000000000', '2017-10-29T09:00:00.000000000',
           '2017-10-29T10:00:00.000000000', '2017-10-29T11:00:00.000000000',
           '2017-10-29T12:00:00.000000000', '2017-10-29T13:00:00.000000000',
           '2017-10-29T14:00:00.000000000', '2017-10-29T15:00:00.000000000',
           '2017-10-29T16:00:00.000000000', '2017-10-29T17:00:00.000000000',
           '2017-10-29T18:00:00.000000000'], dtype='datetime64[ns]')

  • latitude
    (latitude)
    float32
    55.0 54.5 54.0 ... 51.0 50.5 50.0

    units :

    degrees_north

    long_name :

    latitude

    array([55. , 54.5, 54. , 53.5, 53. , 52.5, 52. , 51.5, 51. , 50.5, 50. ],
          dtype=float32)

  • longitude
    (longitude)
    float32
    -1.0 -0.5 0.0 0.5 ... 7.5 8.0 8.5

    units :

    degrees_east

    long_name :

    longitude

    array([-1. , -0.5,  0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,
            5. ,  5.5,  6. ,  6.5,  7. ,  7.5,  8. ,  8.5], dtype=float32)

• Data variables: (5)
  • mwd
    (time, latitude, longitude)
    float64
    ...

    units :

    Degree true

    long_name :

    Mean wave direction

    [14740 values with dtype=float64]

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

    units :

    s

    long_name :

    Mean wave period

    [14740 values with dtype=float64]

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

    units :

    s

    long_name :

    Mean zero-crossing wave period

    [14740 values with dtype=float64]

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

    units :

    s

    long_name :

    Peak wave period

    [14740 values with dtype=float64]

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

    units :

    m

    long_name :

    Significant height of combined wind waves and swell

    [14740 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2021-06-07 12:25:02 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data8/adaptor.mars.internal-1623068689.815175-27387-9-f370172c-9b84-4760-b5e5-32e00e3964bd.nc /cache/tmp/f370172c-9b84-4760-b5e5-32e00e3964bd-adaptor.mars.internal-1623068689.81571-27387-2-tmp.grib

mrERA5 = ms.model_result(fn, item="swh", name='ERA5')

mrERA5

<GridModelResult>: ERA5
Time: 2017-10-27 00:00:00 - 2017-10-29 18:00:00
Quantity: Significant height of combined wind waves and swell [m]

mrERA5.data  # mr contains the xr.Dataset

<xarray.Dataset> Size: 119kB
Dimensions:  (time: 67, y: 11, x: 20)
Coordinates:
  * time     (time) datetime64[ns] 536B 2017-10-27 ... 2017-10-29T18:00:00
  * y        (y) float32 44B 55.0 54.5 54.0 53.5 53.0 ... 51.5 51.0 50.5 50.0
  * x        (x) float32 80B -1.0 -0.5 0.0 0.5 1.0 1.5 ... 6.5 7.0 7.5 8.0 8.5
Data variables:
    swh      (time, y, x) float64 118kB ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-06-07 12:25:02 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • time: 67
  • y: 11
  • x: 20
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2017-10-27 ... 2017-10-29T18:00:00

    long_name :

    time

    array(['2017-10-27T00:00:00.000000000', '2017-10-27T01:00:00.000000000',
           '2017-10-27T02:00:00.000000000', '2017-10-27T03:00:00.000000000',
           '2017-10-27T04:00:00.000000000', '2017-10-27T05:00:00.000000000',
           '2017-10-27T06:00:00.000000000', '2017-10-27T07:00:00.000000000',
           '2017-10-27T08:00:00.000000000', '2017-10-27T09:00:00.000000000',
           '2017-10-27T10:00:00.000000000', '2017-10-27T11:00:00.000000000',
           '2017-10-27T12:00:00.000000000', '2017-10-27T13:00:00.000000000',
           '2017-10-27T14:00:00.000000000', '2017-10-27T15:00:00.000000000',
           '2017-10-27T16:00:00.000000000', '2017-10-27T17:00:00.000000000',
           '2017-10-27T18:00:00.000000000', '2017-10-27T19:00:00.000000000',
           '2017-10-27T20:00:00.000000000', '2017-10-27T21:00:00.000000000',
           '2017-10-27T22:00:00.000000000', '2017-10-27T23:00:00.000000000',
           '2017-10-28T00:00:00.000000000', '2017-10-28T01:00:00.000000000',
           '2017-10-28T02:00:00.000000000', '2017-10-28T03:00:00.000000000',
           '2017-10-28T04:00:00.000000000', '2017-10-28T05:00:00.000000000',
           '2017-10-28T06:00:00.000000000', '2017-10-28T07:00:00.000000000',
           '2017-10-28T08:00:00.000000000', '2017-10-28T09:00:00.000000000',
           '2017-10-28T10:00:00.000000000', '2017-10-28T11:00:00.000000000',
           '2017-10-28T12:00:00.000000000', '2017-10-28T13:00:00.000000000',
           '2017-10-28T14:00:00.000000000', '2017-10-28T15:00:00.000000000',
           '2017-10-28T16:00:00.000000000', '2017-10-28T17:00:00.000000000',
           '2017-10-28T18:00:00.000000000', '2017-10-28T19:00:00.000000000',
           '2017-10-28T20:00:00.000000000', '2017-10-28T21:00:00.000000000',
           '2017-10-28T22:00:00.000000000', '2017-10-28T23:00:00.000000000',
           '2017-10-29T00:00:00.000000000', '2017-10-29T01:00:00.000000000',
           '2017-10-29T02:00:00.000000000', '2017-10-29T03:00:00.000000000',
           '2017-10-29T04:00:00.000000000', '2017-10-29T05:00:00.000000000',
           '2017-10-29T06:00:00.000000000', '2017-10-29T07:00:00.000000000',
           '2017-10-29T08:00:00.000000000', '2017-10-29T09:00:00.000000000',
           '2017-10-29T10:00:00.000000000', '2017-10-29T11:00:00.000000000',
           '2017-10-29T12:00:00.000000000', '2017-10-29T13:00:00.000000000',
           '2017-10-29T14:00:00.000000000', '2017-10-29T15:00:00.000000000',
           '2017-10-29T16:00:00.000000000', '2017-10-29T17:00:00.000000000',
           '2017-10-29T18:00:00.000000000'], dtype='datetime64[ns]')

  • y
    (y)
    float32
    55.0 54.5 54.0 ... 51.0 50.5 50.0

    units :

    degrees_north

    long_name :

    latitude

    array([55. , 54.5, 54. , 53.5, 53. , 52.5, 52. , 51.5, 51. , 50.5, 50. ],
          dtype=float32)

  • x
    (x)
    float32
    -1.0 -0.5 0.0 0.5 ... 7.5 8.0 8.5

    units :

    degrees_east

    long_name :

    longitude

    array([-1. , -0.5,  0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,
            5. ,  5.5,  6. ,  6.5,  7. ,  7.5,  8. ,  8.5], dtype=float32)

• Data variables: (1)
  • swh
    (time, y, x)
    float64
    ...

    units :

    m

    long_name :

    Significant height of combined wind waves and swell

    [14740 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2021-06-07 12:25:02 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data8/adaptor.mars.internal-1623068689.815175-27387-9-f370172c-9b84-4760-b5e5-32e00e3964bd.nc /cache/tmp/f370172c-9b84-4760-b5e5-32e00e3964bd-adaptor.mars.internal-1623068689.81571-27387-2-tmp.grib

Test extract from XArray

1. Extract point
2. Extract track

mrERA5.extract(o1, spatial_method="nearest").data.head()

<xarray.Dataset> Size: 104B
Dimensions:  (time: 5)
Coordinates:
  * time     (time) datetime64[ns] 40B 2017-10-27 ... 2017-10-27T04:00:00
    x        float64 8B 4.242
    y        float64 8B 52.69
    z        object 8B None
Data variables:
    ERA5     (time) float64 40B 1.22 1.347 1.466 1.612 1.793
Attributes:
    gtype:               point
    modelskill_version:  1.4.dev0

xarray.Dataset

• Dimensions:
  • time: 5
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2017-10-27 ... 2017-10-27T04:00:00

    array(['2017-10-27T00:00:00.000000000', '2017-10-27T01:00:00.000000000',
           '2017-10-27T02:00:00.000000000', '2017-10-27T03:00:00.000000000',
           '2017-10-27T04:00:00.000000000'], dtype='datetime64[ns]')

  • x
    ()
    float64
    4.242

    array(4.242)

  • y
    ()
    float64
    52.69

    array(52.6887)

  • z
    ()
    object
    None

    array(None, dtype=object)

• Data variables: (1)
  • ERA5
    (time)
    float64
    1.22 1.347 1.466 1.612 1.793

    long_name :

    Significant height of combined wind waves and swell

    units :

    m

    is_directional :

    0

    kind :

    model

    color :

    #ebdc78

    array([1.22033748, 1.34657334, 1.46574705, 1.61175692, 1.79290097])

• Attributes: (2)

  gtype :

  point

  modelskill_version :

  1.4.dev0

mrERA5.extract(o3).data.head()

<xarray.Dataset> Size: 168B
Dimensions:  (time: 5)
Coordinates:
  * time     (time) datetime64[ns] 40B 2017-10-27T12:52:52.337000 ... 2017-10...
    x        (time) float64 40B 2.423 2.414 2.405 2.396 2.387
    y        (time) float64 40B 51.25 51.31 51.37 51.42 51.48
    z        float64 8B nan
Data variables:
    ERA5     (time) float64 40B 1.439 1.464 1.489 1.514 1.538
Attributes:
    gtype:               track
    modelskill_version:  1.4.dev0

xarray.Dataset

• Dimensions:
  • time: 5
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2017-10-27T12:52:52.337000 ... 2...

    array(['2017-10-27T12:52:52.337000000', '2017-10-27T12:52:53.280000000',
           '2017-10-27T12:52:54.224000000', '2017-10-27T12:52:55.167000000',
           '2017-10-27T12:52:56.111000000'], dtype='datetime64[ns]')

  • x
    (time)
    float64
    2.423 2.414 2.405 2.396 2.387

    array([2.42285442, 2.4137888 , 2.40471053, 2.39561939, 2.38651562])

  • y
    (time)
    float64
    51.25 51.31 51.37 51.42 51.48

    array([51.25335312, 51.3102684 , 51.36718369, 51.42409897, 51.48101425])

  • z
    ()
    float64
    nan

    array(nan)

• Data variables: (1)
  • ERA5
    (time)
    float64
    1.439 1.464 1.489 1.514 1.538

    long_name :

    Significant height of combined wind waves and swell

    units :

    m

    kind :

    model

    color :

    #ebdc78

    array([1.43880867, 1.46419095, 1.48919941, 1.51383323, 1.53809179])

• Attributes: (2)

  gtype :

  track

  modelskill_version :

  1.4.dev0

Multi-file ModelResult

Use mfdataset to load multiple files as a single ModelResult.

fn = "../data/SW/CMEMS_DutchCoast_*.nc"
mrCMEMS = ms.model_result(fn, item="VHM0", name='CMEMS')
mrCMEMS

<GridModelResult>: CMEMS
Time: 2017-10-28 00:00:00 - 2017-10-29 18:00:00
Quantity: Spectral significant wave height (Hm0) [m]

Connect multiple models and observations and extract

ms.plotting.temporal_coverage(obs=[o1,o2,o3], mod=[mrERA5, mrCMEMS, mrMIKE])

# o1 is slightly outside the model domain of mrERA5, 
# we therefore use "nearest" instead of the default spatial interpolation method  
cc = ms.match(
    obs=[o1, o2, o3], 
    mod=[mrERA5, mrCMEMS, mrMIKE], 
    spatial_method='nearest',
)

Analysis and plotting

Which model is better?

sk = cc.skill()
sk.swaplevel().sort_index(level="observation").style()

		n	bias	rmse	urmse	mae	cc	si	r2
observation	model
EPL	CMEMS	43	-0.441	0.519	0.273	0.443	0.920	0.090	0.445
	ERA5	43	-0.247	0.335	0.226	0.269	0.948	0.074	0.769
	MIKE21SW	43	-0.078	0.205	0.189	0.174	0.973	0.062	0.913
HKNA	CMEMS	242	-0.742	0.882	0.476	0.742	0.903	0.128	0.222
	ERA5	242	-0.551	0.654	0.352	0.556	0.954	0.094	0.572
	MIKE21SW	242	-0.230	0.411	0.341	0.296	0.949	0.092	0.831
c2	CMEMS	36	-0.200	0.412	0.360	0.339	0.930	0.093	0.807
	ERA5	36	-0.540	0.653	0.368	0.573	0.951	0.095	0.513
	MIKE21SW	36	0.315	0.405	0.254	0.349	0.962	0.065	0.813

sk["urmse"].plot.bar(figsize=(6,3));

cc.mean_skill().style()

	n	bias	rmse	urmse	mae	cc	si	r2
model
CMEMS	321	-0.461	0.604	0.370	0.508	0.918	0.103	0.491
ERA5	321	-0.446	0.547	0.315	0.466	0.951	0.088	0.618
MIKE21SW	321	0.002	0.340	0.262	0.273	0.962	0.073	0.852

cc.plot.taylor(figsize=(6,6))