#!/usr/bin/env python
# coding: utf-8

# # Multi model comparison

# In[1]:


import numpy as np
from fmskill.model import ModelResult
from fmskill import PointObservation, TrackObservation, Connector
get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')
get_ipython().run_line_magic('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


# In[5]:


con.modelresults


# 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


# ## 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()


# In[11]:


cc.skill(observation="c2")


# In[12]:


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


# 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')


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


# 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()


# In[24]:


s.style(columns='rmse')


# In[25]:


s.plot_bar('rmse');


# In[26]:


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


# In[27]:


s.style()


# 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()


# In[31]:


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


# In[32]:


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


# In[33]:


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


# In[34]:


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


# In[ ]: