Clear sky¶
This section reviews the clear sky modeling capabilities of pvlib-python.
pvlib-python supports two ways to generate clear sky irradiance:
- A
Location
object’sget_clearsky()
method. - The functions contained in the
clearsky
module, includingineichen()
andsimplified_solis()
.
Users that work with simple time series data may prefer to use
get_clearsky()
, while users
that want finer control, more explicit code, or work with
multidimensional data may prefer to use the basic functions in the
clearsky
module.
The Location subsection demonstrates the easiest way to obtain a time series of clear sky data for a location. The Ineichen and Perez and Simplified Solis subsections detail the clear sky algorithms and input data. The Detect Clearsky subsection demonstrates the use of the clear sky detection algorithm.
The bird_hulstrom80_aod_bb()
, and
kasten96_lt()
functions are useful for
calculating inputs to the clear sky functions.
We’ll need these imports for the examples below.
In [1]: import itertools
In [2]: import matplotlib.pyplot as plt
In [3]: import pandas as pd
# seaborn makes the plots look nicer
In [4]: import seaborn as sns
In [5]: sns.set_color_codes()
In [6]: import pvlib
In [7]: from pvlib import clearsky, atmosphere
In [8]: from pvlib.location import Location
Location¶
The easiest way to obtain a time series of clear sky irradiance is to use a
Location
object’s
get_clearsky()
method. The
get_clearsky()
method does the dirty
work of calculating solar position, extraterrestrial irradiance,
airmass, and atmospheric pressure, as appropriate, leaving the user to
only specify the most important parameters: time and atmospheric
attenuation. The time input must be a pandas.DatetimeIndex
,
while the atmospheric attenuation inputs may be constants or arrays.
The get_clearsky()
method always
returns a pandas.DataFrame
.
In [9]: tus = Location(32.2, -111, 'US/Arizona', 700, 'Tucson')
In [10]: times = pd.DatetimeIndex(start='2016-07-01', end='2016-07-04', freq='1min', tz=tus.tz)
In [11]: cs = tus.get_clearsky(times) # ineichen with climatology table by default
In [12]: cs.plot();
In [13]: plt.ylabel('Irradiance $W/m^2$');
In [14]: plt.title('Ineichen, climatological turbidity');
The get_clearsky()
method accepts a
model keyword argument and propagates additional arguments to the
functions that do the computation.
In [15]: cs = tus.get_clearsky(times, model='ineichen', linke_turbidity=3)
In [16]: cs.plot();
In [17]: plt.title('Ineichen, linke_turbidity=3');
In [18]: plt.ylabel('Irradiance $W/m^2$');
In [19]: cs = tus.get_clearsky(times, model='simplified_solis', aod700=0.2, precipitable_water=3)
In [20]: cs.plot();
In [21]: plt.title('Simplfied Solis, aod700=0.2, precipitable_water=3');
In [22]: plt.ylabel('Irradiance $W/m^2$');
See the sections below for more detail on the clear sky models.
Ineichen and Perez¶
The Ineichen and Perez clear sky model parameterizes irradiance in terms
of the Linke turbidity [Ine02]. pvlib-python implements this model in
the pvlib.clearsky.ineichen()
function.
Turbidity data¶
pvlib includes a file with monthly climatological turbidity values for the globe. The code below creates turbidity maps for a few months of the year. You could run it in a loop to create plots for all months.
In [23]: import calendar
In [24]: import os
In [25]: import scipy.io
In [26]: pvlib_path = os.path.dirname(os.path.abspath(pvlib.clearsky.__file__))
In [27]: filepath = os.path.join(pvlib_path, 'data', 'LinkeTurbidities.mat')
In [28]: mat = scipy.io.loadmat(filepath)
# data is in units of 20 x turbidity
In [29]: linke_turbidity_table = mat['LinkeTurbidity'] # / 20. # crashes on rtd
In [30]: def plot_turbidity_map(month, vmin=1, vmax=100):
....: plt.figure();
....: plt.imshow(linke_turbidity_table[:, :, month-1], vmin=vmin, vmax=vmax);
....: plt.title('Linke turbidity x 20, ' + calendar.month_name[month]);
....: plt.colorbar(shrink=0.5);
....: plt.tight_layout();
....:
In [31]: plot_turbidity_map(1)
In [32]: plot_turbidity_map(7)
The lookup_linke_turbidity()
function takes a
time, latitude, and longitude and gets the corresponding climatological
turbidity value for that time at those coordinates. By default, the
lookup_linke_turbidity()
function will linearly
interpolate turbidity from month to month, assuming that the raw data is
valid on 15th of each month. This interpolation removes discontinuities
in multi-month PV models. Here’s a plot of a few locations in the
Southwest U.S. with and without interpolation. We chose points that are
relatively close so that you can get a better sense of the spatial noise
and variability of the data set. Note that the altitude of these sites
varies from 300 m to 1500 m.
In [33]: times = pd.DatetimeIndex(start='2015-01-01', end='2016-01-01', freq='1D')
In [34]: sites = [(32, -111, 'Tucson1'), (32.2, -110.9, 'Tucson2'),
....: (33.5, -112.1, 'Phoenix'), (35.1, -106.6, 'Albuquerque')]
....:
In [35]: plt.figure();
In [36]: for lat, lon, name in sites:
....: turbidity = pvlib.clearsky.lookup_linke_turbidity(times, lat, lon, interp_turbidity=False)
....: turbidity.plot(label=name)
....:
In [37]: plt.legend();
In [38]: plt.title('Raw data (no interpolation)');
In [39]: plt.ylabel('Linke Turbidity');
In [40]: plt.figure();
In [41]: for lat, lon, name in sites:
....: turbidity = pvlib.clearsky.lookup_linke_turbidity(times, lat, lon)
....: turbidity.plot(label=name)
....:
In [42]: plt.legend();
In [43]: plt.title('Interpolated to the day');
In [44]: plt.ylabel('Linke Turbidity');
Examples¶
A clear sky time series using only basic pvlib functions.
In [45]: latitude, longitude, tz, altitude, name = 32.2, -111, 'US/Arizona', 700, 'Tucson'
In [46]: times = pd.date_range(start='2014-01-01', end='2014-01-02', freq='1Min', tz=tz)
In [47]: solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
In [48]: apparent_zenith = solpos['apparent_zenith']
In [49]: airmass = pvlib.atmosphere.relativeairmass(apparent_zenith)
In [50]: pressure = pvlib.atmosphere.alt2pres(altitude)
In [51]: airmass = pvlib.atmosphere.absoluteairmass(airmass, pressure)
In [52]: linke_turbidity = pvlib.clearsky.lookup_linke_turbidity(times, latitude, longitude)
In [53]: dni_extra = pvlib.irradiance.extraradiation(times.dayofyear)
# an input is a pandas Series, so solis is a DataFrame
In [54]: ineichen = clearsky.ineichen(apparent_zenith, airmass, linke_turbidity, altitude, dni_extra)
In [55]: plt.figure();
In [56]: ax = ineichen.plot()
In [57]: ax.set_ylabel('Irradiance $W/m^2$');
In [58]: ax.set_title('Ineichen Clear Sky Model');
In [59]: ax.legend(loc=2);
In [60]: plt.show();
The input data types determine the returned output type. Array input results in an OrderedDict of array output, and Series input results in a DataFrame output. The keys are ‘ghi’, ‘dni’, and ‘dhi’.
Grid with a clear sky irradiance for a few turbidity values.
In [61]: times = pd.date_range(start='2014-09-01', end='2014-09-02', freq='1Min', tz=tz)
In [62]: solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
In [63]: apparent_zenith = solpos['apparent_zenith']
In [64]: airmass = pvlib.atmosphere.relativeairmass(apparent_zenith)
In [65]: pressure = pvlib.atmosphere.alt2pres(altitude)
In [66]: airmass = pvlib.atmosphere.absoluteairmass(airmass, pressure)
In [67]: linke_turbidity = pvlib.clearsky.lookup_linke_turbidity(times, latitude, longitude)
In [68]: print('climatological linke_turbidity = {}'.format(linke_turbidity.mean()))
climatological linke_turbidity = 3.329496820286459
In [69]: dni_extra = pvlib.irradiance.extraradiation(times.dayofyear)
In [70]: linke_turbidities = [linke_turbidity.mean(), 2, 4]
In [71]: fig, axes = plt.subplots(ncols=3, nrows=1, sharex=True, sharey=True, squeeze=True, figsize=(12, 4))
In [72]: axes = axes.flatten()
In [73]: for linke_turbidity, ax in zip(linke_turbidities, axes):
....: ineichen = clearsky.ineichen(apparent_zenith, airmass, linke_turbidity, altitude, dni_extra)
....: ineichen.plot(ax=ax, title='Linke turbidity = {:0.1f}'.format(linke_turbidity));
....:
In [74]: ax.legend(loc=1);
In [75]: plt.show();
Validation¶
Will Holmgren compared pvlib’s Ineichen model and climatological turbidity to SoDa’s McClear service in Arizona. Here are links to an ipynb notebook and its html rendering.
Simplified Solis¶
The Simplified Solis model parameterizes irradiance in terms of
precipitable water and aerosol optical depth [Ine08ss]. pvlib-python
implements this model in the pvlib.clearsky.simplified_solis()
function.
Aerosol and precipitable water data¶
There are a number of sources for aerosol and precipitable water data
of varying accuracy, global coverage, and temporal resolution.
Ground based aerosol data can be obtained from
Aeronet. Precipitable water can be obtained
from radiosondes,
ESRL GPS-MET, or
derived from surface relative humidity using functions such as
pvlib.atmosphere.gueymard94_pw()
.
Numerous gridded products from satellites, weather models, and climate models
contain one or both of aerosols and precipitable water. Consider data
from the ECMWF
and SoDa.
Aerosol optical depth is a function of wavelength, and the Simplified Solis model requires AOD at 700 nm. Models exist to convert AOD between different wavelengths, as well as convert Linke turbidity to AOD and PW [Ine08con], [Ine16].
Examples¶
A clear sky time series using only basic pvlib functions.
In [76]: latitude, longitude, tz, altitude, name = 32.2, -111, 'US/Arizona', 700, 'Tucson'
In [77]: times = pd.date_range(start='2014-01-01', end='2014-01-02', freq='1Min', tz=tz)
In [78]: solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
In [79]: apparent_elevation = solpos['apparent_elevation']
In [80]: aod700 = 0.1
In [81]: precipitable_water = 1
In [82]: pressure = pvlib.atmosphere.alt2pres(altitude)
In [83]: dni_extra = pvlib.irradiance.extraradiation(times.dayofyear)
# an input is a Series, so solis is a DataFrame
In [84]: solis = clearsky.simplified_solis(apparent_elevation, aod700, precipitable_water,
....: pressure, dni_extra)
....:
In [85]: ax = solis.plot();
In [86]: ax.set_ylabel('Irradiance $W/m^2$');
In [87]: ax.set_title('Simplified Solis Clear Sky Model');
In [88]: ax.legend(loc=2);
In [89]: plt.show();
The input data types determine the returned output type. Array input results in an OrderedDict of array output, and Series input results in a DataFrame output. The keys are ‘ghi’, ‘dni’, and ‘dhi’.
Irradiance as a function of solar elevation.
In [90]: apparent_elevation = pd.Series(np.linspace(-10, 90, 101))
In [91]: aod700 = 0.1
In [92]: precipitable_water = 1
In [93]: pressure = 101325
In [94]: dni_extra = 1364
In [95]: solis = clearsky.simplified_solis(apparent_elevation, aod700,
....: precipitable_water, pressure, dni_extra)
....:
In [96]: ax = solis.plot();
In [97]: ax.set_xlabel('Apparent elevation (deg)');
In [98]: ax.set_ylabel('Irradiance $W/m^2$');
In [99]: ax.set_title('Irradiance vs Solar Elevation')
Out[99]: <matplotlib.text.Text at 0x7f83d176a898>
In [100]: ax.legend(loc=2);
Grid with a clear sky irradiance for a few PW and AOD values.
In [101]: times = pd.date_range(start='2014-09-01', end='2014-09-02', freq='1Min', tz=tz)
In [102]: solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
In [103]: apparent_elevation = solpos['apparent_elevation']
In [104]: pressure = pvlib.atmosphere.alt2pres(altitude)
In [105]: dni_extra = pvlib.irradiance.extraradiation(times.dayofyear)
In [106]: aod700 = [0.01, 0.1]
In [107]: precipitable_water = [0.5, 5]
In [108]: fig, axes = plt.subplots(ncols=2, nrows=2, sharex=True, sharey=True, squeeze=True)
In [109]: axes = axes.flatten()
In [110]: for (aod, pw), ax in zip(itertools.chain(itertools.product(aod700, precipitable_water)), axes):
.....: cs = clearsky.simplified_solis(apparent_elevation, aod, pw, pressure, dni_extra)
.....: cs.plot(ax=ax, title='aod700={}, pw={}'.format(aod, pw))
.....:
In [111]: plt.show();
Contour plots of irradiance as a function of both PW and AOD.
In [112]: aod700 = np.linspace(0, 0.5, 101)
In [113]: precipitable_water = np.linspace(0, 10, 101)
In [114]: apparent_elevation = 70
In [115]: pressure = 101325
In [116]: dni_extra = 1364
In [117]: aod700, precipitable_water = np.meshgrid(aod700, precipitable_water)
# inputs are arrays, so solis is an OrderedDict
In [118]: solis = clearsky.simplified_solis(apparent_elevation, aod700,
.....: precipitable_water, pressure,
.....: dni_extra)
.....:
In [119]: cmap = plt.get_cmap('viridis')
In [120]: n = 15
In [121]: vmin = None
In [122]: vmax = None
In [123]: def plot_solis(key):
.....: irrad = solis[key]
.....: fig, ax = plt.subplots()
.....: im = ax.contour(aod700, precipitable_water, irrad[:, :], n, cmap=cmap, vmin=vmin, vmax=vmax)
.....: imf = ax.contourf(aod700, precipitable_water, irrad[:, :], n, cmap=cmap, vmin=vmin, vmax=vmax)
.....: ax.set_xlabel('AOD')
.....: ax.set_ylabel('Precipitable water (cm)')
.....: ax.clabel(im, colors='k', fmt='%.0f')
.....: fig.colorbar(imf, label='{} (W/m**2)'.format(key))
.....: ax.set_title('{}, elevation={}'.format(key, apparent_elevation))
.....:
In [124]: plot_solis('ghi')
In [125]: plt.show()
In [126]: plot_solis('dni')
In [127]: plt.show()
In [128]: plot_solis('dhi')
In [129]: plt.show()
Validation¶
See [Ine16].
We encourage users to compare the pvlib implementation to Ineichen’s Excel tool.
Detect Clearsky¶
The detect_clearsky()
function implements the
[Ren16] algorithm to detect the clear and cloudy points of a time
series. The algorithm was designed and validated for analyzing GHI time
series only. Users may attempt to apply it to other types of time series
data using different filter settings, but should be skeptical of the
results.
The algorithm detects clear sky times by comparing statistics for a measured time series and an expected clearsky time series. Statistics are calculated using a sliding time window (e.g., 10 minutes). An iterative algorithm identifies clear periods, uses the identified periods to estimate bias in the clearsky data, scales the clearsky data and repeats.
Clear times are identified by meeting 5 criteria. Default values for these thresholds are appropriate for 10 minute windows of 1 minute GHI data.
Next, we show a simple example of applying the algorithm to synthetic GHI data. We first generate and plot the clear sky and measured data.
In [130]: abq = Location(35.04, -106.62, altitude=1619)
In [131]: times = pd.DatetimeIndex(start='2012-04-01 10:30:00', tz='Etc/GMT+7', periods=30, freq='1min')
In [132]: cs = abq.get_clearsky(times)
# scale clear sky data to account for possibility of different turbidity
In [133]: ghi = cs['ghi']*.953
# add a cloud event
In [134]: ghi['2012-04-01 10:42:00':'2012-04-01 10:44:00'] = [500, 300, 400]
# add an overirradiance event
In [135]: ghi['2012-04-01 10:56:00'] = 950
In [136]: fig, ax = plt.subplots()
In [137]: ghi.plot(label='input');
In [138]: cs['ghi'].plot(label='ineichen clear');
In [139]: ax.set_ylabel('Irradiance $W/m^2$');
In [140]: plt.legend(loc=4);
In [141]: plt.show();
Now we run the synthetic data and clear sky estimate through the
detect_clearsky()
function.
In [142]: clear_samples = clearsky.detect_clearsky(ghi, cs['ghi'], cs.index, 10)
In [143]: fig, ax = plt.subplots()
In [144]: clear_samples.plot();
In [145]: ax.set_ylabel('Clear (1) or Cloudy (0)');
The algorithm detected the cloud event and the overirradiance event.
References¶
[Ine02] | (1, 2) P. Ineichen and R. Perez, “A New airmass independent formulation for the Linke turbidity coefficient”, Solar Energy, 73, pp. 151-157, 2002. |
[Ine08ss] | P. Ineichen, “A broadband simplified version of the Solis clear sky model,” Solar Energy, 82, 758-762 (2008). |
[Ine16] | (1, 2) P. Ineichen, “Validation of models that estimate the clear sky global and beam solar irradiance,” Solar Energy, 132, 332-344 (2016). |
[Ine08con] | P. Ineichen, “Conversion function between the Linke turbidity and the atmospheric water vapor and aerosol content”, Solar Energy, 82, 1095 (2008). |
[Ren12] | M. Reno, C. Hansen, and J. Stein, “Global Horizontal Irradiance Clear Sky Models: Implementation and Analysis”, Sandia National Laboratories, SAND2012-2389, 2012. |
[Ren16] | Reno, M.J. and C.W. Hansen, “Identification of periods of clear sky irradiance in time series of GHI measurements” Renewable Energy, v90, p. 520-531, 2016. |