First, let's consider the effects of terrain on radiation:
Direct Solar Irradiance
The direct solar irradiance received on the surface is modulated by the cosine of the illumination angle and by the local horizon.The illumination angle is determined by the position of the sun, and by the slope and aspect of the terrain. The position of the sun at any given time can be computed easily, and since we have a DEM, we can compute slopes and aspects for all the grid elements of the DEM.
The local horizon influences how much of the sky a given grid element "sees". In this case, we're concerned about whether the sun is above or below the horizon. Since we have the DEM, we can calculate the position of the horizon in the direction of the sun.
Diffuse Irradiance
Diffuse irradiance is that portion of solar radiaiton that is scattered and therefore is not part of the direct solar beam. It is diffuse irradiance that permits us to see in shaded areas. Under overcast skies, there is no direct irradiance received at the surface, and all light available is diffuse.Diffuse irradiance is influenced by atmospheric constituents that scatter light. So, we have to know something about the nature of the atmosphere, and we also have to know how much of the atmosphere is "seen" by a grid element. Using the DEM once again, we can compute the "sky-view factor", which is the fraction of the overlying hemisphere not obscured for any given grid element.
Terrain Reflection
In rugged terrain, we also have to consider light reflected from surrounding terrain (remember that when we see those pretty mountain views, we're seeing reflected solar energy). The amount of light reflected depends on the albedo of the surface of the surrounding terrain. So, we have to know something about the albedo, and we have to know what surrounding terrain is "seen" by a given grid element. The latter - called the "terrain configuration factor" - is another item we can derive from our DEM.So - let's review the pieces we need to compute the spatial distribution of solar radiation:
- Digital Elevation Model
- Slope and Aspect Map (from DEM)
- Horizon in any direction (from DEM)
- Sky View Factor (from DEM)
- Terrain Configuration Factor (from DEM)
- Solar zenith and azimuth
- Cosine of Solar Illumination Angle
- Atmospheric Parameters
- Surface Albedo
There are models developed specifically to compute atmospheric conditions; two of these are called LOWTRAN and MODTRAN. Some remote sensing algorithms can also be used to determine atmospheric parameters.
There are three atmospheric parameters of interest here - tau (t), omega (w), and gamma (g). Tau is called the optical depth, omega is the single-scattering albedo (of atmospheric particles, essentially), and gamma is the scattering phase function. For now, you don't need to worry about what these mean, only about what values of them to use in TOPORAD or TOPQUAD.
One important thing we haven't mentioned yet is that the optical properties (t,w,g) are different for different wavelengths. The atmosphere may be quite transparent for one wavelength region, but almost opaque for another. Therefore, we need to consider this in determining which parameters to use in TOPORAD and TOPQUAD.
With the help of the model LOWTRAN7, here are some values for the atmospheric parameters for 27 bands (width = 0.1 um) within the visible and near-Infrared portion of the electromagnetic spectrum, for very clear skies:
wavelength tau omeg gamma
0.3-0.4 5.03 1.00 0.67 0.4-0.5 3.56 0.95 0.73 0.5-0.6 2.78 0.95 0.77 0.6-0.7 2.25 0.95 0.79 0.7-0.8 1.91 1.00 0.82 0.8-0.9 1.63 0.95 0.81 0.9-1.0 1.46 1.00 0.83 1.0-1.1 1.16 0.87 0.79 1.1-1.2 1.18 0.99 0.82 1.2-1.3 0.94 0.49 0.49 1.3-1.4 1.67 0.58 0.70 1.4-1.5 1.17 0.98 0.82 1.5-1.6 0.65 0.98 0.82 1.6-1.7 0.59 0.95 0.81 1.7-1.8 0.58 0.14 0.50 1.8-1.9 1.98 0.46 0.63 1.9-2.0 1.10 0.67 0.68 2.0-2.1 0.62 0.97 0.81 2.1-2.2 0.39 0.96 0.80 2.2-2.3 0.40 0.91 0.80 2.3-2.4 0.44 0.91 0.80 2.4-2.5 0.53 0.82 0.79 2.5-2.6 4.07 0.12 0.43 2.6-2.7 7.01 0.01 0.27 2.7-2.8 7.75 0.01 0.23 2.8-2.9 2.27 0.09 0.44 2.9-3.0 0.71 0.36 0.80 Mean 1.99 0.71 0.69
Note how the values change for different wavelengths. For example, the optical depth (t) in the UV (0.3-0.4) is very high (due largely to the presence of ozone, which absorbs UV effectively), which means less radiation in that wavelength region makes it through (which is good).
It is also important to consider Planck's Law and Wien's Displacement Law concerning the distribution of solar radiation from a body at a given temperature, such as the sun. At 6000 K, the most intense solar radiation occurs down in the visible region, with the magnitude of radiation tapering off rapidly out near 3.0 um. So the high optical depths between 2.5 and 2.8 um occur in a part of the EM spectrum that isn't receiving a large magnitude of radiation in the first place.
To model solar irradiance at the surface, do you think it would be better to
model the irradiance in several different wavelength regions and integrate the
total, or do you think it would be just as effective to use the mean parameters
and model the irradiance for the entire spectrum at once?
First, we'll generate the slope and aspect images from the DEM, using the IPW command gradient.
gradient -- slope and aspect of image
Usage: gradient [-s] [-a] [-i bits[,bits]] [-d delta[,delta]] [image] > [newimage]
Options:
s: compute slope image only (default both)
a: compute aspect image only
i: # bits in output linear quantization (default CHAR_BIT)
d: grid spacing (line & samp, normally gotten from GEOH header)
Operands:
image: input elevation image
newimage: output elevation image
gradient DEM.ipw > slas.ipw
Unless you specify either -s or -a, gradient outputs a 2-band image, the first band is slope (band 0), and the second band is aspect (band 1). The xv software doesn't know what to do with 2 bands, so if you want to look at them you have to split them apart. We'll go ahead and split them, but we'll have to put them back together later. To split them, we'll use the program demux:
demux -- extract image bands
Usage: demux -b band,... [image] > [newimage]
Options:
b: bands to extract
Operands:
image: input image file
newimage: output image file
Note that above I said the first band is band 0, the second is band 1, etc. This is the convention in IPW, and is worth remembering so that in the future you don't go out and grab the wrong band. To demux the slas.ipw file into two seperate files, you'll run demux twice:
demux -b 0 slas.ipw > slope.ipw
demux -b 1 slas.ipw > aspect.ipw
Now, take a look at each image using xv2 (Note: for convenience, you can open up two windows and run xv2 in each to look at both images at the same time).
xv2 slope.ipw
xv2 aspect.ipw
The slope image should appear to be very dark. Using the color editing function of xv, adjust the intensity curve until you can see features in the image.
Now it's time to compute the sky view and terrain configuration factors, using the single IPW command viewf:
viewf -- sky view and terrain configuration
Usage: viewf elev_image > newimage
This also will be a two band file, where band 0 is the sky-view factor and
band 1 is the terrain configuration factor. The sky view factor is defined
as the fraction (0..1) of a pixels overlying hemisphere (centered at the
pixel's zenith) that is subtended by sky (as opposed to surrounding terrain).
The terrain configuration factor is defined as
1 + cos(slope) -------------- - sky view factor 2
So, by these definitions, you would expect these two images look sort of opposite of each other - if a cell sees a lot of terrain, it doesn't see a lot of sky, and vice versa. Maybe you should check it out!
viewf DEM.ipw > skyterr.ipw
demux -b 0 skyterr.ipw > skyview.ipw
demux -b 1 skyterr.ipw > terrain.ipw
Now take a look at them as you did before, using xv2.
Do they look like they are "opposite" images?
Take a look at the headers of each image using the UNIX more command. Look for the linear quantization header (!
OK. Now you have to come up with an albedo image, and you'll be all set to run TOPORAD. If you had some remote sensing data, you might be able to come up with a map of surface albedo. For now, we'll fake it, and assume that the albedo
is constant throughout the area. The examples in the exercises following this tutorial are for snowy conditions, so let's assume that the albedo of the region
is 0.6. How can we make an albedo image with 0.6's everywhere?
One easy way is to go back to the method we used to import the ASCII elevation
values. Using the UNIX editor vi, we can change our file of 74025
elevation values to 1.0's. This is not something you want to do manually, so
here's a suggestion for how to do it. Use the UNIX cut -c command to cut out the first column:
Now you should have a file with either 3's or 4's in it (all the elevations wereeither in the 3 thousands or 4 thousands meters). Then use vi's global search
and replace function to replace all the 3's and 4's with 1's:
Now you have your ASCII file with 74025 1's, and you can follow the
instructions from the first tutorial to read them in, with the following modifications:
Once you've done that, come back to here and you're ready to proceed with TOPORAD!
I said before that we'd have to put some files back together. To do that, we'll use the IPW mux command (if you can demux, you must be able to mux!):
mux -- multiplex (band-interleave) images
Usage: mux image ... > newimage
Operands:
What we want to do here is put all the images that we've created together into
one multi-band image:
Let's take a look at the usage for TOPQUAD:
topquad -- daily integrated radiation over topographic grid
Usage: topquad [-n] -z elev -t tau -w omega -g g -r R0 [-s S0] [-x w1,w2] \
-d y,m,d -b d,m,s -l d,m,s [image]
Options:
Let's test it out on Niwot Ridge, for May 1, 1996. We'll say that z = 3600 m, and we'll use the mean atmospheric parameters noted above for very clear skies,
running it once for the broad wavelength range 0.3 - 3.0 um, even though we
know it should be run for narrower wavelength regions.
When it's finished, view the output using xv2. You should have an image that
distinctly reveals the terrain of Niwot Ridge. The output data are in units of meagajoules (MJ), but you will need to look at the LQ header again to see
what the range of units actually is (should be ~116 - 147 MJ).
This output is pretty low for clear sky conditions at this time of year. The
reason is that we used the mean atmospheric parameters, which aren't really
appropriate for what we're trying to do. If we ran TOPQUAD 27 times, using
the atmopspheric parameters for each 0.1 um spectral band, then added all 27
together, we'd come much closer to the total daily irradiance.
toporad -- topographic distribution of solar radiation
Usage: toporad [-n] [image] > [newimage]
Options:
Operands:
If you recall, TOPQUAD calls TOPORAD several times between sunrise and sunset.
Not surprisingly, then, the inputs to the two models are similar. However, we have to mux together a slightly different input file for TOPORAD. The reason is that TOPQUAD calls a few IPW programs (sunang, solar, shade, and elevrad), which return
information needed to run TOPORAD. If you run TOPORAD by itself, you have to
run these programs yourself.
First, let's look at the IPW program sunang:
Usage: sunang -t i,i,i,i[,i,i] -b lat[,lat,lat] -l lon[,lon,lon] [-s slope]
[-a azm] [-r] [-z zone] [-y] [-d]
Options:
So let's try sunang for Niwot, at say 10:00 a.m.
You should see:
where -z is the zenith angle, -u is the cosine of the zenith angle, and -a is the solar azimuth. This is a handy little program to remember, in case you
suddenly need to know the solar zenith and azimuth!
Now let's look at the IPW program solar, which calculates the exoatmospheric solar radiation:
solar -- exoatmospheric direct solar irradiance
Usage: solar [-w um[,um]] [-d d,d,d] [-a]
Options:
We want to be sure we specify the same wavelength range that we plan to use
for the model run, so for now we'll stick to 0.3-3.0 um.
You should get the answer 1267.09 W m^-2.
Now we'll look at the IPW program shade. You supply shade with the solar zenith and azimuth, and using the slope/aspect image you generated previously it will compute the cosine of the local illumination angle for all of the
cells in the image. Let's look at it's usage:
shade -- cosine of illumination angle
Usage: shade [-z zenith] [-u cos] -a azimuth [image] > [newimage]
Options:
Operands:
Let's run shade:
Now if you view your results using xv2, you should see an image that looks like a shaded relief image illuminated from the south-southeast. In fact, you can use shade to make very nice shaded relief images if you select a solar azimuth to the northwest. The larger you make the zenith angle, the greater the effect.
The last IPW program to run before running TOPORAD is elevrad, which calculates the beam and diffuse irradiance that would be incident on a flat surface.After that, TOPORAD will adjust the irradiances for all the slopes and aspects.
Let's look at elevrad:
elevrad -- beam and diffuse radiation (preprocessor for toporad)
Usage: elevrad -z elev -u mu0 -t tau -w omega -g gfact -r R0 -s S0 [-n bits[,bits]] [image] > [newimage]
Options:
Operands:
Look familiar? It is very similar to the input switches used for TOPQUAD. Let's give it a try:
You'll get a 2-band, 12-bit image; you can demux and requantize to view in xv2 if you like.
Now, you need to mux together the multiband input file for TOPORAD, again in the specified order:
Then go ahead and run TOPORAD (finally, eh?):
The output file will be a single band, 12-bit file with instantaneous irradiance in W m^-2 for 10:00 a.m. on May 1, 1996, under "clear" skies. To view the result, you need to do the exact same procedure that you followed previously to view your DEM:
That's it!
cut -c 1 elevs.asc > temp.asc
vi temp.asc
:1,74025s/3/1/g
:1,74025s/4/1/g
:w! Mask.asc
:q!
Now you have an ipw image with every pixel equal to 1, but we want an albedo
image with every cell equal to 0.6. The ascii to binary conversion cannot
convert a floating point number (decimal), so we had to read in ones instead.
We can use a linear quantization header on the ipw image to tell ipw that the
ones are to be interpreted as 0.6's:
mklqh -m 0,0,1,.6 Mask.ipw > Albedo.ipw
To return to the IPW data importing tutorial,
Click Here!
3. Running TOPQUAD
OK - let's do it already!
image: input image file
mux DEM.ipw slas.ipw skyterr.ipw Albedo.ipw > TquadIN.ipw
It is important that you mux the files in this order!!
-n: net radiation instead of incoming
-z: elevation of optical depth measurement
-t: optical depth at z
-w: single scattering albedo
-r: reflectance of substrate
-s: exoatmospheric solar irradiance
-x: wavelength range, micrometers
-d: date (year, month, day)
-b: latitude (degrees, minutes, seconds)
-l: longitude (degrees, minutes, seconds)
topquad -z 3600 -t 1.99 -w .71 -g .69 -r .60 -x .3,3.0 -d 1996,5,1 -b 40,2,30 -l 105,35,0 TquadIN.ipw > TquadOUT.ipw
Type it on 1 line; it will take a minute or two to run, so be patient!
4. Running TOPORAD
Now let's take a look at the usage for TOPORAD:
n: net radiation instead of incoming
image: beam & diffuse irrad, mu_s, Vf, Ct, albedo
newimage: output image of instantaneous irradiance, in W/m^2
sunang -- calculates sun angles
t: date & time: year, month, day, hour [, min, sec]
b: latitude: deg, min, sec unless -r specified
l: longitude: deg, min, sec unless -r specified
s: slope of surface: deg unless -r specified
a: azimuth of surface: deg unless -r specified
r: lat/lon in radians
z: time zone, minutes west of Greenwich (GMT if omitted)
y: daylight savings flag
d: print Earth-Sun radius vector also
sunang -t 1996,5,1,10 -z -420 -b 40,2,30 -l 105,35,0
GMT Wed May 1 03:00:00 1996
-z 35.289 -u 0.816251 -a 53.339
w: wavelength range in um (otherwise from stdin)
d: date (year, month, day)
a: abbreviated output (no annotation, just numbers)
solar -w 0.3,3.0 -d 1996,5,1
z: solar zenith angle, degrees or
u: cosine solar zenith angle
a: solar azimuth, degrees_ccw from south (+ E, - W)
image: input slope/azimuth file
newimage: output file
shade -z 35.289 -a 53.339 slas.ipw > shade.ipw
z: elevation of optical depth measurement
u: cosine solar zenith angle
t: optical depth at z
w: single-scattering albedo
g: scattering asymmetry parameter
r: mean surface albedo
s: exoatmospheric solar irradiance
n: # output bits (default same as input)
image: input elevation image
newimage: output image
elevrad -z 3600 -u 0.816251 -t 1.99 -w .71 -g .69 -r .60 -s 1267.09 DEM.ipw > erad.ipw
mux erad.ipw shade.ipw skyterr.ipw Albedo.ipw > TradIN.ipw
toporad TradIN.ipw > TradOUT.ipw
interp | mklut -i 12 | lutx -i TradOUT.ipw > Trad8bit.ipw
0 0
4095 255
Return
Return
Don't worry about the error message.
If you want to return to the Overview Page,
Click Here!
If you want to return to the IPW data importing tutorial,
Click Here!
If you want to return to the Viewing IPW images tutorial,
Click Here!
If you want to go on to the Exercises tutorial,
Click Here!