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

# # Forecast Tutorial
# 
# This tutorial will walk through forecast data from Unidata forecast model data using the forecast.py module within pvlib.
# 
# Table of contents:
# 1. [Setup](#Setup)
# 2. [Intialize and Test Each Forecast Model](#Instantiate-GFS-forecast-model)
# 
# This tutorial has been tested against the following package versions:
# * Python 3.5.2
# * IPython 5.0.0
# * pandas 0.18.0
# * matplotlib 1.5.1
# * netcdf4 1.2.1
# * siphon 0.4.0
# 
# It should work with other Python and Pandas versions. It requires pvlib >= 0.3.0 and IPython >= 3.0.
# 
# Authors:
# * Derek Groenendyk (@moonraker), University of Arizona, November 2015
# * Will Holmgren (@wholmgren), University of Arizona, November 2015, January 2016, April 2016, July 2016

# ## Setup

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
try:
    import seaborn as sns
    sns.set(rc={"figure.figsize": (12, 6)})
except ImportError:
    print('We suggest you install seaborn using conda or pip and rerun this cell')

# built in python modules
import datetime
import os

# python add-ons
import numpy as np
import pandas as pd
try:
    import netCDF4
    from netCDF4 import num2date
except ImportError:
    print('We suggest you install netCDF4 using conda rerun this cell')

# for accessing UNIDATA THREDD servers
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS

import pvlib
from pvlib.forecast import GFS, HRRR_ESRL, NAM, NDFD, HRRR, RAP


# In[2]:


# Choose a location and time.
# Tucson, AZ
latitude = 32.2
longitude = -110.9 
tz = 'US/Arizona'

start = pd.Timestamp(datetime.date.today(), tz=tz) # today's date
end = start + pd.Timedelta(days=7) # 7 days from today
print(start, end)


# ## GFS (0.5 deg)

# In[4]:


from pvlib.forecast import GFS, HRRR_ESRL, NAM, NDFD, HRRR, RAP 


# In[5]:


# GFS model, defaults to 0.5 degree resolution
fm = GFS()


# In[6]:


# retrieve data
data = fm.get_data(latitude, longitude, start, end)


# In[7]:


data


# In[8]:


data = fm.process_data(data)


# In[9]:


data[['ghi', 'dni', 'dhi']].plot()


# In[10]:


cs = fm.location.get_clearsky(data.index)


# In[11]:


fig, ax = plt.subplots()
cs['ghi'].plot(ax=ax, label='ineichen')
data['ghi'].plot(ax=ax, label='gfs+liujordan')
ax.set_ylabel('ghi')
ax.legend()


# In[12]:


fig, ax = plt.subplots()
cs['dni'].plot(ax=ax, label='ineichen')
data['dni'].plot(ax=ax, label='gfs+liujordan')
ax.set_ylabel('ghi')
ax.legend()


# In[13]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[14]:


data


# In[15]:


data['temp_air'].plot()
plt.ylabel('temperature (%s)' % fm.units['temp_air'])


# In[16]:


cloud_vars = ['total_clouds', 'low_clouds', 'mid_clouds', 'high_clouds']


# In[17]:


for varname in cloud_vars:
    data[varname].plot()
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('GFS 0.5 deg')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[18]:


total_cloud_cover = data['total_clouds']


# In[19]:


total_cloud_cover.plot(color='r', linewidth=2)
plt.ylabel('Total cloud cover' + ' (%s)' % fm.units['total_clouds'])
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('GFS 0.5 deg')


# ## GFS (0.25 deg)

# In[20]:


# GFS model at 0.25 degree resolution
fm = GFS(resolution='quarter')


# In[22]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[23]:


for varname in cloud_vars:
    data[varname].plot(ls='-', linewidth=2)
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('GFS 0.25 deg')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[24]:


data


# ## NAM

# In[25]:


fm = NAM()


# In[26]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[27]:


for varname in cloud_vars:
    data[varname].plot(ls='-', linewidth=2)
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('NAM')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[28]:


data['ghi'].plot(linewidth=2, ls='-')
plt.ylabel('GHI W/m**2')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')


# In[29]:


data


# ## NDFD

# In[30]:


fm = NDFD()


# In[31]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[33]:


total_cloud_cover = data['total_clouds']
temp = data['temp_air']
wind = data['wind_speed']


# In[34]:


total_cloud_cover.plot(color='r', linewidth=2)
plt.ylabel('Total cloud cover' + ' (%s)' % fm.units['total_clouds'])
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('NDFD')
plt.ylim(0,100)


# In[36]:


temp.plot(color='r', linewidth=2)
plt.ylabel('Temperature' + ' (%s)' % fm.units['temp_air'])
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')    


# In[37]:


wind.plot(color='r', linewidth=2)
plt.ylabel('Wind Speed' + ' (%s)' % fm.units['wind_speed'])
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')   


# In[38]:


data


# ## RAP

# In[39]:


fm = RAP(resolution=20)


# In[40]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[41]:


cloud_vars = ['total_clouds', 'high_clouds', 'mid_clouds', 'low_clouds']


# In[42]:


for varname in cloud_vars:
    data[varname].plot(ls='-', linewidth=2)
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('RAP')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[43]:


data


# ## HRRR

# In[44]:


fm = HRRR()


# In[46]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[47]:


cloud_vars = ['total_clouds', 'high_clouds', 'mid_clouds', 'low_clouds']


# In[48]:


for varname in cloud_vars:
    data[varname].plot(ls='-', linewidth=2)
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time (' + str(data.index.tz) + ')')
plt.title('HRRR')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[49]:


data


# ## HRRR (ESRL)

# In[50]:


fm = HRRR_ESRL()


# In[51]:


# retrieve data
data = fm.get_processed_data(latitude, longitude, start, end)


# In[ ]:


cloud_vars = ['total_clouds','high_clouds','mid_clouds','low_clouds']


# In[ ]:


for varname in cloud_vars:
    data[varname].plot(ls='-', linewidth=2)
plt.ylabel('Cloud cover' + ' %')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')
plt.title('HRRR_ESRL')
plt.legend(bbox_to_anchor=(1.18,1.0))


# In[ ]:


data['ghi'].plot(linewidth=2, ls='-')
plt.ylabel('GHI W/m**2')
plt.xlabel('Forecast Time ('+str(data.index.tz)+')')


# ## Quick power calculation

# In[54]:


from pvlib.pvsystem import PVSystem, retrieve_sam
from pvlib.modelchain import ModelChain

sandia_modules = retrieve_sam('SandiaMod')
sapm_inverters = retrieve_sam('cecinverter')
module = sandia_modules['Canadian_Solar_CS5P_220M___2009_']
inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_']

system = PVSystem(module_parameters=module,
                  inverter_parameters=inverter)

# fx is a common abbreviation for forecast
fx_model = GFS()
fx_data = fx_model.get_processed_data(latitude, longitude, start, end)

# use a ModelChain object to calculate modeling intermediates
mc = ModelChain(system, fx_model.location,
                orientation_strategy='south_at_latitude_tilt')

# extract relevant data for model chain
irradiance = fx_data[['ghi', 'dni', 'dhi']]
weather = fx_data[['wind_speed', 'temp_air']]
mc.run_model(fx_data.index, irradiance=irradiance, weather=weather)


# In[55]:


mc.total_irrad.plot()


# In[56]:


mc.temps.plot()


# In[57]:


mc.ac.plot()


# In[ ]: