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
  2. Load Forecast data
  3. Calculate modeling intermediates
  4. DC power using SAPM
  5. 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
%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


/Users/holmgren/git_repos/pvlib-python/pvlib/forecast.py:22: UserWarning: The forecast module algorithms and features are highly experimental. The API may change, the functionality may be consolidated into an io module, or the module may be separated into its own package.
  'module, or the module may be separated into its own package.')

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


Out[7]:
temp_air wind_speed ghi dni dhi total_clouds low_clouds mid_clouds high_clouds
2016-07-27 06:00:00-06:00 30.149994 1.916377 0.000000 0.000000 0.000000 7.0 0.0 3.0 5.0
2016-07-27 09:00:00-06:00 28.450012 0.990202 374.978175 468.494996 150.358998 20.0 0.0 7.0 17.0
2016-07-27 12:00:00-06:00 27.450012 3.230511 765.140377 427.065650 375.447762 25.0 0.0 6.0 22.0
2016-07-27 15:00:00-06:00 32.350006 1.501533 325.739565 18.348873 308.994555 99.0 0.0 77.0 96.0
2016-07-27 18:00:00-06:00 46.050018 4.222464 240.470023 87.506983 198.519719 68.0 0.0 44.0 64.0

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


Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x11612c048>

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


Out[9]:
<matplotlib.text.Text at 0x115b33550>

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


Out[17]:
<matplotlib.text.Text at 0x11f43acf8>

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


Out[18]:
<matplotlib.text.Text at 0x11fd53828>

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


Out[20]:
Vintage                                                          2009
Area                                                            1.701
Material                                                         c-Si
Cells_in_Series                                                    96
Parallel_Strings                                                    1
Isco                                                          5.09115
Voco                                                          59.2608
Impo                                                          4.54629
Vmpo                                                          48.3156
Aisc                                                         0.000397
Aimp                                                         0.000181
C0                                                            1.01284
C1                                                         -0.0128398
Bvoco                                                        -0.21696
Mbvoc                                                               0
Bvmpo                                                       -0.235488
Mbvmp                                                               0
N                                                              1.4032
C2                                                           0.279317
C3                                                           -7.24463
A0                                                           0.928385
A1                                                           0.068093
A2                                                         -0.0157738
A3                                                          0.0016606
A4                                                          -6.93e-05
B0                                                                  1
B1                                                          -0.002438
B2                                                          0.0003103
B3                                                         -1.246e-05
B4                                                           2.11e-07
B5                                                          -1.36e-09
DTC                                                                 3
FD                                                                  1
A                                                            -3.40641
B                                                          -0.0842075
C4                                                           0.996446
C5                                                           0.003554
IXO                                                           4.97599
IXXO                                                          3.18803
C6                                                            1.15535
C7                                                          -0.155353
Notes               Source: Sandia National Laboratories Updated 9...
Name: Canadian_Solar_CS5P_220M___2009_, dtype: object

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


Out[21]:
<matplotlib.text.Text at 0x1202e6c50>

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


Out[23]:
Vac          208.000000
Paco         250.000000
Pdco         259.522050
Vdco          40.242603
Pso            1.771614
C0            -0.000025
C1            -0.000090
C2             0.000669
C3            -0.018900
Pnt            0.020000
Vdcmax        65.000000
Idcmax        10.000000
Mppt_low      20.000000
Mppt_high     50.000000
Name: ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_, dtype: float64

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)


Out[24]:
(0, 200.0)

Plot just a few days.


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


Out[25]:
<matplotlib.axes._subplots.AxesSubplot at 0x12032def0>

Some statistics on the AC power


In [26]:
p_ac.describe()


Out[26]:
count     57.000000
mean      29.373181
std       40.204146
min       -0.020000
25%       -0.020000
50%       -0.020000
75%       54.245105
max      161.882005
dtype: float64

In [27]:
p_ac.index.freq


Out[27]:
<3 * Hours>

In [28]:
# integrate power to find energy yield over the forecast period
p_ac.sum() * 3


Out[28]:
5022.81403518355

In [ ]: