Multi model comparison¶

In [1]:

import numpy as np
from fmskill.model import ModelResult
from fmskill import PointObservation, TrackObservation, Connector
%load_ext autoreload
%autoreload 2
%matplotlib inline

Define observations¶

In [2]:

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")

Define models¶

In [3]:

mr1 = ModelResult('../tests/testdata/SW/HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
mr2 = ModelResult('../tests/testdata/SW/HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)

Connect observations and model results¶

In [4]:

con = Connector([o1, o2, o3], [mr1, mr2])
con

Out[4]:

<Connector> with 
 -<PointConnector> obs=HKNA(n=564) :: 2 models=[SW_1, SW_2]
 -<PointConnector> obs=EPL(n=95) :: 2 models=[SW_1, SW_2]
 -<TrackConnector> obs=c2(n=298) :: 2 models=[SW_1, SW_2]

In [5]:

con.modelresults

Out[5]:

{'SW_1': <DfsModelResultItem> 'SW_1'
 File: ../tests/testdata/SW/HKZN_local_2017_DutchCoast.dfsu
 - Item: 0: Sign. Wave Height <Significant wave height> (meter),
 'SW_2': <DfsModelResultItem> 'SW_2'
 File: ../tests/testdata/SW/HKZN_local_2017_DutchCoast_v2.dfsu
 - Item: 0: Sign. Wave Height <Significant wave height> (meter)}

In [6]:

con.plot_observation_positions();

In [7]:

con.plot_temporal_coverage();

In [8]:

cc = con.extract()    # returns a collection of comparisons

In [9]:

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

Out[9]:

<PointComparer>
Observation: EPL, n_points=67
 Model: SW_1, rmse=0.224
 Model: SW_2, rmse=0.232

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 [10]:

cc.skill()

Out[10]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	EPL	67	-0.066597	0.223597	0.213449	0.188513	0.969846	0.078296	0.932596
	HKNA	386	-0.194260	0.351964	0.293499	0.251839	0.971194	0.088669	0.905300
	c2	113	-0.001210	0.351796	0.351794	0.294585	0.974335	0.118511	0.899507
SW_2	EPL	67	-0.000199	0.232479	0.232479	0.198294	0.969846	0.085276	0.927134
	HKNA	386	-0.100426	0.293033	0.275287	0.214422	0.971194	0.083167	0.934358
	c2	113	0.081431	0.430268	0.422492	0.357138	0.974335	0.142327	0.849675

In [11]:

cc.skill(observation="c2")

Out[11]:

	observation	n	bias	rmse	urmse	mae	cc	si	r2
model
SW_1	c2	113	-0.001210	0.351796	0.351794	0.294585	0.974335	0.118511	0.899507
SW_2	c2	113	0.081431	0.430268	0.422492	0.357138	0.974335	0.142327	0.849675

In [12]:

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

Out[12]:

	n	bias	rmse	urmse	mae	cc	si	r2
model
SW_1	499	-0.097735	0.35188	0.322647	0.273212	0.972764	0.10359	0.902404

In [13]:

cc.scatter(model='SW_1', cmap='OrRd')

In [14]:

cc.taylor()

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 [15]:

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 [16]:

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

Out[16]:

	model	n	bias	rmse	urmse	mae	cc	si	r2
observation
EPL	SW_1	48	-0.100116	0.238728	0.216720	0.203709	0.919574	0.097755	0.812073
HKNA	SW_1	281	-0.097075	0.203241	0.178559	0.164563	0.968106	0.067642	0.916126
c2	SW_1	72	-0.188193	0.313787	0.251089	0.258661	0.478554	0.122423	-1.420894

In [17]:

cc.scatter(model='SW_2', start='2017-10-28', 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 [18]:

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

In [19]:

ax = con.plot_observation_positions();
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 [20]:

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

Out[20]:

	model	n	bias	rmse	urmse	mae	cc	si	r2
observation
HKNA	SW_1	386	-0.19426	0.351964	0.293499	0.251839	0.971194	0.088669	0.905300
c2	SW_1	42	-0.05587	0.388404	0.384365	0.336023	0.952688	0.139796	0.749645

In [21]:

cc.scatter(model="SW_2", area=polygon, 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 [22]:

s = cc.skill()

In [23]:

s.style()

Out[23]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	EPL	67	-0.067	0.224	0.213	0.189	0.970	0.078	0.933
	HKNA	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
	c2	113	-0.001	0.352	0.352	0.295	0.974	0.119	0.900
SW_2	EPL	67	-0.000	0.232	0.232	0.198	0.970	0.085	0.927
	HKNA	386	-0.100	0.293	0.275	0.214	0.971	0.083	0.934
	c2	113	0.081	0.430	0.422	0.357	0.974	0.142	0.850

In [24]:

s.style(columns='rmse')

Out[24]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	EPL	67	-0.067	0.224	0.213	0.189	0.970	0.078	0.933
	HKNA	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
	c2	113	-0.001	0.352	0.352	0.295	0.974	0.119	0.900
SW_2	EPL	67	-0.000	0.232	0.232	0.198	0.970	0.085	0.927
	HKNA	386	-0.100	0.293	0.275	0.214	0.971	0.083	0.934
	c2	113	0.081	0.430	0.422	0.357	0.974	0.142	0.850

In [25]:

s.plot_bar('rmse');

In [26]:

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

In [27]:

s.style()

Out[27]:

		n	bias	rmse	si
model	datetime
SW_1	2017-10-27 00:00:00	84	-0.096	0.257	0.114
	2017-10-27 12:00:00	156	-0.166	0.253	0.094
	2017-10-28 00:00:00	77	-0.111	0.178	0.056
	2017-10-28 12:00:00	84	-0.037	0.204	0.058
	2017-10-29 00:00:00	82	-0.421	0.596	0.088
	2017-10-29 12:00:00	83	0.007	0.419	0.105
SW_2	2017-10-27 00:00:00	84	-0.070	0.236	0.108
	2017-10-27 12:00:00	156	-0.146	0.247	0.098
	2017-10-28 00:00:00	77	-0.056	0.150	0.056
	2017-10-28 12:00:00	84	0.091	0.234	0.062
	2017-10-29 00:00:00	82	-0.226	0.477	0.088
	2017-10-29 12:00:00	83	0.173	0.472	0.110

In [28]:

s.plot_line('rmse', title='Hm0 rmse [m]');

In [29]:

s.plot_grid('si', fmt='0.1%', title='Hm0 Scatter index');

The sel() method can subset the skill object¶

A new skill object will be returned

In [30]:

s = cc.skill()
s.style()

Out[30]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	EPL	67	-0.067	0.224	0.213	0.189	0.970	0.078	0.933
	HKNA	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
	c2	113	-0.001	0.352	0.352	0.295	0.974	0.119	0.900
SW_2	EPL	67	-0.000	0.232	0.232	0.198	0.970	0.085	0.927
	HKNA	386	-0.100	0.293	0.275	0.214	0.971	0.083	0.934
	c2	113	0.081	0.430	0.422	0.357	0.974	0.142	0.850

In [31]:

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

Out[31]:

	model	n	bias	rmse	urmse	mae	cc	si	r2
observation
EPL	SW_1	67	-0.067	0.224	0.213	0.189	0.970	0.078	0.933
HKNA	SW_1	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
c2	SW_1	113	-0.001	0.352	0.352	0.295	0.974	0.119	0.900

In [32]:

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

Out[32]:

	observation	n	bias	rmse	urmse	mae	cc	si	r2
model
SW_1	HKNA	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
SW_2	HKNA	386	-0.100	0.293	0.275	0.214	0.971	0.083	0.934

In [33]:

s.sel('rmse>0.25').style()

Out[33]:

		n	bias	rmse	urmse	mae	cc	si	r2
model	observation
SW_1	HKNA	386	-0.194	0.352	0.293	0.252	0.971	0.089	0.905
SW_1	c2	113	-0.001	0.352	0.352	0.295	0.974	0.119	0.900
SW_2	HKNA	386	-0.100	0.293	0.275	0.214	0.971	0.083	0.934
SW_2	c2	113	0.081	0.430	0.422	0.357	0.974	0.142	0.850

In [34]:

s.sel('rmse>0.3', columns=['rmse','mae']).style()

Out[34]:

		rmse	mae
model	observation
SW_1	HKNA	0.352	0.252
SW_1	c2	0.352	0.295
SW_2	c2	0.430	0.357

In [ ]: