from fmskill import ModelResult, PointObservation, TrackObservation, Connector
%matplotlib inline
o1 = PointObservation('../tests/testdata/SW/HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
o2 = PointObservation("../tests/testdata/SW/eur_Hm0.dfs0", item=0, x=3.2760, y=51.9990, name="EPL")
o3 = TrackObservation("../tests/testdata/SW/Alti_c2_Dutch.dfs0", item=3, name="c2")
mrMIKE = ModelResult('../tests/testdata/SW/HKZN_local_2017_DutchCoast.dfsu', name='MIKE21SW', item=0)
fn = "../tests/testdata/SW/ERA5_DutchCoast.nc"
mr = ModelResult(fn, name='ERA5')
mr
<XArrayModelResult> 'ERA5' - Item: 0: mwd - Item: 1: mwp - Item: 2: mp2 - Item: 3: pp1d - Item: 4: swh
mr.ds # mr contains the xr.Dataset
<xarray.Dataset>
Dimensions: (time: 67, x: 20, y: 11)
Coordinates:
* x (x) float32 -1.0 -0.5 0.0 0.5 1.0 1.5 ... 6.0 6.5 7.0 7.5 8.0 8.5
* y (y) float32 55.0 54.5 54.0 53.5 53.0 52.5 52.0 51.5 51.0 50.5 50.0
* time (time) datetime64[ns] 2017-10-27 ... 2017-10-29T18:00:00
Data variables:
mwd (time, y, x) float32 332.3 330.0 327.9 325.3 ... nan nan nan nan
mwp (time, y, x) float32 4.44 4.744 5.034 5.258 ... nan nan nan nan
mp2 (time, y, x) float32 3.381 3.684 3.921 4.097 ... nan nan nan nan
pp1d (time, y, x) float32 5.221 5.343 5.766 6.099 ... nan nan nan nan
swh (time, y, x) float32 0.8702 1.222 1.393 1.519 ... nan nan nan nan
Attributes:
Conventions: CF-1.6
history: 2021-06-07 12:25:02 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...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)array([55. , 54.5, 54. , 53.5, 53. , 52.5, 52. , 51.5, 51. , 50.5, 50. ],
dtype=float32)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]')array([[[332.30542 , 329.98718 , ..., 300.44882 , 296.3397 ],
[ nan, nan, ..., 300.48727 , 293.57098 ],
...,
[254.97388 , 260.30255 , ..., nan, nan],
[262.857 , 266.2904 , ..., nan, nan]],
[[337.59564 , 334.65665 , ..., 302.95383 , 298.92163 ],
[ nan, nan, ..., 302.70663 , 295.9002 ],
...,
[259.29724 , 265.43344 , ..., nan, nan],
[266.2904 , 270.53137 , ..., nan, nan]],
...,
[[ 6.086761, 0.428497, ..., 324.37836 , 321.93927 ],
[ nan, nan, ..., 326.2901 , 320.60986 ],
...,
[328.73468 , 346.539 , ..., nan, nan],
[336.1399 , 347.01147 , ..., nan, nan]],
[[ 6.388901, 0.813034, ..., 324.5871 , 322.2249 ],
[ nan, nan, ..., 326.4439 , 320.9779 ],
...,
[345.50073 , 357.07 , ..., nan, nan],
[351.181 , 356.75687 , ..., nan, nan]]], dtype=float32)array([[[4.439693, 4.744041, ..., 6.447595, 5.899815],
[ nan, nan, ..., 6.212597, 5.65404 ],
...,
[6.205334, 5.801177, ..., nan, nan],
[8.043139, 7.669205, ..., nan, nan]],
[[4.692028, 4.95186 , ..., 6.47407 , 5.90614 ],
[ nan, nan, ..., 6.262033, 5.689652],
...,
[5.881071, 5.423494, ..., nan, nan],
[7.568693, 7.100104, ..., nan, nan]],
...,
[[8.567489, 8.596307, ..., 8.454559, 7.229434],
[ nan, nan, ..., 8.422227, 7.11299 ],
...,
[4.306613, 4.284824, ..., nan, nan],
[4.577692, 4.576755, ..., nan, nan]],
[[8.497904, 8.516882, ..., 8.453856, 7.236698],
[ nan, nan, ..., 8.429021, 7.125174],
...,
[4.224845, 4.223673, ..., nan, nan],
[4.499672, 4.534816, ..., nan, nan]]], dtype=float32)array([[[3.380553, 3.684037, ..., 5.08779 , 4.612505],
[ nan, nan, ..., 4.909118, 4.441819],
...,
[3.484005, 3.248498, ..., nan, nan],
[4.819781, 4.322208, ..., nan, nan]],
[[3.580028, 3.852309, ..., 5.113235, 4.611577],
[ nan, nan, ..., 4.946449, 4.466707],
...,
[3.281373, 3.106972, ..., nan, nan],
[4.192011, 3.778203, ..., nan, nan]],
...,
[[6.337016, 6.533519, ..., 6.508817, 5.483768],
[ nan, nan, ..., 6.475199, 5.262192],
...,
[3.198537, 3.282859, ..., nan, nan],
[3.567398, 3.636304, ..., nan, nan]],
[[6.29504 , 6.475757, ..., 6.519774, 5.503084],
[ nan, nan, ..., 6.495072, 5.28968 ],
...,
[3.162505, 3.262986, ..., nan, nan],
[3.538053, 3.628318, ..., nan, nan]]], dtype=float32)array([[[ 5.22063 , 5.343474, ..., 7.585318, nan],
[ nan, nan, ..., 7.295176, 6.304007],
...,
[ 9.963221, 10.035907, ..., nan, nan],
[10.52278 , 10.567233, ..., nan, nan]],
[[ 5.540807, 5.720419, ..., 7.616554, nan],
[ nan, nan, ..., 7.36666 , 6.315721],
...,
[ 9.935287, 10.010077, ..., nan, nan],
[10.460908, 10.50536 , ..., nan, nan]],
...,
[[10.247355, 10.30983 , ..., 10.341066, nan],
[ nan, nan, ..., 10.370501, 10.155747],
...,
[ 5.663051, 4.05856 , ..., nan, nan],
[ 5.607786, 4.396759, ..., nan, nan]],
[[10.210412, 10.25036 , ..., 10.25967 , nan],
[ nan, nan, ..., 10.299618, 10.117903],
...,
[ 3.90598 , 4.224956, ..., nan, nan],
[ 4.281122, 4.438207, ..., nan, nan]]], dtype=float32)array([[[0.870232, 1.221927, ..., 2.530895, 2.027717],
[ nan, nan, ..., 2.351693, 1.757944],
...,
[0.459568, 0.494349, ..., nan, nan],
[0.565677, 0.532132, ..., nan, nan]],
[[0.957802, 1.355401, ..., 2.620055, 2.093749],
[ nan, nan, ..., 2.434144, 1.798728],
...,
[0.475458, 0.521892, ..., nan, nan],
[0.583685, 0.557202, ..., nan, nan]],
...,
[[2.931142, 3.54908 , ..., 4.048197, 2.953388],
[ nan, nan, ..., 3.946856, 2.480754],
...,
[0.793961, 0.929201, ..., nan, nan],
[1.092866, 1.187499, ..., nan, nan]],
[[2.875528, 3.490111, ..., 4.120231, 3.031778],
[ nan, nan, ..., 4.013946, 2.543254],
...,
[0.79449 , 0.937322, ..., nan, nan],
[1.098163, 1.188381, ..., nan, nan]]], dtype=float32)We can select an item from the NetCDF file be using square brackets refering either to the item name or number.
mr['swh']
<XArrayModelResultItem> 'ERA5' - Item: swh
mr['swh']._extract_point(o1).head()
| ERA5 | |
|---|---|
| time | |
| 2017-10-27 00:00:00 | 1.220337 |
| 2017-10-27 01:00:00 | 1.346573 |
| 2017-10-27 02:00:00 | 1.465747 |
| 2017-10-27 03:00:00 | 1.611757 |
| 2017-10-27 04:00:00 | 1.792901 |
mr['swh'].extract_observation(o1) # will call _extract_point
<PointComparer> Observation: HKNA, n_points=386 Model: ERA5, rmse=0.545
mr['swh']._extract_track(o3).head()
| x | y | ERA5 | |
|---|---|---|---|
| time | |||
| 2017-10-27 12:52:52.337 | 2.422854 | 51.253353 | 1.438809 |
| 2017-10-27 12:52:53.280 | 2.413789 | 51.310268 | 1.464191 |
| 2017-10-27 12:52:54.224 | 2.404711 | 51.367184 | 1.489200 |
| 2017-10-27 12:52:55.167 | 2.395619 | 51.424099 | 1.513833 |
| 2017-10-27 12:52:56.111 | 2.386516 | 51.481014 | 1.538092 |
Use mfdataset to load multiple files as a single ModelResult.
fn = "../tests/testdata/SW/CMEMS_DutchCoast_*.nc"
mr2 = ModelResult(fn, name='CMEMS')
mr2
<XArrayModelResult> 'CMEMS' - Item: 0: VHM0 - Item: 1: VMDR - Item: 2: VTPK - Item: 3: VTM02 - Item: 4: VHM0_WW - Item: 5: VHM0_SW1 - Item: 6: VHM0_SW2
fn = r"../tests/testdata/SW/NWW3_hs_201710.grib"
#mr3 = ModelResult(fn, name='WW3', engine='cfgrib') # not yet supported
mrERA5 = mr['swh']
mrCMEMS = mr2['VHM0']
con = Connector([o1,o2,o3], [mrERA5, mrCMEMS, mrMIKE])
con.plot_temporal_coverage();
cc = con.extract()
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.087 | 0.445 |
| ERA5 | 43 | -0.247 | 0.335 | 0.226 | 0.269 | 0.948 | 0.072 | 0.769 | |
| MIKE21SW | 43 | -0.078 | 0.205 | 0.189 | 0.174 | 0.973 | 0.061 | 0.913 | |
| HKNA | CMEMS | 242 | -0.742 | 0.882 | 0.476 | 0.742 | 0.903 | 0.124 | 0.222 |
| ERA5 | 242 | -0.551 | 0.654 | 0.352 | 0.556 | 0.954 | 0.091 | 0.572 | |
| MIKE21SW | 242 | -0.230 | 0.411 | 0.341 | 0.296 | 0.949 | 0.088 | 0.831 | |
| c2 | CMEMS | 41 | -0.311 | 0.475 | 0.358 | 0.402 | 0.937 | 0.087 | 0.741 |
| ERA5 | 41 | -0.412 | 0.711 | 0.580 | 0.620 | 0.887 | 0.141 | 0.417 | |
| MIKE21SW | 41 | 0.327 | 0.410 | 0.247 | 0.358 | 0.964 | 0.060 | 0.806 |
sk.plot_bar('urmse', figsize=(6,3));
cc.mean_skill().style()
| n | bias | rmse | urmse | mae | cc | si | r2 | |
|---|---|---|---|---|---|---|---|---|
| model | ||||||||
| CMEMS | 326 | -0.498 | 0.625 | 0.369 | 0.529 | 0.920 | 0.099 | 0.469 |
| ERA5 | 326 | -0.403 | 0.567 | 0.386 | 0.482 | 0.930 | 0.101 | 0.586 |
| MIKE21SW | 326 | 0.006 | 0.342 | 0.259 | 0.276 | 0.962 | 0.070 | 0.850 |
cc.taylor(figsize=6)