Cosmogenic isotope production rates over large areas by chenshu


									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.
        If you have previously delineated your basins and have already generated these
grids prior to starting this production rate calculation, take a moment to ensure all of your
grids have the same extent and cell size. MATLAB will not know how to reference
matrices with different extents or cell sizes.

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
/* extract relevant elevations
wselvs = con( ^ isnull(wsheds),elv)
/* reproject to geographic
wselvs_g = project(^ isnull(wselvs))
projection geographic
units dd
datum nad27
/* 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)
projection geographic
units dd
projection utm
zone 12
units meters
/* 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:


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,
&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.


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

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,

To top