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.

In [1]:

import xarray as xr
import modelskill as ms

Observations¶

In [2]:

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

MIKE ModelResult¶

In [3]:

mrMIKE = ms.ModelResult('../tests/testdata/SW/HKZN_local_2017_DutchCoast.dfsu', name='MIKE21SW', item=0)

NetCDF ModelResult¶

In [4]:

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

Out[4]:

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

In [5]:

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

In [6]:

mrERA5

Out[6]:

<modelskill.model.grid.GridModelResult at 0x1c1986ee500>

In [7]:

mrERA5.data  # mr contains the xr.Dataset

Test extract from XArray¶

Extract point
Extract track

In [8]:

mrERA5.extract(o1).data.head()

In [9]:

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

Multi-file ModelResult¶

Use mfdataset to load multiple files as a single ModelResult.

In [10]:

fn = "../tests/testdata/SW/CMEMS_DutchCoast_*.nc"
mrCMEMS = ms.ModelResult(fn, item="VHM0", name='CMEMS')
mrCMEMS

Out[10]:

<modelskill.model.grid.GridModelResult at 0x1c19875c610>

Grib file ModelResult¶

In [11]:

fn = r"../tests/testdata/SW/NWW3_hs_201710.grib"

In [12]:

#mr3 = ModelResult(fn, name='WW3', engine='cfgrib') # not yet supported

Connect multiple models and observations and extract¶

In [13]:

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

Out[13]:

<Axes: >

In [14]:

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

Analysis and plotting¶

Which model is better?

In [15]:

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

Out[15]:

		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	33	-0.295	0.480	0.379	0.406	0.903	0.092	0.438
	ERA5	33	-0.529	0.628	0.338	0.547	0.906	0.082	0.038
	MIKE21SW	33	0.325	0.416	0.259	0.363	0.920	0.063	0.578

In [16]:

sk.plot.bar('urmse', figsize=(6,3));

In [17]:

cc.mean_skill().style()

Out[17]:

	n	bias	rmse	urmse	mae	cc	si	r2
model
CMEMS	318	-0.493	0.627	0.376	0.530	0.909	0.103	0.368
ERA5	318	-0.442	0.539	0.305	0.457	0.936	0.084	0.460
MIKE21SW	318	0.006	0.344	0.263	0.277	0.947	0.072	0.774

In [18]:

cc.plot.taylor(figsize=6)

Out[18]: