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:

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)

/* 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 wselvsRight. Now we have a grid called "lat" which contains latitudes for all the pixels in the watersheds we want to analyze.

setcell wselvs

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

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 -asciiThen use the m-file stone2000.m to do the latitude/altitude correction. You will also need the m-file stdatm.m .

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;

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 wshedsOne 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:

setcell wsheds

setwindow wsheds

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)

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

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

This is left as an exercise.

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