4. Example: Using the Engerer2 Separation Model

Alright, let's do some science! Using the Engerer2 separation model.

One of the primary reasons for getting the anusolar rpackage out into the internet is that it allows anyone to use the Engerer 2 separation model, which was recently so well reviewed by Gueymard & Ruiz-Arias 2015.  As a result, I have many scientists asking for the code to compute this model, so let's do it!

Thankfully, I am super keen to share my coding environment with you.  After all, this isn't my model, its our scientific community's model. And I want you to test it, trash it, use it, improve it, and above all, do the science necessary to integrate more solar in more places around the globe. So let's check this snazzy model for separating global radiation measurements from a pyranometer into their direct and diffuse components.

Firstly, you're going to need the rpackage, so go back and follow through steps 1-3. I'll put them right here in case you need them:


Once you have the package repository cloned/forked and initialised, navigate your way to /code/examples/, where you'll find the file "MLB.diffuse_fraction.r".  Open this file in your editor of choice.

The code asks you to enter the path to the code repository in the pathto variable, do that first.

# Enter your path to repository here, end with "/"
pathto = "" # e.g. pathto = "~/anusolar/"

# load environment & initialise
source(paste(pathto,"code/rpkg/init.R",sep=""))
init(pathto)

Next up, we load the macro region. You can read more about how the data environment works in the previous post.  The macro tells the code where you are operating, in this example the macro region is already setup for Melbourne ("MLB"). 

# set the local region ('macro')
mlb = rad.meta("MLB")

Take a look at the 'mlb' metadata array, and you'll see several variables within this list object.  Included is the $loc three-letter identifier, $lat/$lon in decimal degrees, $alt for altitude - all things you can learn a bit more about in the "Working with Data" R manual page.

Take a look at the 'mlb' metadata array, and you'll see several variables within this list object.  Included is the $loc three-letter identifier, $lat/$lon in decimal degrees, $alt for altitude - all things you can learn a bit more about in the "Working with Data" R manual page.

To play around with the Engerer2 model, I've included 1-day of example radiation and weather data for you, from the 15 January 2015, so we'll load that in:

# load in example radiation data, 15 Jan 2015 in Melbourne
indf = read.rad(mlb,2013,1,15) # macro, year, month, day, resolution (sec)

# load in some weather data
wx = read.wx(mlb,2013,1,15)

If you want to load in in your own data, you'll need to feed in data in the same format. So take a look at the indf and wx data.frame objects and their variable names. I'll post more data and explanations in the future.

This example data is from the Bureau of Meteorology radiation site at the Melbourne airport, here you have 1-minute data with UTC ($gtms) and local daylight time ($ltms) timestamps.  The remaining variables are described in the in the "Working with Data" R manual page.

This example data is from the Bureau of Meteorology radiation site at the Melbourne airport, here you have 1-minute data with UTC ($gtms) and local daylight time ($ltms) timestamps.  The remaining variables are described in the in the "Working with Data" R manual page.

Next we need to create the necessary solar geometry calculations, this is accomplished through a function called "SZAQ", which stands for Solar Zenith Angle Quick (meaning that it is vectorised).

# generate solar position info for the time series
sza = SZAQ(indf,mlb)

Plotting some of the output from SZAQ, $zen is the zenith angle in decimal degrees (multiply by 100 to get degrees) and $ltms the local time.  Look at that sexy sun move across the sky, that's looking like some science.

Plotting some of the output from SZAQ, $zen is the zenith angle in decimal degrees (multiply by 100 to get degrees) and $ltms the local time.  Look at that sexy sun move across the sky, that's looking like some science.

With the solar geometry at hand, as well as the weather data we loaded in earlier (a blog post on these items in the near future), we can now use the CSMQ function.  That stands for "Clear Sky Models (Quick)", which features all of the clear-sky models in my 2015 Solar Energy validation publication.  For the purposes of the Engerer 2 model, we'll use the REST2 clear sky model, which is the preferred model for the Engerer2 computation.

# load up clear-sky model output, using REST2 model
csm = CSMQ(sza,wx,'rest2',mlb)

Here it is, warts and all (we can ignore the low sun angle behaviour, I'll fix that), a clear sky radiation estimate from the REST2 model, where the $Egc represents the global horizontal clear sky radiation in kW/m^2.

Here it is, warts and all (we can ignore the low sun angle behaviour, I'll fix that), a clear sky radiation estimate from the REST2 model, where the $Egc represents the global horizontal clear sky radiation in kW/m^2.

Next up, we can calculate the Engerer2 model to separate our global radiation measurements from the MLB radiation station.

# calculate the diffuse fraction via Engerer 2
Edh = engerer2(indf$Egh,sza,csm) # inputs from global horizontal irradiance (Egh), SZAQ, CSMQ
Edf_e = Edh/indf$Egh # calculated diffuse fraction

# note: here, Egh values should match up, row for row with sza$gtms/ltms

Lucky you, you can also compute the DIRINT model of Richard Perez...

# calculate the diffuse fraction via Perez model
Edh = perez.sep.wrapper(indf$Egh,sza,wx,mlb)
Edf_p = Edh/indf$Egh

And we can also take a look at what we just did. Let's use the classic Kt (clearness index) and DF (diffuse fraction) method of visualising the data:

# visualise:
Eext = spencer(sza)
EextH = Eext
cos(sza$zenr)
Kt = indf$Egh/EextH # clearness index

DF = indf$Edh/indf$Egh # observed diffuse fraction

plot(Kt,DF,xlim=c(0,1),ylim=c(0,1),pch=19,cex=.30) # plot obs
points(Kt,Edf_e,col='blue',pch=19,cex=0.3) # Engerer 2
points(Kt,Edf_p,col='red',pch=19,cex=0.3) # Perez*

NOTE
# radiation data should always be processed through a quality control algorithm
# the data used here are not filtered, and are raw data

Black is observed data, blue is the Engerer2 model and Red is the DIRINT model.  Looks like Richard wins this time - argh! :-)

Black is observed data, blue is the Engerer2 model and Red is the DIRINT model.  Looks like Richard wins this time - argh! :-)

And there you have it, our first functional example!

Now, for this to be truly useful, I need to cover a few more items:

  • More information about setting up your own macro
  • More information about loading in your own data
  • Providing some additional datasets
  • Examples of all the functions included in this entry (and more!)

So I've got some more work to do, but trust that I'll get around to it soon :-)

For now, I'm glad to introduce you to how to use this model!

 

/*