Metocean multi-model comparison¶

In this example, two metocean wave models "SW_1" and "SW_2" are compared to three observations: two points and one track.

In [1]:

import numpy as np
import modelskill as ms

Define observations¶

In [2]:

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

Define models¶

In [3]:

mr1 = ms.model_result(fldr + 'HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
mr2 = ms.model_result(fldr + 'HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)

Get spatial and temporal overview¶

In [4]:

ms.plotting.spatial_overview(obs=[o1, o2, o3], mod=[mr1, mr2]);

In [5]:

ms.plotting.temporal_coverage(obs=[o1, o2, o3], mod=[mr1, mr2]);

Compare¶

In [6]:

cc = ms.match(obs=[o1, o2, o3], mod=[mr1, mr2])    # returns a collection of comparisons

In [7]:

cc["EPL"]   # select a single comparer from the collection like this

Out[7]:

<Comparer>
Quantity: Significant wave height [m]
Observation: EPL, n_points=67
Model(s):
0: SW_1
1: SW_2

Perform analysis¶

You can perform simple filtering on specific observation or specific model. You can refer to observations and models using their name or index.

The main analysis methods are:

skill()
mean_skill()
scatter()
taylor()

In [8]:

cc.skill()

Out[8]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	HKNA	386	-0.202413	0.355195	0.291877	0.255866	0.971708	0.093967	0.903554
	EPL	67	-0.071238	0.224923	0.213344	0.189455	0.969760	0.082482	0.931793
	c2	113	-0.004701	0.352470	0.352439	0.294758	0.975050	0.128010	0.899121
SW_2	HKNA	386	-0.109149	0.294151	0.273151	0.215417	0.971708	0.087938	0.933856
	EPL	67	-0.005165	0.231782	0.231724	0.197612	0.969760	0.089588	0.927570
	c2	113	0.077695	0.431332	0.424276	0.359415	0.975050	0.154102	0.848931

In [9]:

cc.sel(observation="c2").skill()

Out[9]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	c2	113	-0.004701	0.352470	0.352439	0.294758	0.97505	0.128010	0.899121
SW_2	c2	113	0.077695	0.431332	0.424276	0.359415	0.97505	0.154102	0.848931

In [10]:

cc.sel(model=0, observation=[0,"c2"]).mean_skill()

Out[10]:

	n	bias	rmse	urmse	mae	cc	si	r2
model
SW_1	499	-0.103557	0.353833	0.322158	0.275312	0.973379	0.110988	0.901338
SW_2	499	-0.015727	0.362741	0.348713	0.287416	0.973379	0.121020	0.891393

In [11]:

cc.sel(model='SW_1').plot.scatter();

In [12]:

cc.plot.taylor(normalize_std=True, aggregate_observations=False)

Out[12]:

Time series plot (specifically for point comparisons)¶

If you select an comparison from the collection which is a PointComparer, you can do a time series plot

In [13]:

cc['EPL'].plot.timeseries(figsize=(12,4));

Filtering on time¶

Use the start and end arguments to do your analysis on part of the time series

In [14]:

cc.sel(model="SW_1", end='2017-10-28').skill()

Out[14]:

	n	bias	rmse	urmse	mae	cc	si	r2
observation
HKNA	281	-0.105743	0.206592	0.177479	0.168153	0.968267	0.069743	0.913338
EPL	48	-0.104505	0.239904	0.215946	0.204544	0.920404	0.100557	0.810216
c2	72	-0.192928	0.318751	0.253735	0.260704	0.475245	0.124315	-1.498094

In [15]:

cc.sel(model='SW_2', start='2017-10-28').plot.scatter(cmap='OrRd', figsize=(6,7));

Filtering on area¶

You can do you analysis in a specific area by providing a bounding box or a closed polygon

In [16]:

bbox = np.array([0.5,52.5,5,54])
polygon = np.array([[6,51],[0,55],[0,51],[6,51]])

In [17]:

ax = ms.plotting.spatial_overview(obs=[o1, o2, o3], mod=[mr1, mr2]);
ax.plot([bbox[0],bbox[2],bbox[2],bbox[0],bbox[0]],[bbox[1],bbox[1],bbox[3],bbox[3],bbox[1]]);
ax.plot(polygon[:,0],polygon[:,1]);

In [18]:

cc.sel(model="SW_1", area=bbox).skill()

Out[18]:

	n	bias	rmse	urmse	mae	cc	si	r2
observation
HKNA	386	-0.202413	0.355195	0.291877	0.255866	0.971708	0.093967	0.903554
c2	42	-0.050795	0.390762	0.387447	0.336623	0.954156	0.146893	0.746596

In [19]:

cc.sel(model="SW_2", area=polygon).plot.scatter(backend='plotly')

Skill object¶

The skill() and mean_skill() methods return a skill object that can visualize results in various ways. The primary methods of the skill object are:

style()
plot.bar()
plot.line()
plot.grid()
sel()

In [20]:

sk = cc.skill()

In [21]:

sk.style()

Out[21]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	HKNA	386	-0.202	0.355	0.292	0.256	0.972	0.094	0.904
	EPL	67	-0.071	0.225	0.213	0.189	0.970	0.082	0.932
	c2	113	-0.005	0.352	0.352	0.295	0.975	0.128	0.899
SW_2	HKNA	386	-0.109	0.294	0.273	0.215	0.972	0.088	0.934
	EPL	67	-0.005	0.232	0.232	0.198	0.970	0.090	0.928
	c2	113	0.078	0.431	0.424	0.359	0.975	0.154	0.849

In [22]:

sk.style(columns='rmse')

Out[22]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	HKNA	386	-0.202	0.355	0.292	0.256	0.972	0.094	0.904
	EPL	67	-0.071	0.225	0.213	0.189	0.970	0.082	0.932
	c2	113	-0.005	0.352	0.352	0.295	0.975	0.128	0.899
SW_2	HKNA	386	-0.109	0.294	0.273	0.215	0.972	0.088	0.934
	EPL	67	-0.005	0.232	0.232	0.198	0.970	0.090	0.928
	c2	113	0.078	0.431	0.424	0.359	0.975	0.154	0.849

In [23]:

sk["rmse"].plot.bar();

In [24]:

sk = cc.skill(by=['model','freq:12h'], metrics=['bias','rmse','si'])

In [25]:

sk.style()

Out[25]:

		n	bias	rmse	si
model	time
SW_1	2017-10-27 00:00:00	84	-0.100	0.258	0.116
	2017-10-27 12:00:00	156	-0.171	0.257	0.095
	2017-10-28 00:00:00	77	-0.120	0.184	0.056
	2017-10-28 12:00:00	84	-0.051	0.206	0.058
	2017-10-29 00:00:00	82	-0.426	0.599	0.088
	2017-10-29 12:00:00	83	0.002	0.421	0.107
SW_2	2017-10-27 00:00:00	84	-0.074	0.236	0.110
	2017-10-27 12:00:00	156	-0.152	0.251	0.099
	2017-10-28 00:00:00	77	-0.065	0.153	0.056
	2017-10-28 12:00:00	84	0.076	0.227	0.062
	2017-10-29 00:00:00	82	-0.232	0.478	0.088
	2017-10-29 12:00:00	83	0.167	0.473	0.112

In [26]:

sk["rmse"].plot.line(title='Hm0 rmse [m]');

In [27]:

sk["si"].plot.grid(fmt='0.1%', title='Hm0 Scatter index');

The sel() method can subset the skill object¶

A new skill object will be returned

In [28]:

sk = cc.skill()
sk.style()

Out[28]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	HKNA	386	-0.202	0.355	0.292	0.256	0.972	0.094	0.904
	EPL	67	-0.071	0.225	0.213	0.189	0.970	0.082	0.932
	c2	113	-0.005	0.352	0.352	0.295	0.975	0.128	0.899
SW_2	HKNA	386	-0.109	0.294	0.273	0.215	0.972	0.088	0.934
	EPL	67	-0.005	0.232	0.232	0.198	0.970	0.090	0.928
	c2	113	0.078	0.431	0.424	0.359	0.975	0.154	0.849

In [29]:

sk.sel(model='SW_1').style()

Out[29]:

	model	n	bias	rmse	urmse	mae	cc	si	r2
observation
HKNA	SW_1	386	-0.202	0.355	0.292	0.256	0.972	0.094	0.904
EPL	SW_1	67	-0.071	0.225	0.213	0.189	0.970	0.082	0.932
c2	SW_1	113	-0.005	0.352	0.352	0.295	0.975	0.128	0.899

In [30]:

sk.sel(observation='HKNA').style()

Out[30]:

	observation	n	bias	rmse	urmse	mae	cc	si	r2
model
SW_1	HKNA	386	-0.202	0.355	0.292	0.256	0.972	0.094	0.904
SW_2	HKNA	386	-0.109	0.294	0.273	0.215	0.972	0.088	0.934

In [32]:

sk[['rmse','mae']].query('rmse>0.3').style()

Out[32]:

		rmse	mae
model	observation
SW_1	HKNA	0.355	0.256
SW_1	c2	0.352	0.295
SW_2	c2	0.431	0.359