1
|
#!/usr/bin/env python
|
2
|
|
3
|
############################################################################
|
4
|
# MODULE: r.noisemag3sec 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: Calculates a measure of surface noisiness
|
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: Noise standard deviation output raster map
|
29
|
#% required : yes
|
30
|
#%end
|
31
|
#%option
|
32
|
#% key: region
|
33
|
#% type: string
|
34
|
#% gisprompt: old,cell,raster
|
35
|
#% description: Map used to set extent; default is to use current region
|
36
|
#% required : no
|
37
|
#%end
|
38
|
|
39
|
# NOTES
|
40
|
# - make GRASS be quiet with this:
|
41
|
# gs.os.environ['GRASS_VERBOSE'] = '0'
|
42
|
|
43
|
import atexit
|
44
|
import sys
|
45
|
import math
|
46
|
import grass.script as gs
|
47
|
import tempfile
|
48
|
|
49
|
# create set to store names of temporary maps to be deleted upon exit
|
50
|
tmp_rast = set()
|
51
|
|
52
|
def cleanup():
|
53
|
gs.message("Removing temporary files...", flag='i')
|
54
|
for rast in tmp_rast:
|
55
|
gs.run_command("g.remove", rast=rast, quiet=True)
|
56
|
|
57
|
def coarsen_region(factor=3):
|
58
|
gs.run_command('g.region',
|
59
|
rows=gs.region()['rows']/factor,
|
60
|
cols=gs.region()['cols']/factor)
|
61
|
|
62
|
def refine_region(factor=3):
|
63
|
gs.run_command('g.region',
|
64
|
rows=gs.region()['rows']*factor,
|
65
|
cols=gs.region()['cols']*factor)
|
66
|
|
67
|
def calculate_noise(input, output):
|
68
|
|
69
|
# create temporary annulus weights file
|
70
|
weightfile = tempfile.NamedTemporaryFile()
|
71
|
weightfile.write('0 0 1 0 0\n'
|
72
|
+ '0 1 1 1 0\n'
|
73
|
+ '1 1 0 1 1\n'
|
74
|
+ '0 1 1 1 0\n'
|
75
|
+ '0 0 1 0 0\n')
|
76
|
weightfile.flush()
|
77
|
|
78
|
#todo: do we want flag -a (don't align output with input)?
|
79
|
gs.run_command('r.neighbors', input=input, output='noisemean',
|
80
|
method='average', weight=weightfile.name, size=5, quiet=True)
|
81
|
tmp_rast.add('noisemean')
|
82
|
|
83
|
gs.mapcalc('noisediffmn = %s - noisemean' % input)
|
84
|
tmp_rast.add('noisediffmn')
|
85
|
|
86
|
gs.run_command('r.neighbors', input='noisediffmn',
|
87
|
output='noiseraw', method='stddev', flags='c', size=5,
|
88
|
quiet=True)
|
89
|
tmp_rast.add('noiseraw')
|
90
|
|
91
|
# [jg] 'calculate mean difference from mean'
|
92
|
#noisemndf = focalmean(noisediffmn, circle, 2)
|
93
|
gs.run_command('r.neighbors', input='noisediffmn',
|
94
|
output='noisemndf', method='average', flags='c', size=5,
|
95
|
quiet=True)
|
96
|
tmp_rast.add('noisemndf')
|
97
|
# [jg] 'soft absolute value of that'
|
98
|
#noisemndf_abs = sqrt(1 + sqr(noisemndf)) - 1
|
99
|
gs.mapcalc('noisemndf_abs = sqrt(1 + pow(noisemndf,2)) - 1')
|
100
|
tmp_rast.add('noisemndf_abs')
|
101
|
|
102
|
# [jg] 'smooth, this is the terrain signal'
|
103
|
#noisemndf_mdn = focalmedian(noisemndf_abs, circle, 3)
|
104
|
gs.run_command('r.neighbors', input='noisemndf_abs',
|
105
|
output='noisemndf_mdn', method='median', flags='c', size=7,
|
106
|
quiet=True)
|
107
|
tmp_rast.add('noisemndf_mdn')
|
108
|
|
109
|
# [jg] 'reduce noiseraw by the terrain signal'
|
110
|
#noiseraw_rdc = noiseraw / (1 + 2 * mndiffmn2bfmd)
|
111
|
# jr: looks like a typo in the 2nd input grid name above...
|
112
|
gs.mapcalc('noiseraw_rdc = noiseraw / (1 + 2 * noisemndf_mdn)')
|
113
|
tmp_rast.add('noiseraw_rdc')
|
114
|
|
115
|
#setcell minof
|
116
|
#noisemag_2_2 = focalmedian(aggregate(noiseraw_rdc, 2, median), circle, 2)
|
117
|
coarsen_region(factor=2)
|
118
|
gs.run_command('r.resamp.stats', method='median',
|
119
|
input='noiseraw_rdc', output='tmp', quiet=True)
|
120
|
tmp_rast.add('tmp')
|
121
|
#todo: check if we're doing this next one at the right resolution...
|
122
|
gs.run_command('r.neighbors', input='tmp',
|
123
|
output='noisemag_2_2', method='median', flags='c', size=5,
|
124
|
quiet=True)
|
125
|
tmp_rast.add('noisemag_2_2')
|
126
|
|
127
|
#setcell %DEMgrid%
|
128
|
#noisemag = resample(noisemag_2_2, #, bilinear)
|
129
|
refine_region(factor=2)
|
130
|
gs.run_command('r.resamp.interp', method='bilinear',
|
131
|
input='noisemag_2_2', output=output, quiet=True)
|
132
|
|
133
|
cleanup()
|
134
|
return None
|
135
|
|
136
|
|
137
|
def main():
|
138
|
|
139
|
# process command options
|
140
|
input = options['input']
|
141
|
if not gs.find_file(input)['file']:
|
142
|
gs.fatal(_("Raster map <%s> not found") % input)
|
143
|
|
144
|
output = options['output']
|
145
|
if gs.find_file(output)['file'] and not gs.overwrite():
|
146
|
gs.fatal(_("Output map <%s> already exists") % output)
|
147
|
|
148
|
# set aside region for internal use
|
149
|
gs.use_temp_region()
|
150
|
|
151
|
# subset input if desired
|
152
|
region = options.get('region')
|
153
|
if region:
|
154
|
if not gs.find_file(region)['file']:
|
155
|
gs.fatal(_("Raster map <%s> not found") % region)
|
156
|
gs.message("Setting region to %s" % region, flag='i')
|
157
|
gs.run_command('g.region', rast=region, align=input)
|
158
|
else:
|
159
|
gs.message("Using existing GRASS region", flag='i')
|
160
|
gs.debug('='*50)
|
161
|
gs.debug('\n'.join(gs.parse_command('g.region', 'p').keys()))
|
162
|
gs.debug('='*50)
|
163
|
|
164
|
calculate_noise(input, output)
|
165
|
|
166
|
# restore original region
|
167
|
gs.del_temp_region()
|
168
|
|
169
|
return None
|
170
|
|
171
|
|
172
|
if __name__ == '__main__':
|
173
|
options, flags = gs.parser()
|
174
|
atexit.register(cleanup)
|
175
|
main()
|