1
|
#!/usr/bin/env python
|
2
|
|
3
|
############################################################################
|
4
|
# MODULE: r.multiscalesmooth for GRASS
|
5
|
# AUTHOR(S): Original algorithm and AML implementation by John
|
6
|
# Gallant. Translation to GRASS/Python by Jim Regetz.
|
7
|
# REFERENCES:
|
8
|
# Gallant, John. 2001. Adaptive smoothing for noisy DEMs.
|
9
|
# Geomorphometry 2011.
|
10
|
# http://geomorphometry.org/system/files/Gallant2011geomorphometry.pdf
|
11
|
#
|
12
|
#############################################################################
|
13
|
|
14
|
#%module
|
15
|
#% description: Multiscale raster smoother
|
16
|
#%end
|
17
|
#%option
|
18
|
#% key: input
|
19
|
#% type: string
|
20
|
#% gisprompt: old,cell,raster
|
21
|
#% description: Raster input map
|
22
|
#% required : yes
|
23
|
#%end
|
24
|
#%option
|
25
|
#% key: output
|
26
|
#% type: string
|
27
|
#% gisprompt: new,cell,raster
|
28
|
#% description: Smoothed output raster map
|
29
|
#% required : yes
|
30
|
#%end
|
31
|
#%option
|
32
|
#% key: sd
|
33
|
#% type: string
|
34
|
#% description: Noise standard deviation; can be a raster or a constant
|
35
|
#% required : yes
|
36
|
#%end
|
37
|
#%option
|
38
|
#% key: alpha
|
39
|
#% type: double
|
40
|
#% answer: 0.05
|
41
|
#% description: Alpha used for chi-square test [0-1]
|
42
|
#% required : no
|
43
|
#%end
|
44
|
#%option
|
45
|
#% key: region
|
46
|
#% type: string
|
47
|
#% gisprompt: old,cell,raster
|
48
|
#% description: Map used to set extent; default is to use current region
|
49
|
#% required : no
|
50
|
#%end
|
51
|
|
52
|
# NOTES
|
53
|
# - make GRASS be quiet with this:
|
54
|
# gs.os.environ['GRASS_VERBOSE'] = '0'
|
55
|
|
56
|
import atexit
|
57
|
import sys
|
58
|
import math
|
59
|
import grass.script as gs
|
60
|
|
61
|
# create set to store names of temporary maps to be deleted upon exit
|
62
|
tmp_rast = set()
|
63
|
|
64
|
def cleanup():
|
65
|
gs.message("Removing temporary files...", flag='i')
|
66
|
for rast in tmp_rast:
|
67
|
gs.run_command("g.remove", rast=rast, quiet=True)
|
68
|
|
69
|
def coarsen_region(factor=3):
|
70
|
gs.run_command('g.region',
|
71
|
rows=gs.region()['rows']/3,
|
72
|
cols=gs.region()['cols']/3)
|
73
|
|
74
|
def refine_region(factor=3):
|
75
|
gs.run_command('g.region',
|
76
|
rows=gs.region()['rows']*3,
|
77
|
cols=gs.region()['cols']*3)
|
78
|
|
79
|
def multiscalesmooth(input, smooth, sd, alpha=0.05):
|
80
|
|
81
|
chisqa = (2.807 - 0.6422 * math.log10(alpha) - 3.410 *
|
82
|
math.pow(alpha, 0.3411))
|
83
|
chisqb = (-5.871 - 3.675 * math.log10(alpha) + 4.690 *
|
84
|
math.pow(alpha, 0.3377))
|
85
|
|
86
|
gs.message("Preparing initial grids", flag='i')
|
87
|
|
88
|
# create copy of initial grid of values using current region
|
89
|
gs.mapcalc('z0 = ${input}', input=input, quiet=True)
|
90
|
tmp_rast.add('z0')
|
91
|
|
92
|
# set number of aggregation levels and neighborhood size
|
93
|
NUM_LEVELS = 4
|
94
|
NUM_CELLS = 3
|
95
|
|
96
|
# expand region to accommodate integer number of cells at coarsest
|
97
|
# level, by adding roungly equal number of cells on either side
|
98
|
# (with one extra on top/right if an odd number of cells is needed)
|
99
|
max_size = NUM_CELLS**NUM_LEVELS
|
100
|
region = gs.region()
|
101
|
extra = region['cols'] % max_size
|
102
|
if (0 < extra):
|
103
|
addx = (max_size-extra)/2.0
|
104
|
else:
|
105
|
addx = 0.0
|
106
|
extra = region['rows'] % max_size
|
107
|
if (0 < extra):
|
108
|
addy = (max_size-extra)/2.0
|
109
|
else:
|
110
|
addy = 0.0
|
111
|
gs.run_command('g.region', flags='a',
|
112
|
w = region['w'] - math.floor(addx) * region['ewres'],
|
113
|
e = region['e'] + math.ceil(addx) * region['ewres'],
|
114
|
s = region['s'] - math.floor(addy) * region['nsres'],
|
115
|
n = region['n'] + math.ceil(addy) * region['nsres'])
|
116
|
gs.debug('\n'.join(gs.parse_command('g.region', 'up').keys()))
|
117
|
|
118
|
# create initial grid of variances; sd can be a raster or a constant
|
119
|
gs.mapcalc('v0 = if(isnull(z0), null(), ${sd}^2)', sd=sd, quiet=True)
|
120
|
tmp_rast.add('v0')
|
121
|
# set initial "group variance" to individual msmt variance (noise)
|
122
|
gs.run_command('g.copy', rast='v0,v.g')
|
123
|
tmp_rast.add('v.g')
|
124
|
# weights for aggregation, based on total variance
|
125
|
gs.mapcalc('w = if(isnull(v0), 0, 1/v0)', quiet=True)
|
126
|
tmp_rast.add('w')
|
127
|
# squared weights
|
128
|
gs.mapcalc('wsq = w^2', quiet=True)
|
129
|
tmp_rast.add('wsq')
|
130
|
# effective number of measurements
|
131
|
gs.mapcalc('n = if(isnull(z0), 0, 1)', quiet=True)
|
132
|
tmp_rast.add('n')
|
133
|
|
134
|
# aggregate to broader scales
|
135
|
for j in range(NUM_LEVELS):
|
136
|
i = j + 1
|
137
|
gs.message('Aggregating from %d to %d' % (j, i), flag='i')
|
138
|
gs.debug(_('Region dimensions: %d x %d' % (gs.region()['rows'],
|
139
|
gs.region()['cols'])))
|
140
|
|
141
|
# rename previous (finer scale) weights grid
|
142
|
gs.run_command('g.rename', rast='w,w.finer', overwrite=True,
|
143
|
quiet=True)
|
144
|
tmp_rast.add('w.finer')
|
145
|
|
146
|
# calc neighborhood weights, num cells, effective num cells
|
147
|
coarsen_region()
|
148
|
gs.run_command('r.resamp.stats', method='sum', input='w.finer',
|
149
|
output='w', quiet=True)
|
150
|
gs.run_command('r.resamp.stats', method='sum', input='wsq',
|
151
|
output='wsq', overwrite=True, quiet=True)
|
152
|
gs.run_command('r.resamp.stats', method='sum', input='n',
|
153
|
output='n', overwrite=True, quiet=True)
|
154
|
|
155
|
# calc variance-weighted neighborhood mean
|
156
|
refine_region()
|
157
|
gs.mapcalc('tmp = w.finer * z%d / w' % j, quiet=True)
|
158
|
tmp_rast.add('tmp')
|
159
|
coarsen_region()
|
160
|
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
161
|
output='z%d' % i, quiet=True)
|
162
|
tmp_rast.add('z%d' % i)
|
163
|
|
164
|
# calc between-cell variance, taken over neighborhood
|
165
|
refine_region()
|
166
|
gs.mapcalc('tmp = w.finer * (z%d - z%d)^2 / w' % (j, i),
|
167
|
quiet=True)
|
168
|
tmp_rast.add('tmp')
|
169
|
coarsen_region()
|
170
|
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
171
|
output='v.bg', overwrite=True, quiet=True)
|
172
|
tmp_rast.add('v.bg')
|
173
|
|
174
|
# calc wtd avg of within-cell variance, taken over neighborhood
|
175
|
if (i==1):
|
176
|
gs.mapcalc('v.wg = 0', quiet=True)
|
177
|
tmp_rast.add('v.wg')
|
178
|
else:
|
179
|
refine_region()
|
180
|
gs.mapcalc('tmp = w.finer * v.g / w', quiet=True)
|
181
|
tmp_rast.add('tmp')
|
182
|
coarsen_region()
|
183
|
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
184
|
output='v.wg', overwrite=True, quiet=True)
|
185
|
|
186
|
# calc total group variance
|
187
|
# ~= total variance of cell vals in the underlying neighborhood
|
188
|
gs.mapcalc('v.g = v.bg + v.wg', quiet=True)
|
189
|
tmp_rast.add('v.g')
|
190
|
|
191
|
# calc chisq critical values (where df = n.eff - 1)
|
192
|
gs.mapcalc('chisq = 1 + ${chisqa}/sqrt(${df}) + ${chisqb}/${df}',
|
193
|
chisqa=chisqa, chisqb=chisqb, df='(w^2/wsq - 1)', quiet=True)
|
194
|
tmp_rast.add('chisq')
|
195
|
# set coarsened cell variances: if group variance is small
|
196
|
# relative to noise variance (mean of finer scale variances),
|
197
|
# use variance of neighborhood mean instead of group variance
|
198
|
gs.mapcalc('v%d = if(v.g/chisq < ${mv}, ${vm}, v.g)' % i,
|
199
|
vm = '(1.0/w)', mv = '(n/w)', quiet=True)
|
200
|
tmp_rast.add('v%d' % i)
|
201
|
|
202
|
# get arbitrarily large value to fill null variance cells
|
203
|
bigvar = gs.raster_info('v0')['max'] * 10
|
204
|
|
205
|
# prep for refinement phase
|
206
|
gs.run_command('g.rename', rast='z%d,%s' % (NUM_LEVELS, smooth),
|
207
|
quiet=True)
|
208
|
tmp_rast.remove('z%d' % NUM_LEVELS)
|
209
|
gs.run_command('g.rename', rast='v%d,v.smooth' % NUM_LEVELS,
|
210
|
quiet=True)
|
211
|
tmp_rast.add('v.smooth')
|
212
|
tmp_rast.remove('v%d' % NUM_LEVELS)
|
213
|
|
214
|
# refine, smooth, and combine each layer in turn
|
215
|
for j in reversed(range(NUM_LEVELS)):
|
216
|
|
217
|
i = j + 1
|
218
|
gs.message('Refining from %d to %d' % (i, j), flag='i')
|
219
|
|
220
|
# create smoothed higher resolution versions of z and v
|
221
|
# using weight matrix equivalent to ArcGIS circle with radius 2
|
222
|
refine_region()
|
223
|
gs.run_command('r.neighbors', flags='ac', input=smooth,
|
224
|
output='zs', method='average', size=5, overwrite=True,
|
225
|
quiet=True)
|
226
|
tmp_rast.add('zs')
|
227
|
gs.run_command('r.neighbors', flags='ac', input='v.smooth',
|
228
|
output='vs', method='average', size=5, overwrite=True,
|
229
|
quiet=True)
|
230
|
tmp_rast.add('vs')
|
231
|
|
232
|
# combine two levels using least variance, using no-null
|
233
|
# versions of finer z and v
|
234
|
z_c = 'if(isnull(z%d), 0, z%d)' % (j, j)
|
235
|
v_c = 'if(isnull(v%d), %f, v%d)' % (j, bigvar, j)
|
236
|
gs.mapcalc('${smooth} = (${z_c}/${v_c} + zs/vs) / (1/${v_c} + 1/vs)',
|
237
|
smooth=smooth, z_c=z_c, v_c=v_c, quiet=True)
|
238
|
gs.mapcalc('v.smooth = 1 / (1/${v_c} + 1/vs)', v_c = v_c,
|
239
|
quiet=True)
|
240
|
|
241
|
cleanup()
|
242
|
return None
|
243
|
|
244
|
|
245
|
def main():
|
246
|
|
247
|
# process command options
|
248
|
input = options['input']
|
249
|
if not gs.find_file(input)['file']:
|
250
|
gs.fatal(_("Raster map <%s> not found") % input)
|
251
|
|
252
|
smooth = options['output']
|
253
|
if gs.find_file(smooth)['file'] and not gs.overwrite():
|
254
|
gs.fatal(_("Output map <%s> already exists") % smooth)
|
255
|
|
256
|
sd = options['sd']
|
257
|
try:
|
258
|
sd = float(sd)
|
259
|
except ValueError:
|
260
|
if not gs.find_file(sd)['file']:
|
261
|
gs.fatal(_("Raster map <%s> not found") % sd)
|
262
|
|
263
|
alpha = float(options['alpha'])
|
264
|
|
265
|
# set aside region for internal use
|
266
|
gs.use_temp_region()
|
267
|
|
268
|
# subset input if desired
|
269
|
region = options.get('region')
|
270
|
if region:
|
271
|
if not gs.find_file(sd)['file']:
|
272
|
gs.fatal(_("Raster map <%s> not found") % region)
|
273
|
gs.message("Setting region to %s" % region, flag='i')
|
274
|
gs.run_command('g.region', rast=region, align=input)
|
275
|
else:
|
276
|
gs.message("Using existing GRASS region", flag='i')
|
277
|
gs.debug('='*50)
|
278
|
gs.debug('\n'.join(gs.parse_command('g.region', 'p').keys()))
|
279
|
gs.debug('='*50)
|
280
|
|
281
|
multiscalesmooth(input, smooth, sd, alpha)
|
282
|
|
283
|
# restore original region
|
284
|
gs.del_temp_region()
|
285
|
|
286
|
return None
|
287
|
|
288
|
|
289
|
if __name__ == '__main__':
|
290
|
options, flags = gs.parser()
|
291
|
atexit.register(cleanup)
|
292
|
main()
|