Cosmogenic isotope production rates over large areas...
UW Cosmogenic Isotope Lab
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
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)
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.
bgrid = pointgrid(basinbottoms,basinbottoms-id)
wsheds = watershed(fd2,bgrid)
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))
/* dump latitudes
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)
/* partial cleanup
Right. Now we have a grid called "lat" which contains latitudes for
all the pixels in the watersheds we want to analyze.
lat = con( ^ isnull(wselvs),latgrid_utm)
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))
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:
wsheds.txt = gridascii(wsheds)
lat.txt = gridascii(lat)
load elv.txt -ascii
Then use the m-file stone2000.m
to do the latitude/altitude correction. You will also need the m-file
load wsheds.txt -ascii
load lat.txt -ascii
To do this, issue the following in MATLAB:
t_elv = reshape(elv,size(elv,1)*size(elv,2),1);
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.
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;
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:
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
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:
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
ygrid = (%grd$ymax% - (%grd$dy% / 2)) - ( $$rowmap * %grd$dy% )
xgrid = (%grd$xmin% + (%grd$dx% / 2)) + ( $$colmap * %grd$dx% )
fa2 = flowaccumulation(fd2)
Then export them:
screen = con((fa2 EQ 0),1,0)
xgrid.txt = gridascii(int(xgrid))
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.
ygrid.txt = gridascii(int(ygrid))
screen.txt = gridascii(screen)
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);
Here is an example:
p = 5.1; % Stone et al production rate
p_all = totalcorr.*p;
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,