|
1
|
/* 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)
|
