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