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

# built in python modules
import datetime
import os

# python add-ons
import numpy as np
import pandas as pd

# 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 = 'America/Phoenix'

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


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


# In[4]:


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


# In[5]:


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


# In[6]:


data[sorted(data.columns)]


# In[7]:


data = fm.process_data(data)


# In[8]:


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


# In[9]:


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


# In[10]:


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


# In[11]:


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


# In[12]:


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


# In[13]:


data[sorted(data.columns)]


# In[14]:


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


# In[15]:


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


# In[16]:


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


total_cloud_cover = data['total_clouds']


# In[18]:


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


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


# In[20]:


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


# In[21]:


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


data[sorted(data.columns)]


# ## NAM

# In[23]:


fm = NAM()


# In[24]:


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


# In[25]:


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


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


# In[27]:


data[sorted(data.columns)]


# ## NDFD

# In[28]:


fm = NDFD()


# In[29]:


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


# In[30]:


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


# In[31]:


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


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


# In[33]:


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


# In[34]:


data[sorted(data.columns)]


# ## RAP

# In[35]:


fm = RAP(resolution=20)


# In[36]:


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


# In[37]:


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


# In[38]:


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


data[sorted(data.columns)]


# ## HRRR

# In[40]:


fm = HRRR()


# In[41]:


data_raw = fm.get_data(latitude, longitude, start, end)


# In[42]:


# The HRRR model pulls in u, v winds for 2 layers above ground (10 m, 80 m)
# They are labeled as _0, _1 in the raw data
data_raw[sorted(data_raw.columns)]


# In[43]:


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


# In[44]:


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


# In[45]:


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


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


# In[47]:


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


# In[48]:


data[sorted(data.columns)]


# ## HRRR (ESRL)

# In[49]:


# NBVAL_SKIP
fm = HRRR_ESRL()


# In[ ]:


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


# In[ ]:


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


# In[ ]:


# NBVAL_SKIP
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[ ]:


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


# ## Quick power calculation

# In[51]:


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_']

system = PVSystem(module_parameters=module,
                  inverter_parameters=inverter,
                  surface_tilt=latitude,
                  surface_azimuth=180)

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

# extract relevant data for model chain
mc.run_model(weather=fx_data)


# In[52]:


mc.total_irrad.plot();


# In[53]:


mc.cell_temperature.plot();


# In[54]:


mc.ac.plot();


# In[ ]: