1
|
/* WRR paper is based on vd5r.aml. Prior history in vd5r.aml
|
2
|
/*
|
3
|
/* mrvbf6a - vd5r modified to incorporate min (fuzzy and) instead of
|
4
|
/* multiplication when combining flatness and lowness. The transformation
|
5
|
/* following combination of flatness and lowness is no longer used, so
|
6
|
/* tvf and rvf layers are not computed, vf and rf are used directly.
|
7
|
/*
|
8
|
/* mrvbf6b - mrvbf6a modified to use hard thresholds when recombining
|
9
|
/* resolutions, rather than the smooth approach that resulted in
|
10
|
/* values bleeding into lower resolution classes
|
11
|
/*
|
12
|
/* -2 changed threshold for converting pctl to low from 0.4 to 0.5
|
13
|
/*
|
14
|
/* mrvbf6c - tried using standard deviation - not continued
|
15
|
/*
|
16
|
/* mrvbf6d - another attempt to remove small hilltops from inclusion in
|
17
|
/* valley bottom classes by including a fraction of the low indicator
|
18
|
/* from the previous scale
|
19
|
/*
|
20
|
/* -2 changed weighting of previous low/high from 0.3 to 0.5
|
21
|
/*
|
22
|
/* mrvbf6f - (ignoring the 6e experiment)
|
23
|
/* -1 Following results from Linda Gregory's transcription of the
|
24
|
/* method to Python in ArcGIS, removed the use of nibble to
|
25
|
/* extend the DEM beyond its original borders. Seems to be
|
26
|
/* unnecessary.
|
27
|
/* Included extension of grids by one pixel at each step
|
28
|
/*
|
29
|
/* mrvbf6g - all of 6f, plus std dev
|
30
|
/* Use sd thresholds with _3 applying to level 3, _9 to level 4 etc
|
31
|
/* -2 Reduce standard deviation threshold by a factor of 2, so first
|
32
|
/* threshold is 4 instead of 8. This follows analysis of -1 that
|
33
|
/* shows that ridge tops are still being left in that should be
|
34
|
/* eliminated.
|
35
|
/* -3 Reduce pctl radius from 6 cells to 3 at all steps. Radius is the
|
36
|
/* same for first two steps.
|
37
|
/* -7 As for -3, but change carrying forward of low indicator to use
|
38
|
/* raw pctl from previous scale, not low weight. This means that the
|
39
|
/* information only comes from one scale finer, not all finer scales.
|
40
|
/* This is to attempt to remedy the defects introduced by the -3 version
|
41
|
/* while retaining its ability to remove some problems.
|
42
|
/* -10 As for -8, but pctl radius at 6 cells from step 2 onwards.
|
43
|
/*
|
44
|
/* -a2 Added code to automatically detect geographic coords
|
45
|
/* Calculates sdthres from resolution, same as slopethres
|
46
|
/* -a3 Incorporate mssd calculations directly in the AML, using
|
47
|
/* correct methods developed in SRTM project; changed name of
|
48
|
/* sd layers to match slope naming
|
49
|
|
50
|
|
51
|
&args dem
|
52
|
/* dem: name of dem to work from
|
53
|
|
54
|
&sv keep 1
|
55
|
/* higher keep values result in more intermediate results being kept
|
56
|
|
57
|
&describe %dem%
|
58
|
&sv mindim = [min %grd$ncols% %grd$nrows%]
|
59
|
&sv maxresmult = 1
|
60
|
&sv nsteps = 2
|
61
|
&do &while [calc %mindim% / %maxresmult%] > [calc 5 * 3]
|
62
|
&sv maxresmult = [calc %maxresmult% * 3]
|
63
|
&sv nsteps = [calc %nsteps% + 1]
|
64
|
&end
|
65
|
&type Minimum DEM dimension = %mindim%, max resolution = %maxresmult%
|
66
|
&type %nsteps% steps
|
67
|
|
68
|
|
69
|
|
70
|
&if [null %prj$units%] &then
|
71
|
&return &error No horizontal units defined for DEM %dem%
|
72
|
|
73
|
|
74
|
&select %prj$units%
|
75
|
&when METERS
|
76
|
&do
|
77
|
&type DEM horizontal units is metres, good
|
78
|
&sv geographic = .false.
|
79
|
&end
|
80
|
|
81
|
&when DD
|
82
|
&do
|
83
|
&type DEM horizontal units is decimal degrees, scale
|
84
|
&sv geographic = .true.
|
85
|
&end
|
86
|
|
87
|
&otherwise
|
88
|
&return &error Unrecognised horizontal units for DEM %dem%: %prj$units%
|
89
|
|
90
|
&end
|
91
|
|
92
|
/* set up scaling factors for geographic data
|
93
|
&if %geographic% &then &do
|
94
|
&type Calculating geographic scaling factors
|
95
|
&sv xyfactor [calc 110000 * [cos [calc ( %grd$ymin% + %grd$ymax% ) / 2 / 57.2958] ] ]
|
96
|
&sv zfactor [calc 1 / %xyfactor% ]
|
97
|
&end
|
98
|
&else &do
|
99
|
&sv xyfactor 1
|
100
|
&sv zfactor 1
|
101
|
&end
|
102
|
|
103
|
&sv baseres = %grd$dx%
|
104
|
&sv resmetres = [calc %baseres% * %xyfactor%]
|
105
|
|
106
|
&type DEM resolution = %resmetres% metres
|
107
|
|
108
|
&sv minvbf = 1
|
109
|
&sv vbfres = 4 /* maximum resolution for first vbf = 1
|
110
|
&sv slopethres = 64 /* 64% = 33 degrees, about steepest for any deposit
|
111
|
&sv sdthres = 0.375 /* developed with 1.5 at 25 m
|
112
|
&sv pctlradius = [calc %vbfres% * 3]
|
113
|
&do &while %vbfres% < %resmetres%
|
114
|
&sv minvbf = [calc %minvbf% + 1]
|
115
|
&sv vbfres = [calc %vbfres% * 3]
|
116
|
&sv slopethres = [calc %slopethres% / 2]
|
117
|
&sv sdthres = [calc %sdthres% * 2]
|
118
|
&sv pctlradius = [calc %pctlradius% * 3]
|
119
|
&end
|
120
|
&type CHOICE setting pctlradius = 3 * dem resolution
|
121
|
&sv pctlradius = [calc %resmetres% * 3.02]
|
122
|
/* Modified 18/10/2006 to overcome numerical issues in pctl's cell inclusion calcs
|
123
|
|
124
|
&type TEMPORARY setting minvbf = 1
|
125
|
&sv minvbf = 1
|
126
|
|
127
|
&type VBF for base resolution = %minvbf%
|
128
|
&type Slope threshold = %slopethres%, percentile radius = %pctlradius%
|
129
|
|
130
|
setwindow %dem%
|
131
|
setcell minof
|
132
|
|
133
|
/* set weighting for pctl value from previous step
|
134
|
&sv prevpctlwt = 0.5
|
135
|
|
136
|
|
137
|
&sv l = 0
|
138
|
&sv v = %minvbf%
|
139
|
|
140
|
&sv l = 0
|
141
|
&sv l1 = 0
|
142
|
&sv res = %baseres%
|
143
|
|
144
|
&do &while %l% < %nsteps%
|
145
|
&if %l% = 0 &then &sv l = 1
|
146
|
&else &sv l = [calc %l% + 1]
|
147
|
&type === Step %l% ====
|
148
|
&if %l% > 2 &then &do
|
149
|
&sv res = [calc %res% * 3]
|
150
|
&sv resmetres = [calc %res% * %xyfactor%]
|
151
|
&end
|
152
|
&type Resolution = %res% (%resmetres% metres)
|
153
|
&type Slope threshold = %slopethres%, pctl radius = %pctlradius%
|
154
|
&type Std dev threshold = %sdthres%
|
155
|
|
156
|
&type Checkpoint 1
|
157
|
|
158
|
&if %l% > 2 &then &do
|
159
|
&type Checkpoint 1.1
|
160
|
/* create a smoothed version of the DEM from the previous resolution
|
161
|
&if [exists dem_%l%_%l1% -grid] &then kill dem_%l%_%l1%
|
162
|
&describe %dem%
|
163
|
setwindow [calc %grd$xmin% - %res%] [calc %grd$ymin% - %res%] [calc %grd$xmax% + %res%] [calc %grd$ymax% + %res%]
|
164
|
&if %l% = 3 &then &do
|
165
|
dem_%l%_%l1% = focalsum(%dem%, weight, gauss3, data) / focalsum(con(%dem% > -999, 1, 0), weight, gauss3, data)
|
166
|
&end
|
167
|
&else &do
|
168
|
dem_%l%_%l1% = focalsum(dem_%l1%_%l1%, weight, gauss3, data) / focalsum(con(dem_%l1%_%l1% > -999, 1, 0), weight, gauss3, data)
|
169
|
&end
|
170
|
|
171
|
&type Checkpoint 1.2
|
172
|
/* calculate standard deviation
|
173
|
&if [exists sd_%l% -grid] &then kill sd_%l%
|
174
|
&if %l% = 3 &then &do
|
175
|
sd_%l% = resample(blockstd(%dem%, rectangle, 3, 3), %res%, nearest)
|
176
|
&end
|
177
|
&else &do
|
178
|
&if [exists sswg_%l% -grid] &then kill sswg_%l%
|
179
|
&if [exists ssbg_%l% -grid] &then kill ssbg_%l%
|
180
|
sswg_%l% = aggregate(pow(sd_%l1%, 2), 3, sum)
|
181
|
ssbg_%l% = 9 * pow(resample(blockstd(dem_%l1%_%l1%, rectangle, 3, 3), %res%, nearest), 2)
|
182
|
sd_%l% = sqrt((sswg_%l% + ssbg_%l%) / 9)
|
183
|
&if %keep% <= 4 &then &do
|
184
|
kill ssbg_%l%
|
185
|
kill sswg_%l%
|
186
|
&end
|
187
|
&end
|
188
|
|
189
|
&type Checkpoint 1.3
|
190
|
&if %keep% <= 2 &then &do
|
191
|
/* get rid of DEM from previous step
|
192
|
&if %l1% > 2 &then kill dem_%l1%_%l1%
|
193
|
&if %l1% > 2 &then kill sd_%l1%
|
194
|
&end
|
195
|
&end
|
196
|
|
197
|
&type Checkpoint 2
|
198
|
|
199
|
&if [exists s_%l%_1 -grid] &then kill s_%l%_1
|
200
|
&if %l% = 1 &then &do
|
201
|
/* first step, compute slope from base resolution
|
202
|
&type Compute slope from DEM
|
203
|
s_1_1 = slope(%dem%, %zfactor%, percentrise)
|
204
|
&end
|
205
|
&else &if %l% = 2 &then &do
|
206
|
/* second step, slope is same as from first step
|
207
|
&if %keep% <= 4 &then rename s_1_1 s_2_1
|
208
|
&else copy s_1_1 s_2_1
|
209
|
&end
|
210
|
&else &do
|
211
|
/* compute slope from smoothed DEM
|
212
|
&if [exists s_%l%_%l1% -grid] &then kill s_%l%_%l1%
|
213
|
s_%l%_%l1% = slope(dem_%l%_%l1%, %zfactor%, percentrise)
|
214
|
&if %l% > 2 &then &do
|
215
|
/* resample back to base resolution
|
216
|
s_%l%_1 = resample(s_%l%_%l1%, %baseres%, bilinear)
|
217
|
&end
|
218
|
&if %keep% <= 5 &then kill s_%l%_%l1%
|
219
|
|
220
|
&if [exists dem_%l%_%l% -grid] &then kill dem_%l%_%l%
|
221
|
dem_%l%_%l% = resample(dem_%l%_%l1%, %res%)
|
222
|
&if %keep% <= 5 &then kill dem_%l%_%l1%
|
223
|
&end
|
224
|
&if [exists f_%l% -grid] &then kill f_%l%
|
225
|
f_%l% = 1 / (1 + pow(s_%l%_1 / %slopethres%, 4))
|
226
|
&if %keep% <= 4 &then &do
|
227
|
/* get rid of s_%l%_1, except when l = 1 and we need it for s_2_1
|
228
|
&if %l% ^= 1 &then kill s_%l%_1
|
229
|
&end
|
230
|
|
231
|
&type Checkpoint 3
|
232
|
|
233
|
&if [exists pctl_%l%_1 -grid] &then kill pctl_%l%_1
|
234
|
&if [exists cpctl_%l%_1 -grid] &then kill cpctl_%l%_1
|
235
|
&if [exists cf_%l% -grid] &then kill cf_%l%
|
236
|
&if [exists low_%l% -grid] &then kill low_%l%
|
237
|
&if [exists vf_%l% -grid] &then kill vf_%l%
|
238
|
&if [exists high_%l% -grid] &then kill high_%l%
|
239
|
&if [exists rf_%l% -grid] &then kill rf_%l%
|
240
|
&if %l% <= 2 &then &do
|
241
|
&r chunkproc pctl %dem% pctl_%l%_1 %pctlradius% %geographic%
|
242
|
&if %l% = 1 &then &do
|
243
|
low_%l% = 1 / (1 + pow(pctl_%l%_1 / 0.5, 3))
|
244
|
&end
|
245
|
&else &do
|
246
|
cpctl_%l%_1 = [calc 1 - %prevpctlwt%] * pctl_%l%_1 + %prevpctlwt% * pctl_%l1%_1
|
247
|
low_%l% = 1 / (1 + pow(cpctl_%l%_1 / 0.5, 3))
|
248
|
&end
|
249
|
vf_%l% = min(f_%l%, low_%l%)
|
250
|
&if %l% = 1 &then &do
|
251
|
high_%l% = 1 / (1 + pow((1 - pctl_%l%_1) / 0.5, 3))
|
252
|
&end
|
253
|
&else &do
|
254
|
high_%l% = 1 / (1 + pow((1 - cpctl_%l%_1) / 0.5, 3))
|
255
|
&end
|
256
|
rf_%l% = min(f_%l%, high_%l%)
|
257
|
&if %l% = 2 &then &do
|
258
|
cf_2 = f_2 * f_1
|
259
|
&if %keep% <= 3 &then kill f_1
|
260
|
&end
|
261
|
&end
|
262
|
&else &do
|
263
|
&if [exists pctl_%l%_%l% -grid] &then kill pctl_%l%_%l%
|
264
|
&r chunkproc pctl dem_%l%_%l% pctl_%l%_%l% %pctlradius% %geographic%
|
265
|
pctl_%l%_1 = resample(pctl_%l%_%l%, %baseres%, bilinear)
|
266
|
cpctl_%l%_1 = [calc 1 - %prevpctlwt%] * pctl_%l%_1 + %prevpctlwt% * pctl_%l1%_1
|
267
|
&if %keep% <= 5 &then kill pctl_%l%_%l%
|
268
|
cf_%l% = f_%l% * cf_%l1%
|
269
|
&if %keep% <= 3 &then kill cf_%l1%
|
270
|
low_%l% = 1 / (1 + pow(cpctl_%l%_1 / 0.5, 3))
|
271
|
&if [exists lowrelief_%l% -grid] &then kill lowrelief_%l%
|
272
|
lowrelief_%l% = 1 / (1 + pow(sd_%l% / %sdthres%, 4))
|
273
|
|
274
|
vf_%l% = min(cf_%l%, low_%l%, lowrelief_%l%)
|
275
|
high_%l% = 1 / (1 + pow((1 - cpctl_%l%_1) / 0.5, 3))
|
276
|
rf_%l% = min(cf_%l%, high_%l%, lowrelief_%l%)
|
277
|
&if %keep% <= 3 &then kill lowrelief_%l%
|
278
|
&end
|
279
|
|
280
|
&type Checkpoint 4
|
281
|
|
282
|
&if %l% > 1 &then &do
|
283
|
&if %keep% <= 4 &then kill pctl_%l1%_1
|
284
|
&if %keep% <= 4 &then kill cpctl_%l%_1
|
285
|
&if %keep% <= 3 &then kill f_%l%
|
286
|
&end
|
287
|
&if [exists mrvbf_%l% -grid] &then kill mrvbf_%l%
|
288
|
&if [exists mrrtf_%l% -grid] &then kill mrrtf_%l%
|
289
|
&if %l% = 2 &then &do
|
290
|
mrvbf_%l% = con(vf_2 >= 0.5, vf_2 + %v% - 1, vf_1)
|
291
|
mrrtf_%l% = con(rf_2 >= 0.5, rf_2 + %v% - 1, rf_1)
|
292
|
&if %keep% <= 3 &then kill vf_1
|
293
|
&if %keep% <= 3 &then kill rf_1
|
294
|
&end
|
295
|
&else &if %l% > 2 &then &do
|
296
|
mrvbf_%l% = con(vf_%l% >= 0.5, vf_%l% + %v% - 1, mrvbf_%l1%)
|
297
|
mrrtf_%l% = con(rf_%l% >= 0.5, rf_%l% + %v% - 1, mrrtf_%l1%)
|
298
|
&if %keep% <= 3 &then kill mrvbf_%l1%
|
299
|
&if %keep% <= 3 &then kill mrrtf_%l1%
|
300
|
&end
|
301
|
&if %l% >= 2 &then &do
|
302
|
&if %keep% <= 3 &then kill vf_%l%
|
303
|
&if %keep% <= 3 &then kill rf_%l%
|
304
|
&end
|
305
|
|
306
|
&type Checkpoint 5
|
307
|
|
308
|
&sv v = [calc %v% + 1]
|
309
|
|
310
|
&sv l1 = %l%
|
311
|
&sv slopethres = [calc %slopethres% / 2]
|
312
|
&if %l% > 2 &then &do
|
313
|
&sv sdthres = [calc %sdthres% * 2]
|
314
|
&end
|
315
|
&if %l% = 1 &then &sv pctlradius = [calc %pctlradius% * 2]
|
316
|
&else &sv pctlradius = [calc %pctlradius% * 3]
|
317
|
|
318
|
&type Checkpoint 6
|
319
|
&end
|
320
|
|
321
|
&if %keep% <= 4 &then kill cf_%l%
|
322
|
&if %keep% <= 4 &then kill pctl_%l%_1
|
323
|
&if %keep% <= 2 &then kill dem_%l%_%l%
|
324
|
|
325
|
&if [exists mrvbf6g-a -grid] &then kill mrvbf6g-a
|
326
|
&if [exists mrrtf6g-a -grid] &then kill mrrtf6g-a
|
327
|
mrvbf6g-a = con(%dem% > -999, con(mrvbf_%l% > [calc %minvbf% - 0.5], mrvbf_%l%, 0))
|
328
|
mrrtf6g-a = con(%dem% > -999, con(mrrtf_%l% > [calc %minvbf% - 0.5], mrrtf_%l%, 0))
|
329
|
&if %keep% <= 2 &then kill mrvbf_%l%
|
330
|
&if %keep% <= 2 &then kill mrrtf_%l%
|
331
|
|