1 |
c87139a7
|
Jim Regetz
|
#!/usr/bin/env python
|
2 |
|
|
|
3 |
|
|
############################################################################
|
4 |
7a06885b
|
Jim Regetz
|
# MODULE: r.multiscalesmooth for GRASS
|
5 |
c87139a7
|
Jim Regetz
|
# 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 |
7a06885b
|
Jim Regetz
|
# set number of aggregation levels and neighborhood size
|
93 |
|
|
NUM_LEVELS = 4
|
94 |
|
|
NUM_CELLS = 3
|
95 |
|
|
|
96 |
c87139a7
|
Jim Regetz
|
# expand region to accommodate integer number of cells at coarsest
|
97 |
7a06885b
|
Jim Regetz
|
# 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 |
c87139a7
|
Jim Regetz
|
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 |
7a06885b
|
Jim Regetz
|
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 |
c87139a7
|
Jim Regetz
|
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 |
7a06885b
|
Jim Regetz
|
# get arbitrarily large value to fill null variance cells
|
203 |
c87139a7
|
Jim Regetz
|
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 |
3abda63e
|
Jim Regetz
|
# refine, smooth, and combine each layer in turn
|
215 |
c87139a7
|
Jim Regetz
|
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 |
7a06885b
|
Jim Regetz
|
# using weight matrix equivalent to ArcGIS circle with radius 2
|
222 |
6c9ca794
|
Jim Regetz
|
refine_region()
|
223 |
|
|
gs.run_command('r.neighbors', flags='ac', input=smooth,
|
224 |
c87139a7
|
Jim Regetz
|
output='zs', method='average', size=5, overwrite=True,
|
225 |
|
|
quiet=True)
|
226 |
|
|
tmp_rast.add('zs')
|
227 |
6c9ca794
|
Jim Regetz
|
gs.run_command('r.neighbors', flags='ac', input='v.smooth',
|
228 |
c87139a7
|
Jim Regetz
|
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()
|