Tutorial 3: Running TOPORAD and TOPQUAD

1. Background

The procedures for running TOPORAD and TOPQUAD will seem more clear if you take a minute to consider the components of radiative transfer through the atmosphere and the interaction of radiation with the surface.

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:

  1. Digital Elevation Model
  2. Slope and Aspect Map (from DEM)
  3. Horizon in any direction (from DEM)
  4. Sky View Factor (from DEM)
  5. Terrain Configuration Factor (from DEM)
  6. Solar zenith and azimuth
  7. Cosine of Solar Illumination Angle
  8. Atmospheric Parameters
  9. Surface Albedo

Atmospheric Parameters

There are several ways to figure out what's going on in the atmosphere that might affect solar radiation. One way is to measure the atmosphere directly, such as is routinely done by NOAA. However, often those measurements can be far away from the site we're interested in, so we look for other ways to get the information needed.

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?

2. Preparing to Run TOPQUAD

To begin, you have to generate some of the DEM-based products we talked about earlier. To do so, you'll need to learn some new IPW commands.

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]


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)


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]


b: bands to extract


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

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 (!

lq). Can you figure out why the two images look similar, when it seems they should look opposite?

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:

cut -c 1 elevs.asc > temp.asc

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:

vi temp.asc
:w! Mask.asc

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:

  1. These data will only require 1 bit, and 1 byte.
  2. For consistency below, name your IPW albedo image Mask.ipw
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

Once you've done that, come back to here and you're ready to proceed with TOPORAD!

To return to the IPW data importing tutorial, Click Here!

3. Running TOPQUAD

OK - let's do it already!

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


image: input image file

What we want to do here is put all the images that we've created together into one multi-band image:

mux DEM.ipw slas.ipw skyterr.ipw Albedo.ipw > TquadIN.ipw

It is important that you mux the files in this order!!

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]


-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)

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.

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!

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.

4. Running TOPORAD

Now let's take a look at the usage for TOPORAD:

toporad -- topographic distribution of solar radiation

Usage: toporad [-n] [image] > [newimage]


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

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:

sunang -- calculates sun angles

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]


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

So let's try sunang for Niwot, at say 10:00 a.m.

sunang -t 1996,5,1,10 -z -420 -b 40,2,30 -l 105,35,0

You should see:
GMT Wed May 1 03:00:00 1996
-z 35.289 -u 0.816251 -a 53.339

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]


w: wavelength range in um (otherwise from stdin)
d: date (year, month, day)
a: abbreviated output (no annotation, just numbers)

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.

solar -w 0.3,3.0 -d 1996,5,1

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]


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

Let's run shade:

shade -z 35.289 -a 53.339 slas.ipw > shade.ipw

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]


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

Look familiar? It is very similar to the input switches used for TOPQUAD. Let's give it a try:

elevrad -z 3600 -u 0.816251 -t 1.99 -w .71 -g .69 -r .60 -s 1267.09 DEM.ipw > erad.ipw

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:

mux erad.ipw shade.ipw skyterr.ipw Albedo.ipw > TradIN.ipw

Then go ahead and run TOPORAD (finally, eh?):

toporad TradIN.ipw > TradOUT.ipw

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:

interp | mklut -i 12 | lutx -i TradOUT.ipw > Trad8bit.ipw
0 0
4095 255

Don't worry about the error message.

That's it!

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!