location

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

query the data

You’ll need your own API key to run this, which you can get for free through NREL’s website https://nsrdb.nrel.gov/download-instructions. For this example I’m just going to grab a random year of data.

# information on dataset: https://nsrdb.nrel.gov/current-version 
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'

utc='true'
year<-2010
interval='60'
local_tzone<-'America/New_York'
url = paste('http://developer.nrel.gov/api/solar/nsrdb_psm3_download.csv?wkt=POINT(',lon,'%20',lat,')&names=',year,'&leap_day=false&interval=',interval,'&utc=',utc,
                '&full_name=',name,'&email=',email,'&affiliation=',affiliation,'&mailing_list=',mailing_list,'&reason=',reason,'&api_key=',api_key,'&attributes=',attributes,sep='')
    
#download data as a dataframe
df<-read_csv(url,skip=2,col_types = cols())
df$timestamp_utc=as.POSIXct(paste(df$Year,df$Month,df$Day,df$Hour,df$Minute,'00',sep='-'),format='%Y-%m-%d-%H-%M-%S',tz='UTC')
# tidy names
names(df)<-gsub(pattern = ' ',replacement = '',x = names(df))

here’s what the data looks like.

print(head(df))

## # 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.

pvlib<-import("pvlib",convert=F)
np<-import("numpy")
pd<-import("pandas")

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)

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,
                                           surface_azimuth=pvsystem_fixed$surface_azimuth,
                                           apparent_zenith= solpos$apparent_zenith,
                                           albedo=pvsystem_fixed$albedo,
                                           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,
                                           model='perez',
                                           solar_zenith= solpos$zenith,
                                           airmass= airmass)

df<-cbind(df,py_to_r(total_irrad))
df<-replace_na(df,list(poa_global=0,poa_direct=0,poa_diffuse=0,poa_sky_diffuse=0,poa_ground_diffuse=0)) %>%
  dplyr::rename(poa_global_fixed=poa_global)

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,
                                         solpos$azimuth,
                                         axis_tilt=pvsystem_tracking$axis_tilt,
                                         axis_azimuth=pvsystem_tracking$axis_azimuth,
                                         max_angle=52,
                                         backtrack=TRUE, gcr=.5)

perez_diffuse = pvlib$irradiance$perez(tracker_data$surface_tilt,
                                        tracker_data$surface_azimuth, 
                                       pd$Series(df$DHI,python_timestamp), 
                                       pd$Series(df$DNI, python_timestamp),
                                       dni_extra,
                                       solpos$apparent_zenith,
                                       solpos$azimuth,
                                       airmass=airmass)

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

df<-cbind(df,poa_global_tracking)%>%
  replace_na(list(poa_global_tracking=0))

All data

# convert to local time
df<-mutate(df,timestamp_local=timestamp_utc)
attributes(df$timestamp_local)$tzone<-local_tzone

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) %>%
  ggplot(aes(timestamp_local,`W/m^2`,color=type))+
  geom_line(lwd=1.2)+
  theme_minimal()+
  scale_color_viridis(discrete = T)+
  labs(title='Fixed Tilt and Tracking Sample Clear Day')

Total Irradiance

aggregates<-df %>%
  group_by(Month) %>%
  summarise(GHI=sum(GHI)*.001,
    POA_fixed=sum(poa_global_fixed)*.001,
            POA_tracking=sum(poa_global_tracking)*.001
  )

  kable(aggregates,digits = 1,align = 'c',caption='Monthly kWh/m^2')

  kable(colSums(select(aggregates,-Month)),digits=1,align='c',caption = 'Annual kWh/m^2')

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) %>%
  ggplot(aes(x=`W/m^2`,y=month,fill=type))+
  geom_density_ridges(alpha=.25,color='black')+
  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')