1 |
86eb951b
|
Jim Regetz
|
/* 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 |
|
|
|