This tutorial will walk through the process of going from Unidata forecast model data to AC power using the SAPM.
Table of contents:
This tutorial has been tested against the following package versions:
It should work with other Python and Pandas versions. It requires pvlib >= 0.2.0 and IPython >= 3.0.
Authors:
These are just your standard interactive scientific python imports that you'll get very used to using.
In [1]:
# built-in python modules
from datetime import datetime,timedelta
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
%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
pvlib forecast module only includes several models. To see the full list of forecast models visit the Unidata website:
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 = datetime.now() # today's date
end = start + timedelta(days=7) # 7 days from today
timerange = pd.date_range(start, end, tz=tz)
In [5]:
# Define forecast model
fm = GFS()
#fm = NAM()
#fm = NDFD()
#fm = RAP()
#fm = HRRR()
In [6]:
# Retrieve data
forecast_data = fm.get_query_data(latitude, longitude, timerange)
Let's look at the downloaded version of the forecast data.
In [7]:
forecast_data.head()
Out[7]:
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['temperature'].plot()
Out[8]:
Plot the GHI data
In [9]:
ghi = forecast_data['ghi']
ghi.plot()
plt.ylabel('Irradiance ($W/m^{-2}$)')
Out[9]:
Before we can calculate power for all the forecast times, we will need to calculate:
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.
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 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 = solarposition.get_solarposition(time, a_point)
#solpos.plot()
The funny looking jump in the azimuth is just due to the coarse time sampling in the TMY file.
In [12]:
dni_extra = irradiance.extraradiation(fm.time)
dni_extra = pd.Series(dni_extra, index=fm.time)
#dni_extra.plot()
#plt.ylabel('Extra terrestrial radiation ($W/m^{-2}$)')
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.
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}$)')
In [15]:
poa_ground_diffuse = irradiance.grounddiffuse(surface_tilt, ghi, albedo=albedo)
#poa_ground_diffuse.plot()
#plt.ylabel('Irradiance ($W/m^{-2}$)')
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.
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')
Out[17]:
In [18]:
temperature = forecast_data['temperature']
wnd_spd = forecast_data['wind_speed']
pvtemps = pvsystem.sapm_celltemp(poa_irrad['poa_global'], wnd_spd, temperature)
pvtemps.plot()
plt.ylabel('Temperature (C)')
Out[18]:
Get module data from the web.
In [19]:
sandia_modules = pvsystem.retrieve_sam(name='SandiaMod')
Choose a particular module
In [20]:
sandia_module = sandia_modules.Canadian_Solar_CS5P_220M___2009_
sandia_module
Out[20]:
Run the SAPM using the parameters we calculated above.
In [21]:
sapm_out = pvsystem.sapm(sandia_module, poa_irrad.poa_direct, poa_irrad.poa_diffuse,
pvtemps['temp_cell'], airmass, aoi)
#print(sapm_out.head())
sapm_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')
Out[21]:
In [22]:
cec_modules = pvsystem.retrieve_sam(name='CECMod')
cec_module = cec_modules.Canadian_Solar_CS5P_220M
In [23]:
photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth = (
pvsystem.calcparams_desoto(poa_irrad.poa_global,
temp_cell=pvtemps['temp_cell'],
alpha_isc=cec_module['alpha_sc'],
module_parameters=cec_module,
EgRef=1.121,
dEgdT=-0.0002677) )
In [24]:
single_diode_out = pvsystem.singlediode(cec_module, photocurrent, saturation_current,
resistance_series, resistance_shunt, nNsVth)
In [25]:
single_diode_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')
Out[25]:
Get the inverter database from the web
In [26]:
sapm_inverters = pvsystem.retrieve_sam('sandiainverter')
Choose a particular inverter
In [27]:
sapm_inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_']
sapm_inverter
Out[27]:
In [28]:
p_acs = pd.DataFrame()
p_acs['sapm'] = pvsystem.snlinverter(sapm_inverter, sapm_out.v_mp, sapm_out.p_mp)
p_acs['sd'] = pvsystem.snlinverter(sapm_inverter, single_diode_out.v_mp, single_diode_out.p_mp)
p_acs.plot()
plt.ylabel('AC Power (W)')
Out[28]:
In [29]:
#diff = p_acs['sapm'] - p_acs['sd']
#diff.plot()
#plt.ylabel('SAPM - SD Power (W)')
Plot just a few days.
In [30]:
p_acs[start:start+timedelta(days=2)].plot()
Out[30]:
Some statistics on the AC power
In [31]:
p_acs.describe()
Out[31]:
In [32]:
p_acs.sum()
Out[32]: