import xarray as xr
import modelskill as ms
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.
Observations
= ms.PointObservation('../data/SW/HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
o1 = ms.PointObservation("../data/SW/eur_Hm0.dfs0", item=0, x=3.2760, y=51.9990, name="EPL")
o2 = ms.TrackObservation("../data/SW/Alti_c2_Dutch.dfs0", item=3, name="c2") o3
MIKE ModelResult
= ms.model_result('../data/SW/HKZN_local_2017_DutchCoast.dfsu', name='MIKE21SW', item=0) mrMIKE
NetCDF ModelResult
= "../data/SW/ERA5_DutchCoast.nc"
fn xr.open_dataset(fn)
<xarray.Dataset> Size: 590kB Dimensions: (longitude: 20, latitude: 11, time: 67) Coordinates: * longitude (longitude) float32 80B -1.0 -0.5 0.0 0.5 1.0 ... 7.0 7.5 8.0 8.5 * latitude (latitude) float32 44B 55.0 54.5 54.0 53.5 ... 51.0 50.5 50.0 * time (time) datetime64[ns] 536B 2017-10-27 ... 2017-10-29T18:00:00 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...
= ms.model_result(fn, item="swh", name='ERA5') mrERA5
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]
# mr contains the xr.Dataset mrERA5.data
<xarray.Dataset> Size: 119kB Dimensions: (time: 67, y: 11, x: 20) Coordinates: * 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 * y (y) float32 44B 55.0 54.5 54.0 53.5 53.0 ... 51.5 51.0 50.5 50.0 * time (time) datetime64[ns] 536B 2017-10-27 ... 2017-10-29T18:00:00 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...
Test extract from XArray
- Extract point
- Extract track
="nearest").data.head() mrERA5.extract(o1, spatial_method
<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.2.dev0
mrERA5.extract(o3).data.head()
<xarray.Dataset> Size: 168B Dimensions: (time: 5) Coordinates: 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 * time (time) datetime64[ns] 40B 2017-10-27T12:52:52.337000 ... 2017-10... 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.2.dev0
Multi-file ModelResult
Use mfdataset to load multiple files as a single ModelResult.
= "../data/SW/CMEMS_DutchCoast_*.nc"
fn = ms.model_result(fn, item="VHM0", name='CMEMS')
mrCMEMS 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
=[o1,o2,o3], mod=[mrERA5, mrCMEMS, mrMIKE]) ms.plotting.temporal_coverage(obs
# o1 is slightly outside the model domain of mrERA5,
# we therefore use "nearest" instead of the default spatial interpolation method
= ms.match(
cc =[o1, o2, o3],
obs=[mrERA5, mrCMEMS, mrMIKE],
mod='nearest',
spatial_method )
Analysis and plotting
Which model is better?
= cc.skill()
sk ="observation").style() sk.swaplevel().sort_index(level
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 |
"urmse"].plot.bar(figsize=(6,3)); sk[
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 |
=6) cc.plot.taylor(figsize