1 |
86eb951b
|
Jim Regetz
|
/* using srtmv41 3x3 degree patch smoothed to srtmv41_smth
|
2 |
|
|
/* need to run GRID
|
3 |
|
|
|
4 |
|
|
&describe srtm_or_smth
|
5 |
|
|
&sv cellsize %grd$dx% * 111000
|
6 |
|
|
dzdx = (srtm_or_smth(1,0) - srtm_or_smth(-1,0)) / (2 * %cellsize% * cos($$ymap / deg))
|
7 |
|
|
dzdy = (srtm_or_smth(0,-1) - srtm_or_smth(0,1)) / (2 * %cellsize%)
|
8 |
|
|
slope = sqrt(sqr(dzdx) + sqr(dzdy))
|
9 |
|
|
slopedeg = atan(slope) * deg
|
10 |
|
|
/* use scaling factor that is average of x and y scaling factors, 0.707 and 1 ~= 0.85
|
11 |
|
|
/* and 0.85 / 111000 ~ 0.0000107
|
12 |
|
|
arcslope = slope(srtm_or_smth, 0.0000107, degree)
|
13 |
|
|
hillshade = hillshade(srtm_or_smth, 300, 30, shade, 0.00001)
|
14 |
|
|
|
15 |
|
|
aspect1 = atan2(0 - dzdy, 0 - dzdx) * deg
|
16 |
|
|
aspect = con(aspect1 > 90, 450 - aspect1, 90 - aspect1)
|
17 |
|
|
kill aspect1
|
18 |
|
|
|
19 |
|
|
/*use elevation
|
20 |
|
|
flowdir = flowdirection(srtm_or_smth)
|
21 |
|
|
|
22 |
|
|
/*cell area to account for latitudinal area effects, smallest area near pole
|
23 |
|
|
setwindow srtm_or_smth
|
24 |
|
|
setcell srtm_or_smth
|
25 |
|
|
cellarea = sqr(%cellsize%) * cos($$ymap / deg)
|
26 |
|
|
flowacc = flowaccumulation(flowdir, cellarea)
|
27 |
|
|
|
28 |
|
|
/*account for sinks
|
29 |
|
|
/*change # defaulto z limit
|
30 |
|
|
/*caution fill is a command not function
|
31 |
|
|
/*fill srtm_or_smth srtm_or_fill sink # flowdir_fill
|
32 |
|
|
|
33 |
|
|
flowacc_fill = flowaccumulation(flowdir_fill, cellarea)
|
34 |
|
|
|
35 |
|
|
/*scatch_fill is specific catchment area in m
|
36 |
|
|
/*wrong filling? flowdir_fill produces artefacts
|
37 |
|
|
scatch_fill = flowacc_fill / sqrt(cellarea) + sqrt(cellarea) / 2
|
38 |
|
|
/*setting the min slope to 0.1
|
39 |
|
|
adjslope = tan(max(slopedeg, 0.1) / deg)
|
40 |
|
|
|
41 |
|
|
twi_fill = ln(scatch_fill / adjslope)
|