Revision f3799e26
Added by Jim Regetz almost 12 years ago
| terrain/research/oregon/multiscalesmooth.py | ||
|---|---|---|
| 154 | 154 |
|
| 155 | 155 |
# calc variance-weighted neighborhood mean |
| 156 | 156 |
refine_region() |
| 157 |
gs.mapcalc('tmp = w.finer * z%d / w' % j, quiet=True)
|
|
| 157 |
gs.mapcalc('tmp = w.finer * z%d / w' % j, quiet=True,
|
|
| 158 |
overwrite=True) |
|
| 158 | 159 |
tmp_rast.add('tmp')
|
| 159 | 160 |
coarsen_region() |
| 160 | 161 |
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
| ... | ... | |
| 164 | 165 |
# calc between-cell variance, taken over neighborhood |
| 165 | 166 |
refine_region() |
| 166 | 167 |
gs.mapcalc('tmp = w.finer * (z%d - z%d)^2 / w' % (j, i),
|
| 167 |
quiet=True) |
|
| 168 |
quiet=True, overwrite=True)
|
|
| 168 | 169 |
tmp_rast.add('tmp')
|
| 169 | 170 |
coarsen_region() |
| 170 | 171 |
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
| ... | ... | |
| 173 | 174 |
|
| 174 | 175 |
# calc wtd avg of within-cell variance, taken over neighborhood |
| 175 | 176 |
if (i==1): |
| 176 |
gs.mapcalc('v.wg = 0', quiet=True)
|
|
| 177 |
gs.mapcalc('v.wg = 0', quiet=True, overwrite=True)
|
|
| 177 | 178 |
tmp_rast.add('v.wg')
|
| 178 | 179 |
else: |
| 179 | 180 |
refine_region() |
| 180 |
gs.mapcalc('tmp = w.finer * v.g / w', quiet=True)
|
|
| 181 |
gs.mapcalc('tmp = w.finer * v.g / w', quiet=True,
|
|
| 182 |
overwrite=True) |
|
| 181 | 183 |
tmp_rast.add('tmp')
|
| 182 | 184 |
coarsen_region() |
| 183 | 185 |
gs.run_command('r.resamp.stats', method='sum', input='tmp',
|
| ... | ... | |
| 185 | 187 |
|
| 186 | 188 |
# calc total group variance |
| 187 | 189 |
# ~= total variance of cell vals in the underlying neighborhood |
| 188 |
gs.mapcalc('v.g = v.bg + v.wg', quiet=True)
|
|
| 190 |
gs.mapcalc('v.g = v.bg + v.wg', quiet=True, overwrite=True)
|
|
| 189 | 191 |
tmp_rast.add('v.g')
|
| 190 | 192 |
|
| 191 | 193 |
# calc chisq critical values (where df = n.eff - 1) |
| ... | ... | |
| 234 | 236 |
z_c = 'if(isnull(z%d), 0, z%d)' % (j, j) |
| 235 | 237 |
v_c = 'if(isnull(v%d), %f, v%d)' % (j, bigvar, j) |
| 236 | 238 |
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) |
|
| 239 |
smooth=smooth, z_c=z_c, v_c=v_c, quiet=True, overwrite=True)
|
|
| 238 | 240 |
gs.mapcalc('v.smooth = 1 / (1/${v_c} + 1/vs)', v_c = v_c,
|
| 239 |
quiet=True) |
|
| 241 |
quiet=True, overwrite=True)
|
|
| 240 | 242 |
|
| 241 | 243 |
cleanup() |
| 242 | 244 |
return None |
| ... | ... | |
| 268 | 270 |
# subset input if desired |
| 269 | 271 |
region = options.get('region')
|
| 270 | 272 |
if region: |
| 271 |
if not gs.find_file(sd)['file']:
|
|
| 273 |
if not gs.find_file(region)['file']:
|
|
| 272 | 274 |
gs.fatal(_("Raster map <%s> not found") % region)
|
| 273 | 275 |
gs.message("Setting region to %s" % region, flag='i')
|
| 274 | 276 |
gs.run_command('g.region', rast=region, align=input)
|
Also available in: Unified diff
enabled overwrite and fixed input check typo