%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

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

"""
Explanation: 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
2. Intialize and Test Each 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
End of explanation
"""


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

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

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

data

data = fm.process_data(data)

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

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

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

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

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

data

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

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

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

total_cloud_cover = data['total_clouds']

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

"""
Explanation: GFS (0.5 deg)
End of explanation
"""


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

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

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

data

"""
Explanation: GFS (0.25 deg)
End of explanation
"""


fm = NAM()

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

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

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

data

"""
Explanation: NAM
End of explanation
"""


fm = NDFD()

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

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

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)

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

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

data

"""
Explanation: NDFD
End of explanation
"""


fm = RAP(resolution=20)

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

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

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

data

"""
Explanation: RAP
End of explanation
"""


fm = HRRR()

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

# 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

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

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

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

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

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

data

"""
Explanation: HRRR
End of explanation
"""


fm = HRRR_ESRL()

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

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

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

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

"""
Explanation: HRRR (ESRL)
End of explanation
"""


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
mc.run_model(fx_data.index, weather=fx_data)

mc.total_irrad.plot()

mc.temps.plot()

mc.ac.plot()

"""
Explanation: Quick power calculation
End of explanation
"""


