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

# # Forecast to Power Tutorial
# 
# This tutorial will walk through the process of going from Unidata forecast model data to AC power using the SAPM.
# 
# Table of contents:
# 1. [Setup](#Setup)
# 2. [Load Forecast data](#Load-Forecast-data)
# 2. [Calculate modeling intermediates](#Calculate-modeling-intermediates)
# 2. [DC power using SAPM](#DC-power-using-SAPM)
# 2. [AC power using SAPM](#AC-power-using-SAPM)
# 
# 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

# These are just your standard interactive scientific python imports that you'll get very used to using.

# In[1]:


# built-in python modules
import datetime
import inspect
import os

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

# plotting stuff
# first line makes the plots appear in the notebook
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import matplotlib as mpl
# seaborn makes your plots look better
try:
    import seaborn as sns
    sns.set(rc={"figure.figsize": (12, 6)})
    sns.set_color_codes()
except ImportError:
    print('We suggest you install seaborn using conda or pip and rerun this cell')

# finally, we import the pvlib library
from pvlib import solarposition,irradiance,atmosphere,pvsystem
from pvlib.forecast import GFS, NAM, NDFD, RAP, HRRR


# ## Load Forecast data

# pvlib forecast module only includes several models. To see the full list of forecast models visit the Unidata website:
# 
# http://www.unidata.ucar.edu/data/#tds

# In[2]:


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


# Define some PV system parameters.

# In[3]:


surface_tilt = 30
surface_azimuth = 180 # pvlib uses 0=North, 90=East, 180=South, 270=West convention
albedo = 0.2


# In[4]:


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


# In[5]:


# Define forecast model
fm = GFS()
#fm = NAM()
#fm = NDFD()
#fm = RAP()
#fm = HRRR()


# In[6]:


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


# Let's look at the downloaded version of the forecast data.

# In[7]:


forecast_data.head()


# This is a ``pandas DataFrame`` object. It has a lot of great properties that are beyond the scope of our tutorials.

# In[8]:


forecast_data['temp_air'].plot()


# Plot the GHI data. Most pvlib forecast models derive this data from the weather models' cloud clover data.

# In[9]:


ghi = forecast_data['ghi']
ghi.plot()
plt.ylabel('Irradiance ($W/m^{-2}$)')


# ## Calculate modeling intermediates

# Before we can calculate power for all the forecast times, we will need to calculate:
# * solar position 
# * extra terrestrial radiation
# * airmass
# * angle of incidence
# * POA sky and ground diffuse radiation
# * cell and module temperatures
# 
# The approach here follows that of the pvlib tmy_to_power notebook. You will find more details regarding this approach and the values being calculated in that notebook.

# ### Solar position

# Calculate the solar position for all times in the forecast data. 
# 
# The default solar position algorithm is based on Reda and Andreas (2004). Our implementation is pretty fast, but you can make it even faster if you install [``numba``](http://numba.pydata.org/#installing) and use add  ``method='nrel_numba'`` to the function call below.

# In[10]:


# retrieve time and location parameters
time = forecast_data.index
a_point = fm.location


# In[11]:


solpos = a_point.get_solarposition(time)
#solpos.plot()


# The funny looking jump in the azimuth is just due to the coarse time sampling in the TMY file.

# ### DNI ET
# 
# Calculate extra terrestrial radiation. This is needed for many plane of array diffuse irradiance models.

# In[12]:


dni_extra = irradiance.extraradiation(fm.time)

#dni_extra.plot()
#plt.ylabel('Extra terrestrial radiation ($W/m^{-2}$)')


# ### Airmass
# 
# Calculate airmass. Lots of model options here, see the ``atmosphere`` module tutorial for more details.

# In[13]:


airmass = atmosphere.relativeairmass(solpos['apparent_zenith'])

#airmass.plot()
#plt.ylabel('Airmass')


# The funny appearance is due to aliasing and setting invalid numbers equal to ``NaN``. Replot just a day or two and you'll see that the numbers are right.

# ### POA sky diffuse

# Use the Hay Davies model to calculate the plane of array diffuse sky radiation. See the ``irradiance`` module tutorial for comparisons of different models.

# In[14]:


poa_sky_diffuse = irradiance.haydavies(surface_tilt, surface_azimuth,
                                             forecast_data['dhi'], forecast_data['dni'], dni_extra,
                                             solpos['apparent_zenith'], solpos['azimuth'])
#poa_sky_diffuse.plot()
#plt.ylabel('Irradiance ($W/m^{-2}$)')


# ### POA ground diffuse
# 
# Calculate ground diffuse. We specified the albedo above. You could have also provided a string to the ``surface_type`` keyword argument.

# In[15]:


poa_ground_diffuse = irradiance.grounddiffuse(surface_tilt, ghi, albedo=albedo)

#poa_ground_diffuse.plot()
#plt.ylabel('Irradiance ($W/m^{-2}$)')


# ### AOI
# 
# Calculate AOI

# In[16]:


aoi = irradiance.aoi(surface_tilt, surface_azimuth, solpos['apparent_zenith'], solpos['azimuth'])

#aoi.plot()
#plt.ylabel('Angle of incidence (deg)')


# Note that AOI has values greater than 90 deg. This is ok.

# ### POA total
# 
# Calculate POA irradiance

# In[17]:


poa_irrad = irradiance.globalinplane(aoi, forecast_data['dni'], poa_sky_diffuse, poa_ground_diffuse)

poa_irrad.plot()
plt.ylabel('Irradiance ($W/m^{-2}$)')
plt.title('POA Irradiance')


# ### Cell and module temperature
# 
# Calculate pv cell and module temperature

# In[18]:


temperature = forecast_data['temp_air']
wnd_spd = forecast_data['wind_speed']
pvtemps = pvsystem.sapm_celltemp(poa_irrad['poa_global'], wnd_spd, temperature)

pvtemps.plot()
plt.ylabel('Temperature (C)')


# ## DC power using SAPM

# Get module data from the web.

# In[19]:


sandia_modules = pvsystem.retrieve_sam('SandiaMod')


# Choose a particular module

# In[20]:


sandia_module = sandia_modules.Canadian_Solar_CS5P_220M___2009_
sandia_module


# Run the SAPM using the parameters we calculated above.

# In[21]:


effective_irradiance = pvsystem.sapm_effective_irradiance(poa_irrad.poa_direct, poa_irrad.poa_diffuse, 
                                                          airmass, aoi, sandia_module)

sapm_out = pvsystem.sapm(effective_irradiance, pvtemps['temp_cell'], sandia_module)
#print(sapm_out.head())

sapm_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')


# ## AC power using SAPM

# Get the inverter database from the web

# In[22]:


sapm_inverters = pvsystem.retrieve_sam('sandiainverter')


# Choose a particular inverter

# In[23]:


sapm_inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_']
sapm_inverter


# In[24]:


p_ac = pvsystem.snlinverter(sapm_out.v_mp, sapm_out.p_mp, sapm_inverter)

p_ac.plot()
plt.ylabel('AC Power (W)')
plt.ylim(0, None)


# Plot just a few days.

# In[25]:


p_ac[start:start+pd.Timedelta(days=2)].plot()


# Some statistics on the AC power

# In[26]:


p_ac.describe()


# In[27]:


p_ac.index.freq


# In[28]:


# integrate power to find energy yield over the forecast period
p_ac.sum() * 3


# In[ ]: