GIS production rate calculations

Cosmogenic isotope production rates over large areas...

Greg Balco
UW Cosmogenic Isotope Lab
September, 2001

Note: in order for this document to make any sense you must know quite a lot about cosmogenic isotopes and a fair amount about ARC/INFO and MATLAB. Also, nothing in this document is in any way guaranteed to work. This information is presented as a guide for people who are trying to achieve this themselves and not as a canned method for making production rate calculations. The scripts referenced below are likely to require user modifications to work on the data you want them to. A more user-friendly and wrapped-up version may follow at a later time.

The object of this exercise is to calculate the geographic distribution of cosmogenic-isotope production rates in some area, usually a watershed, using digital elevation model (DEM) data. It's necessary to do this mostly for the purpose of calculating basin-scale erosion rates from cosmognic-isotope measurements on stream sediments, and it's also useful to be able to calculate topographic shielding without actually being present at a particular point, in case someone forgot to make this measurement in the field.

This involves calculating the production rate of whatever isotope for a large number of pixels in a digital elevation model. Thus, for each pixel we need to know the latitude and altitude, to correct for these two factors, and the horizon, to calculate topographic shielding of the incoming cosmic rays.

The latitude/altitude shielding factor is easy to calculate using some combination of ARC/INFO and MATLAB. We achieve this as follows:

1. Get a DEM. Most people will be trying to do this using a USGS 30-meter DEM. The first step is to download the DEMs of the area of interest and assemble them into one ARC/INFO grid. Obviously it's important to make the coverage area big enough to inlcude anything that is likely to significantly shade the pixels of interest. Let's say this elevation grid is called "elv."

2. Define watersheds. Make another grid designating the pixels that need to have calculations done on them. For example, if you had a point coverage of stream-sediment sample locations (make sure the points are actually IN the stream defined on the DEM...) called "basinbottoms", issue the following commands in GRID:

fd2 = flowdirection(elv)
bgrid = pointgrid(basinbottoms,basinbottoms-id)
wsheds = watershed(fd2,bgrid)

Now you have a grid called wsheds in which pixels are numbered according to watersheds defined by the user-ids of the points in the point coverage.

3. Latitudes. Now you need a grid containing the latitudes of all the pixels. Since USGS 30-meter DEM's are in UTM coordinates, this is a fairly elaborate sequence of GRID commands:

/* extract relevant elevations

wselvs = con( ^ isnull(wsheds),elv)

/* reproject to geographic

wselvs_g = project(^ isnull(wselvs))
output
projection geographic
units dd
datum nad27
parameters

end

/* dump latitudes

&describe wselvs_g
setcell wselvs_g
setwindow wselvs_g
latgrid = (%grd$ymax% - (%grd$dy% / 2)) - (( $$rowmap - 1 ) * %grd$dy% )

/* reproject back into UTM
/* make sure to change UTM zone as needed

latgrid_utm = project(latgrid)
input
projection geographic
units dd
parameters
output
projection utm
zone 12
units meters
parameters
end

/* partial cleanup

kill wselvs_g
kill latgrid

/* clip

setwindow wselvs
setcell wselvs
lat = con( ^ isnull(wselvs),latgrid_utm)

Right. Now we have a grid called "lat" which contains latitudes for all the pixels in the watersheds we want to analyze.

4. Latitude/altitude correction. Get out of ARC/INFO and into MATLAB to do the latitude/altitude correction.

It's easiest to get ARC grids into MATLAB by exporting ASCII files:

elv.txt = gridascii(int(elv))
wsheds.txt = gridascii(wsheds)
lat.txt = gridascii(lat)

Then use the MATLAB text editor to remove the first six lines of each text file (header information included by ARC). In MATLAB, load these grids:

load elv.txt -ascii
load wsheds.txt -ascii
load lat.txt -ascii

Then use the m-file stone2000.m to do the latitude/altitude correction. You will also need the m-file stdatm.m .

To do this, issue the following in MATLAB:

t_elv = reshape(elv,size(elv,1)*size(elv,2),1);
t_lat = reshape(lat,size(lat,1)*size(lat,2),1);
tt_elv = elv(find(lat ~= -9999));
tt_lat = lat(find(lat ~= -9999));
tt_p = stdatm(tt_elv);
tt_corr = stone2000(tt_lat,tt_p);
lacorr = zeros(size(elv));
lacorr(find(lat~=-9999)) = tt_corr;

Now lacorr is a MATLAB variable containing the latitude/altitude correction factor for each pixel of interest. Save it and go back to ARC/INFO.

5. Topographic shielding. Next, we need to calculate the topographic shielding for all the pixels we are interested in. If you just care about one pixel, here is an ARC/INFO AML script to calculate the shielding factor for a single point in a DEM:

point_s.aml

If you care about all the pixels in a watershed, it is faster to do this in MATLAB. Keep in mind that this is a very time-consuming operation and may require many hours of machine time for a watershed of even moderate size.

First, you need to get a few more grids out of ARC into MATLAB, in particular grids containing the x and y coordinates of each pixel. This can be done in MATLAB if you know how your elevation grid is georeferenced but is easier to do in ARC. In GRID, issue:

&describe wsheds
setcell wsheds
setwindow wsheds
ygrid = (%grd$ymax% - (%grd$dy% / 2)) - ( $$rowmap * %grd$dy% )
xgrid = (%grd$xmin% + (%grd$dx% / 2)) + ( $$colmap * %grd$dx% )

One more grid is also required. Most pixels in the image are not on ridgelines and will not significantly contribute to topographic shielding of most other points. In the shielding calculation for each pixel we consider two groups of other pixels in the landscape: one, pixels near the pixel of interest, and two, pixels on ridgelines. We blow off all the other pixels. This greatly reduces execution time. The "ridgeline" pixels are most easily obtained by taking only those pixels which have a flow accumulation value of zero:

fa2 = flowaccumulation(fd2)
screen = con((fa2 EQ 0),1,0)

Then export them:

xgrid.txt = gridascii(int(xgrid))
ygrid.txt = gridascii(int(ygrid))
screen.txt = gridascii(screen)

Now back to MATLAB. Load the x and y coordinate grids, the elevation grid, the ridgeline-screening grid, and the grid with watershed identifiers. Here is an m-file to calculate topographic shielding for each pixel.

shielding.m

This results in one or more matrices containing shielding factors for each pixel of interest.

Here is an example:

Topographic shielding in Fisher Valley, Utah

6. Total isotope production rate. The last thing to do is assemble all this into a grid showing production rate in all pixels. Load the shielding and lat/alt correction matrices into MATLAB and do:

totalcorr = lacorr .* (1 - s_factor);
p = 5.1; % Stone et al production rate
p_all = totalcorr.*p;

Here is an example:

10Be production rate in Fisher Valley, Utah

7. Average P in all pixels in watershed.

This is left as an exercise.

Questions, contact Greg Balco, balcs@u.washington.edu