[GRASS-SVN] r57030 - grass/trunk/scripts/r.fillnulls
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jul 6 10:13:14 PDT 2013
Author: marisn
Date: 2013-07-06 10:13:14 -0700 (Sat, 06 Jul 2013)
New Revision: 57030
Modified:
grass/trunk/scripts/r.fillnulls/r.fillnulls.py
Log:
Speed up r.fillnulls by using vector based holes.
A modified version of Stefan Blumentrath patch from #1938
Modified: grass/trunk/scripts/r.fillnulls/r.fillnulls.py
===================================================================
--- grass/trunk/scripts/r.fillnulls/r.fillnulls.py 2013-07-06 13:57:47 UTC (rev 57029)
+++ grass/trunk/scripts/r.fillnulls/r.fillnulls.py 2013-07-06 17:13:14 UTC (rev 57030)
@@ -10,6 +10,7 @@
# Converted to Python by Glynn Clements
# Add support to v.surf.bspline by Luca Delucchi
# Per hole filling with RST by Maris Nartiss
+# Speedup for per hole filling with RST by Stefan Blumentrath
# PURPOSE: fills NULL (no data areas) in raster maps
# The script respects a user mask (MASK) if present.
#
@@ -184,17 +185,21 @@
tmp_rmaps.append(prefix + 'clumped')
if grass.run_command('r.clump', input = prefix + 'grown', output = prefix + 'clumped', quiet = quiet) != 0:
grass.fatal(_("abandoned. Removing temporary map, restoring user mask if needed:"))
-
- # use new IDs to identify holes
+
+ # get a list of unique hole cat's
grass.mapcalc("$out = if(isnull($inp), null(), $clumped)",
out = prefix + 'holes', inp = prefix + 'nulls', clumped = prefix + 'clumped')
tmp_rmaps.append(prefix + 'holes')
+ # use new IDs to identify holes
+ if grass.run_command('r.to.vect', flags = 'v', input = prefix + 'holes', output = prefix + 'holes',
+ type = 'area', quiet = quiet) != 0:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+ tmp_vmaps.append(prefix + 'holes')
+
# get a list of unique hole cat's
- cats_file_name = grass.tempfile(True)
- cats_file = file(cats_file_name, 'w')
- grass.run_command('r.stats', flags = 'n', input = prefix + 'holes', stdout = cats_file, quiet = quiet)
- cats_file.close()
+ cats_file_name = grass.tempfile(False)
+ grass.run_command('v.db.select', flags = 'c', map = prefix + 'holes', columns = 'cat', file = cats_file_name, quiet = quiet)
cat_list = list()
cats_file = file(cats_file_name)
for line in cats_file:
@@ -208,34 +213,33 @@
# GTC Hole is NULL area in a raster map
grass.message(_("Processing %d map holes") % len(cat_list))
first = True
+ hole_n = 1
for cat in cat_list:
holename = prefix + 'hole_' + cat
# GTC Hole is a NULL area in a raster map
- grass.message(_("Filling hole %s of %s") % (cat, len(cat_list)))
-
+ grass.message(_("Filling hole %s of %s") % (hole_n, len(cat_list)))
+ hole_n = hole_n + 1
# cut out only CAT hole for processing
- if grass.run_command('g.region', rast = prefix + 'holes', align = input, quiet = quiet) != 0:
+ if grass.run_command('v.extract', input = prefix + 'holes', output = holename + '_pol',
+ cats = cat, quiet = quiet) != 0:
grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
- grass.mapcalc("$out = if($inp == $catn, $inp, null())",
- out = prefix + 'tmp_hole', inp = prefix + 'holes', catn = cat)
+ tmp_vmaps.append(holename + '_pol')
- # zoom to specific hole to remove rest of data
- if grass.run_command('g.region', zoom = prefix + 'tmp_hole', align = input, quiet = quiet) != 0:
- grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
-
- # increase region a bit to fit edge zone
- if grass.run_command('g.region', align = input,
+ # zoom to specific hole with a buffer of two cells around the hole to remove rest of data
+ if grass.run_command('g.region', vect = holename + '_pol', align = input,
w = 'w-%d' % (edge * 2 * ew_res), e = 'e+%d' % (edge * 2 * ew_res),
n = 'n+%d' % (edge * 2 * ns_res), s = 's-%d' % (edge * 2 * ns_res),
quiet = quiet) != 0:
grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
- # copy only data around hole
- grass.mapcalc("$out = $inp", out = holename, inp = prefix + 'tmp_hole')
-
# remove temporary map to not overfill disk
- if grass.run_command('g.remove', flags = 'f', rast = prefix + 'tmp_hole', quiet = quiet) != 0:
+ if grass.run_command('g.remove', flags = 'f', vect = holename + '_pol', quiet = quiet) != 0:
grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+ tmp_vmaps.remove(holename + '_pol')
+
+ # copy only data around hole
+ grass.mapcalc("$out = if($inp == $catn, $inp, null())",
+ out = holename, inp = prefix + 'holes', catn = cat)
tmp_rmaps.append(holename)
# If here loop is split into two, next part of loop can be run in parallel
@@ -268,6 +272,11 @@
failed_list.append(holename)
continue
+ # Avoid v.surf.rst warnings
+ if pointsnumber < segmax:
+ npmin = pointsnumber + 1
+ segmax = pointsnumber
+
# launch v.surf.rst
tmp_rmaps.append(holename + '_dem')
if grass.run_command('v.surf.rst', quiet = quiet, input = holename, elev = holename + '_dem',
More information about the grass-commit
mailing list