Stephen Lightfoote, blogdown.rmd=TRUE

7 minute read


So, alot is made of the functionality of PVSYST as a standard tool to model PV systems. No arguments here. That said, it does not have a programming interface. Thankfully there are some great alternatives out there developed through Sandia National Laboratories, particularly the PV-LIB Toolbox (, which contains some really excellent functions for modeling solar position, irradiance, and PV systems. The tool was initially created in MATLAB and subsequently migrated to Python (pvlib:

Goals and Objectives

Model Plane of Array Irradiance for a fictitious fixed tilt and tracking PV system scenarios using pvlib scripted through R via the reticulate package. I realize it would probably be a bit smoother to use Python directly, but to be honest I love R and also wanted to test out a hybrid language workflow using a combination of R and Python available through reticulate package Lastly, I wanted to see how much solar irradiance I could expect from an array mounted on my roof at home.

Let’s walk through the steps and see how this all works.

Load R Libraries

First, let’s load the r libraries and define a timestamp convention.

library(tidyverse) # data analysis and visualization
library(reticulate) # loading python libraries in R
library(dygraphs) # timeseries visualization
library(timetk) # tidy timeseries tools
library(leaflet) # mapping
library(knitr) # tables
library(viridis) # color sheme
library(ggridges) # distribution plot

load NSRDB data

To do this, we’ll be querying data from the National Renewable Energy Laboratory (NREL) Physical Solar Model v3 via their developer API ( I’m going to use my own home in MA as an example.


leaflet(data = data.frame(lat=lat,lon=lon,name='home')) %>% 
  addProviderTiles('Esri.WorldImagery') %>%
  addMarkers(lng = ~lon,lat=~lat) %>%

query the data

You’ll need your own API key to run this, which you can get for free through NREL’s website For this example I’m just going to grab a random year of data.

# information on dataset: 
attributes = c('air_temperature,clearsky_dhi,clearsky_dni,clearsky_ghi,cloud_type,dew_point,dhi,dni,fill_flag,ghi,relative_humidity,solar_zenith_angle,surface_albedo,surface_pressure,total_precipitable_water,wind_direction,wind_speed')
leap_year = 'false'

url = paste('',lon,'%20',lat,')&names=',year,'&leap_day=false&interval=',interval,'&utc=',utc,
#download data as a dataframe
df<-read_csv(url,skip=2,col_types = cols())
# tidy names
names(df)<-gsub(pattern = ' ',replacement = '',x = names(df))

here’s what the data looks like.

## # A tibble: 6 x 23
##    Year Month   Day  Hour Minute Temperature ClearskyDHI ClearskyDNI
##   <int> <int> <int> <int>  <int>       <int>       <int>       <int>
## 1  2010     1     1     0     30          -3           0           0
## 2  2010     1     1     1     30          -3           0           0
## 3  2010     1     1     2     30          -3           0           0
## 4  2010     1     1     3     30          -3           0           0
## 5  2010     1     1     4     30          -3           0           0
## 6  2010     1     1     5     30          -3           0           0
## # ... with 15 more variables: ClearskyGHI <int>, CloudType <int>,
## #   DewPoint <int>, DHI <int>, DNI <int>, FillFlag <int>, GHI <int>,
## #   RelativeHumidity <dbl>, SolarZenithAngle <dbl>, SurfaceAlbedo <dbl>,
## #   Pressure <dbl>, PrecipitableWater <dbl>, WindDirection <dbl>,
## #   WindSpeed <dbl>, timestamp_utc <dttm>

Load Python Libraries

Next, I’ll load the python libraries that I’ll need for this example using reticulate, including pvlib. Note, the “convert=F” option for PVLIB in which I’d like to keep the output from each pvlib function in it’s native python format rather than auto-converting it to its R equivalent.


PV system setup

I’ve got pretty decent sun exposure on my roof (170 degrees azimuth, so almost due south) and my roof tilt is almost in lock stop with my latitude at 45 degrees from horizon. So let’s model a fixed tilt system on my roof. Also, since this is somewhat of a tutorial illustrating functionality, for kicks let’s model a tracking system that I’ll put on my front lawn. I’m sure the neighbors would LOVE it.

location = pvlib$location$Location(latitude = lat,longitude = lon,tz = 'UTC',altitude = 100,name = 'home')
pvsystem_fixed = pvlib$pvsystem$PVSystem(surface_tilt=45,surface_azimuth=170,albedo=.2,name='home fixed tilt')
pvsystem_tracking=pvlib$tracking$SingleAxisTracker(axis_tilt = 0,axis_azimuth = 170,max_angle = 52,backtrack = TRUE,gcr = .5)

python_timestamp<-pd$DatetimeIndex(start='2010-01-01 00:30', end='2010-12-31 23:30', freq='1h',tz=location$tz)

Solar Position and Atmospheric Modeling

To convert the three components of irradiance provided from NSRDB (GHI,DHI and DNI) to the plane of array (POA) of my rooftop, I first need to know where the sun is in the sky at any given point in the year, some default assumptions about the amount of extraterrestrial radiation and the atmosphere it is travelling through and ultimately the angle of incidence (AOI) on the collector plane of my rooftop solar modules.

python_timestamp<-pd$DatetimeIndex(start=as.character(min(df$timestamp_utc)),end=as.character(max(df$timestamp_utc)), freq='1h',tz=location$tz)
solpos = location$get_solarposition(python_timestamp)
dni_extra = pvlib$irradiance$extraradiation(python_timestamp)
airmass = pvlib$atmosphere$relativeairmass(solpos$apparent_zenith)
aoi = pvlib$irradiance$aoi(pvsystem_fixed$surface_tilt,pvsystem_fixed$surface_azimuth,solpos$apparent_zenith, solpos$azimuth)

# convert output back to R-equivalent and append to data frame.

Plane of Array Transposition

Calculating POA from GHI, DNI and DHI is known as transposition and standard transposition model is the ‘perez’ model which thankfully pvlib has a handy function for. R. Perez, P. Ineichen, R. Seals, and J. Michalsky, “Modeling daylight availability and irradiance components from direct and global irradiance,” Sol. Energy, vol. 44, pp. 271-289, 1990.

Fixed Tilt Example

We previously calculated the AOI on the module surface for all periods in our dataset as well as other relevant parameters, so we have everything we need to model POA for the fixed tilt scenario.

total_irrad = pvlib$irradiance$total_irrad(surface_tilt=pvsystem_fixed$surface_tilt,
                                           apparent_zenith= solpos$apparent_zenith,
                                           azimuth= solpos$azimuth,
                                           dni= pd$Series(df$DNI,python_timestamp), # pvlib requires a data series with an index
                                           ghi= pd$Series(df$GHI,python_timestamp),
                                           dhi= pd$Series(df$DHI,python_timestamp),
                                           dni_extra= dni_extra,
                                           solar_zenith= solpos$zenith,
                                           airmass= airmass)

df<-replace_na(df,list(poa_global=0,poa_direct=0,poa_diffuse=0,poa_sky_diffuse=0,poa_ground_diffuse=0)) %>%

Single Axis Tracking Example

The single axis tracking example is a bit more complicated than the fixed tilt as we’ll first need to model how the tracker follows the sun such that we can know what the AOI will be on the trackers at any given time. But once we have that, the

tracker_data = pvlib$tracking$singleaxis(solpos$apparent_zenith,
                                         backtrack=TRUE, gcr=.5)

perez_diffuse = pvlib$irradiance$perez(tracker_data$surface_tilt,
                                       pd$Series(df$DNI, python_timestamp),

ground_irrad = pvlib$irradiance$grounddiffuse(tracker_data$surface_tilt, pd$Series(df$GHI,python_timestamp), albedo=.2)
poa_global_tracking = py_to_r(cosd)*df$DNI + py_to_r(perez_diffuse) + py_to_r(ground_irrad)


Results and Comparisons

All data

# convert to local time

df %>%
  select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI,DNI,DHI) %>% 
  tk_xts(date_var=timestamp_local,silent = T) %>%
  dygraph() %>% 
  dyRoller() %>% 
  dyRangeSelector() %>% 
  dyOptions(useDataTimezone = TRUE)

Example Clear Day

You can see on a clear day that the POA irradiance at my latitude in the northern hemisphere is typically a bit higher than it’s GHI equivalent. You can also see a slight skew with higher POA values in the morning for the fixed tilt scenario than for the afternoon and vice versa for the tracker scenario owed to thenon-ideal azimuthal orientation of my home to due south (170 degrees, instead of 180).

df %>%
  filter(timestamp_local>=as.POSIXct('2010-09-21',tz=local_tzone) & timestamp_local<as.POSIXct('2010-09-22',tz=local_tzone)) %>%
  select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI,DNI,DHI) %>%
  gather(type,`W/m^2`,poa_global_fixed:DHI) %>%
  scale_color_viridis(discrete = T)+
  labs(title='Fixed Tilt and Tracking Sample Clear Day')

Total Irradiance

aggregates<-df %>%
  group_by(Month) %>%

  kable(aggregates,digits = 1,align = 'c',caption='Monthly kWh/m^2')
Table 1: Monthly kWh/m^2
Month GHI POA_fixed POA_tracking
1 62.7 125.3 83.6
2 70.3 116.7 92.3
3 101.5 132.6 132.1
4 148.1 159.4 187.5
5 198.8 186.4 250.7
6 179.2 158.7 217.8
7 207.4 187.5 263.5
8 170.0 173.5 215.3
9 130.1 156.4 163.9
10 94.2 141.2 122.7
11 59.2 107.9 77.6
12 52.8 111.1 70.5
  kable(colSums(select(aggregates,-Month)),digits=1,align='c',caption = 'Annual kWh/m^2')
Table 1: Annual kWh/m^2
GHI 1474.3
POA_fixed 1756.7
POA_tracking 1877.5

Regardless of my PV system choice (fixed tilt vs. tracking) I would gain extra insolation over the course of the year as compared to just placing the panels on the ground (GHI). Looking at a distribution of irradiance by month and mounting type, it becomes clear that the tracker mount would pay the most dividends over the summer months (higher irradiance values) and the fixed tilt mounting would win out in the winter months with the crossover point occuring during spring/fall.

df %>%
  select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI) %>%
  gather(type,`W/m^2`,poa_global_fixed:GHI) %>%
  mutate(month=factor(format(timestamp_local,'%m'))) %>%
  filter(`W/m^2`>0) %>%
  theme_ridges(font_size = 12)+
  scale_fill_viridis(discrete = T)+
  theme(legend.position = 'bottom')+
  labs(title='Irradiance Distribution by Month',caption='Filtered for W/m^2 > 0')


OK, so while ultimately mounting my fictitious PV system on a tracker rack would give me ~6% more POA insolation over the course of the year, the juice is probably not worth the squeeze and I should stick with just mounting the panels on my roof like everyone else. In the next iteration I’ll try and walk through modeling energy capture of my fictitous fixed tilt and tracking PV systems in pvlib.

comments powered by Disqus