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

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


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


# In[9]:


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


# In[10]:


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


# In[11]:


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


# In[12]:


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


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


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


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


# In[22]:


sk.style(columns='rmse')


# In[23]:


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


# In[24]:


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


# In[25]:


sk.style()


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


# In[29]:


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


# In[30]:


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


# In[31]:


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


# In[32]:


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