O-Matrix User Stories

Estimating Cloud-free Daily Global Solar Radiation in Australia

The aim is to estimate Long-term Cloud-free Daily Global Radiation from records of the Global Radiation (which are usually affected by cloud). A large amount of satellite-derived observations is required from a grid of sites over Australia.

Australia solar radiation

Firstly, I assume that this Long-term, Cloud-free Daily Global Radiation, G, is represented by the envelope curve of all the "cloud-contaminated" observations plotted in relation to the Julian day. Second figure shows an example from a site in the Gippsland in the state of Victoria.

Secondly, I assume that G follows an annual cycle capable of approximation by a truncated Fourier series. As it turns out, only a few terms are necessary.

Thirdly, for each site a set of Fourier coeficents are chosen to form the envelope curve - the red curve in the second figure.

Using OMatrix, a technique of reiterative multiple linear regression was developed to determine these coefficients. That is to say: regression based upon all of a site's observations is performed to provide a first guess; then those observations less than the calculated first-guess envelope curve are excluded. This process of regression and data exclusion is reiterated until ~ 50 observations remain. The OMatrix rowflag and backslash operators make this very simple to program.

Fourthly, the Fourier coefficients can then be related to latitude and altitude. As an example, the zeroth coefficent a0 is plotted against latitude at third figure. Simple relationships of the Fourier coefficents with latitude and altitude are easily found.

Fifthly, armed with these relationships and the original truncated Fourier series, we arrive at an analytic expression or formula for the cloudfree global radiation for any location in Australia for any day of the year where we need only to know the site's latitude and altitude.

Australia solar radiation

The figure above illustrates Fourier components of 116 envelope curves plotted against latitude. Physical interpretation is that the annual mean of RS0 decreases poleward, the seasonality increases poleward peaking at the summer solstice, and a weak semi-annual cycle, peaking in spring and autumn, decreases poleward.

The following is the complete O-Matrix code for the reiterative multiple linear regression algorithm.
  # Import a  two column Data file. Column 1 is Daily Radiation Y and 
  # column 2 is the Julian day, J.
  # Assume that Y is approximated by a Fourier series on J
  # thus, Y = a0 + a1*cos(2*pi*J/365.25) + b1*sin(2*pi*J/365.25) + ...
  #    a2*cos(4*pi*J/365.25) + b2*cos(4*pi*J/365.25) etc
  # We want to get the 'a' and 'b' coefficients
  # Create X of same length as the Data for multiple linear regression  
  X = fillrows([0.,0.,0.,0.,0.],rowdim(Y))
  X.col(1) = ones(rowdim(Y)) 	  ;# need this to get a0
  X.col(2) = cos(2.*pi*(J)/365.25) ; # this gives a1
  X.col(3) = sin(2.*pi*(J)/365.25) ; # this gives b1
  X.col(4) = cos(4.*pi*(J)/365.25) ; # this gives a2
  X.col(5) = sin(4.*pi*(J)/365.25) ; # this gives b2
  # add more terms if needed; and then also add more columns to X
  rowflag = logical(X.col(1))
  zflag = rowflag
  a = X\Y ;# here a is {a0,a1,b1,a2,b2,...}
  # so X*a is the "predicted" Y by multiple linear regression
  	for cflag  = 1 to 200 begin  ;# 200 is arbitrary
  	if rowdim(Y.row(rowflag)) < 50 then break    ;# retains 50 or more "best" observations
  	rowflag = Y >= X*a    ;# discards observations <  "prediction"
  	Z = psort(Y/(X*a))
  	zflag = Z > 1.*cflag
  	rowflag = rowflag and zflag
  	a = X.row(rowflag)\ Y.row(rowflag) ;# reiterated 'a' coefficents
  # 'a' is now the required set of coefficents for the cloud-free or Maximum Envelope curve
  # you could use this for any observation such as temperature, pressure, 
  # etc to determine the Potential or Possible Maximum.
  # you could use this idea for any sort of cycle, most likely would be a 24 hour cycle  
This formula "predicts" the envelope curve at any of the grid sites, generating curves which visually fit as well as that shown in second figure. Generally we find that G is about 75% of the radiation that arrives at the top of the atmosphere.

Using OMatrix it is easy to manipulate and analyse a vast amount of data (some hundreds of thousands of rows of data representing 15 to 20 years of observations are required in order for any pattern to emerge), vector mean winds for each hour of each month are easily obtained and plotted.

wind plot seasonal wind plot

The hourly-monthly hodographs above show the intriguing daily and seasonal evolution of the mean wind. In general, but not always, winds in the southern hemisphere, turn counter-clockwise through the day. The opposite is true in the northern hemisphere. The wind speeds are in knots, with 5 knot ring-markers, and the diamonds represent the hourly wind, with each three-hourly plot labeled as 00, 03, 06 etc.

This user story was contributed by Warwick Grace, PhD,
Bureau of Meteorology, South Australia

Company |  Products |  Showcase |  Support |  Ordering
Copyright© 1994-2009 Harmonic Software Inc. - All rights reserved.