[GRASS-SVN] r54044 - grass/trunk/scripts/r.fillnulls
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 26 01:51:00 PST 2012
Author: marisn
Date: 2012-11-26 01:50:59 -0800 (Mon, 26 Nov 2012)
New Revision: 54044
Modified:
grass/trunk/scripts/r.fillnulls/r.fillnulls.html
grass/trunk/scripts/r.fillnulls/r.fillnulls.py
Log:
Fill each hole separately when run with RST. Provides significant quality improvement and up to 2x speed up.
Modified: grass/trunk/scripts/r.fillnulls/r.fillnulls.html
===================================================================
--- grass/trunk/scripts/r.fillnulls/r.fillnulls.html 2012-11-26 09:41:47 UTC (rev 54043)
+++ grass/trunk/scripts/r.fillnulls/r.fillnulls.html 2012-11-26 09:50:59 UTC (rev 54044)
@@ -10,7 +10,8 @@
Each area boundary buffer is set to three times the map resolution to get nominally
three points around the edge. This way the algorithm interpolates into the hole with
a trained slope and curvature at the edges, in order to avoid that such a flat plane
-is generated in a hole.
+is generated in a hole. The widht of edge area can be adjusted by
+changing the edge parameter.
<p>During the interpolation following warning may occur when using the RST method:<p>
<tt>
Warning: strip exists with insufficient data<br>
@@ -21,6 +22,9 @@
As the idea of <em>r.fillnulls</em> is to fill such holes, the user may
ignore the warning. The interpolation will be continued. However, the user
may pay attention to below notes.
+<p>
+If interpolation fails, temporary raster and vector maps are left in place to allow
+unfilled map hole (NULL area) identification and manual repair.
<h2>NOTES</h2>
@@ -44,6 +48,10 @@
input map and applying the "differences" color table with <em>r.colors</em>)
and/or <em>d.what.rast</em> to query individual cell values.
+<h2>WARNING</h2>
+RST method stores temporary maps on hard disk. It will require at least as much
+free space as one extra input raster map takes.
+
<h2>EXAMPLE</h2>
In this example, the SRTM elevation map in the
Modified: grass/trunk/scripts/r.fillnulls/r.fillnulls.py
===================================================================
--- grass/trunk/scripts/r.fillnulls/r.fillnulls.py 2012-11-26 09:41:47 UTC (rev 54043)
+++ grass/trunk/scripts/r.fillnulls/r.fillnulls.py 2012-11-26 09:50:59 UTC (rev 54044)
@@ -2,21 +2,22 @@
#
############################################################################
#
-# MODULE: r.fillnulls
-# AUTHOR(S): Markus Neteler
+# MODULE: r.fillnulls
+# AUTHOR(S): Markus Neteler
# Updated to GRASS 5.7 by Michael Barton
# Updated to GRASS 6.0 by Markus Neteler
# Ring and zoom improvements by Hamish Bowman
# Converted to Python by Glynn Clements
# Add support to v.surf.bspline by Luca Delucchi
-# PURPOSE: fills NULL (no data areas) in raster maps
+# Per hole filling with RST by Maris Nartiss
+# PURPOSE: fills NULL (no data areas) in raster maps
# The script respects a user mask (MASK) if present.
#
-# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
+# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
#
#############################################################################
@@ -37,6 +38,7 @@
#% description: Spline tension parameter
#% required : no
#% answer : 40.
+#% guisection: RST options
#%end
#%option
#% key: smooth
@@ -44,8 +46,18 @@
#% description: Spline smoothing parameter
#% required : no
#% answer : 0.1
+#% guisection: RST options
#%end
#%option
+#% key: edge
+#% type: integer
+#% description: Width of hole edge used for interpolation (in cells)
+#% required : no
+#% answer : 3
+#% options : 2-100
+#% guisection: RST options
+#%end
+#%option
#% key: method
#% type: string
#% description: Interpolation method
@@ -57,168 +69,290 @@
import sys
import os
import atexit
+
import grass.script as grass
-vecttmp = None
-tmp1 = None
+tmp_rmaps = list()
+tmp_vmaps = list()
usermask = None
mapset = None
# what to do in case of user break:
def cleanup():
#delete internal mask and any TMP files:
- if tmp1:
- rasts = [tmp1 + ext for ext in ['', '.buf', '_filled']]
- grass.run_command('g.remove', quiet = True, flags = 'f', rast = rasts)
- if vecttmp:
- grass.run_command('g.remove', quiet = True, flags = 'f', vect = vecttmp)
- grass.run_command('g.remove', quiet = True, rast = 'MASK')
+ if len(tmp_vmaps) > 0:
+ grass.run_command('g.remove', quiet = True, flags = 'f', vect = tmp_vmaps)
+ if len(tmp_rmaps) > 0:
+ grass.run_command('g.remove', quiet = True, flags = 'f', rast = tmp_rmaps)
if usermask and mapset:
- if grass.find_file(usermask, mapset = mapset)['file']:
- grass.run_command('g.rename', quiet = True, rast = (usermask, 'MASK'))
+ if grass.find_file(usermask, mapset = mapset)['file']:
+ grass.run_command('g.rename', quiet = True, rast = (usermask, 'MASK'), overwrite = True)
def main():
- global vecttmp, tmp1, usermask, mapset
+ global usermask, mapset, tmp_rmaps, tmp_vmaps
input = options['input']
output = options['output']
tension = options['tension']
smooth = options['smooth']
method = options['method']
-
+ edge = int(options['edge'])
+ quiet = True # FIXME
+
mapset = grass.gisenv()['MAPSET']
- unique = str(os.getpid())
+ unique = str(os.getpid()) # Shouldn't we use temp name?
+ prefix = 'r_fillnulls_%s_' % unique
+ failed_list = list() # a list of failed holes. Caused by issues with v.surf.rst. Connected with #1813
#check if input file exists
if not grass.find_file(input)['file']:
- grass.fatal(_("<%s> does not exist.") % input)
+ grass.fatal(_("Map <%s> does not exist.") % input)
- # check if a MASK is already present:
- usermask = "usermask_mask." + unique
- if grass.find_file('MASK', mapset = mapset)['file']:
- grass.message(_("A user raster mask (MASK) is present. Saving it..."))
- grass.run_command('g.rename', quiet = True, rast = ('MASK',usermask))
-
- #make a mask of NULL cells
- tmp1 = "r_fillnulls_" + unique
-
# save original region
reg_org = grass.region()
+
+ # check if a MASK is already present
+ # and remove it to not interfere with NULL lookup part
+ # as we don't fill MASKed parts!
+ if grass.find_file('MASK', mapset = mapset)['file']:
+ usermask = "usermask_mask." + unique
+ grass.message(_("A user raster mask (MASK) is present. Saving it..."))
+ grass.run_command('g.rename', quiet = quiet, rast = ('MASK',usermask))
#check if method is rst to use v.surf.rst
if method == 'rst':
- # idea: filter all NULLS and grow that area(s) by 3 pixel, then
- # interpolate from these surrounding 3 pixel edge
-
- grass.message(_("Locating and isolating NULL areas..."))
- #creating 0/1 map:
- grass.mapcalc("$tmp1 = if(isnull($input),1,null())",
- tmp1 = tmp1, input = input)
-
- #generate a ring:
- # the buffer is set to three times the map resolution so you get nominally
- # three points around the edge. This way you interpolate into the hole with
- # a trained slope & curvature at the edges, otherwise you just get a flat plane.
- # With just a single row of cells around the hole you often get gaps
- # around the edges when distance > mean (.5 of the time? diagonals? worse
- # when ewres!=nsres).
- # r.buffer broken in trunk for latlon, disabled
-
- #reg = grass.region()
- #res = (float(reg['nsres']) + float(reg['ewres'])) * 3 / 2
-
- #if grass.run_command('r.buffer', input = tmp1, distances = res, out = tmp1 + '.buf') != 0:
-
- # much easier way: use r.grow with radius=3.01
- if grass.run_command('r.grow', input = tmp1, radius = 3.01,
- old = 1, new = 2, out = tmp1 + '.buf') != 0:
- grass.fatal(_("abandoned. Removing temporary map, restoring user mask if needed:"))
-
- grass.mapcalc("MASK = if($tmp1.buf == 2, 1, null())", tmp1 = tmp1)
-
- # now we only see the outlines of the NULL areas if looking at INPUT.
- # Use this outline (raster border) for interpolating the fill data:
- vecttmp = "vecttmp_fillnulls_" + unique
- grass.message(_("Creating interpolation points..."))
- ## use the -b flag to avoid topology building on big jobs?
- ## no, can't, 'g.region vect=' currently wants to see level 2
- if grass.run_command('r.to.vect', input = input, output = vecttmp,
- type = 'point', flags = 'z'):
- grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
-
- # count number of points to control segmax parameter for interpolation:
- pointsnumber = grass.vector_info_topo(map = vecttmp)['points']
-
- grass.message(_("Interpolating %d points") % pointsnumber)
-
- if pointsnumber < 2:
- grass.fatal(_("Not sufficient points to interpolate. Maybe no hole(s) to fill in the current map region?"))
-
- # remove internal MASK first -- WHY???? MN 10/2005
- grass.run_command('g.remove', quiet = True, rast = 'MASK')
-
- # print message is a usermask it was present
- if grass.find_file(usermask, mapset = mapset)['file']:
- grass.message(_("Using user mask while interpolating"))
- maskmap = usermask
- else:
- maskmap = None
-
- grass.message(_("Note: The following 'consider changing' warnings may be ignored."))
-
- # clone current region
+ # idea: filter all NULLS and grow that area(s) by 3 pixel, then
+ # interpolate from these surrounding 3 pixel edge
+ filling = prefix + 'filled'
+
grass.use_temp_region()
- grass.run_command('g.region', vect = vecttmp, align = input)
-
- # set the max number before segmentation
- segmax = 600
- if pointsnumber > segmax:
- grass.message(_("Using segmentation for interpolation..."))
- segmax = None
+ grass.run_command('g.region', align = input, quiet = quiet)
+ region = grass.region()
+ ns_res = region['nsres']
+ ew_res = region['ewres']
+
+ grass.message(_("Using RST interpolation..."))
+ grass.message(_("Locating and isolating NULL areas..."))
+
+ # creating binary (0/1) map
+ if usermask:
+ grass.message(_("Skipping masked raster parts"))
+ grass.mapcalc("$tmp1 = if(isnull($input) && !($mask == 0 || isnull($mask)),1,null())",
+ tmp1 = prefix + 'nulls', input = input, mask = usermask)
else:
- grass.message(_("Using no segmentation for interpolation as not needed..."))
- # launch v.surf.rst
- grass.message(_("Using RST interpolation..."))
- grass.run_command('v.surf.rst', input = vecttmp, elev = tmp1 + '_filled',
- zcol = 'value', tension = tension, smooth = smooth,
- maskmap = maskmap, segmax = segmax, flags = 'z')
+ grass.mapcalc("$tmp1 = if(isnull($input),1,null())",
+ tmp1 = prefix + 'nulls', input = input)
+ tmp_rmaps.append(prefix + 'nulls')
+
+ # restoring user's mask, if present
+ # to ignore MASKed original values
+ if usermask:
+ grass.message(_("Restoring user mask (MASK)..."))
+ if grass.run_command('g.rename', quiet = quiet, rast = (usermask, 'MASK')) != 0:
+ grass.warning(_("Failed to restore user MASK!"))
+ usermask = None
+
+ # grow identified holes by X pixels
+ grass.message(_("Growing NULL areas"))
+ tmp_rmaps.append(prefix + 'grown')
+ if grass.run_command('r.grow', input = prefix + 'nulls', radius = edge + 0.01,
+ old = 1, new = 1, out = prefix + 'grown', quiet = quiet) != 0:
+ grass.fatal(_("abandoned. Removing temporary map, restoring user mask if needed:"))
+
+ # assign unique IDs to each hole or hole system (holes closer than edge distance)
+ grass.message(_("Assigning IDs to NULL areas"))
+ 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
+ grass.mapcalc("$out = if(isnull($inp), null(), $clumped)",
+ out = prefix + 'holes', inp = prefix + 'nulls', clumped = prefix + 'clumped')
+ tmp_rmaps.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()
+ cat_list = list()
+ cats_file = file(cats_file_name)
+ for line in cats_file:
+ cat_list.append(line.rstrip('\n'))
+ cats_file.close()
+ os.remove(cats_file_name)
+
+ if len(cat_list) < 1:
+ grass.fatal(_("Input map has no holes. Check region settings."))
+
+ # GTC Hole is NULL area in a raster map
+ grass.message(_("Processing %d map holes") % len(cat_list))
+ first = True
+ 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)))
+
+ # cut out only CAT hole for processing
+ if grass.run_command('g.region', rast = prefix + 'holes', align = input, 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)
+
+ # 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,
+ 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:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+ tmp_rmaps.append(holename)
+
+ # If here loop is split into two, next part of loop can be run in parallel
+ # (except final result patching)
+ # Downside - on large maps such approach causes large disk usage
+
+ # grow hole border to get it's edge area
+ tmp_rmaps.append(holename + '_grown')
+ if grass.run_command('r.grow', input = holename, radius = edge + 0.01,
+ old = -1, out = holename + '_grown', quiet = quiet) != 0:
+ grass.fatal(_("abandoned. Removing temporary map, restoring user mask if needed:"))
+
+ # no idea why r.grow old=-1 doesn't replace existing values with NULL
+ grass.mapcalc("$out = if($inp == -1, null(), $dem)",
+ out = holename + '_edges', inp = holename + '_grown', dem = input)
+ tmp_rmaps.append(holename + '_edges')
+
+ # convert to points for interpolation
+ tmp_vmaps.append(holename)
+ if grass.run_command('r.to.vect', input = holename + '_edges', output = holename,
+ type = 'point', flags = 'z', quiet = quiet) != 0:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+
+ # count number of points to control segmax parameter for interpolation:
+ pointsnumber = grass.vector_info_topo(map = holename)['points']
+ grass.verbose(_("Interpolating %d points") % pointsnumber)
- grass.message(_("Note: Above warnings may be ignored."))
-
+ if pointsnumber < 2:
+ grass.verbose(_("No points to interpolate"))
+ failed_list.append(holename)
+ continue
+
+ grass.message(_("Note: The following warnings about number of points for interpolation may be ignored."))
+
+ # set the max number before segmentation
+ # npmin and segmax values were choosen by fair guess and thus are superior to all other values ;)
+ segmax = 100
+ npmin = 300
+ if pointsnumber < 600:
+ npmin = pointsnumber
+ 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',
+ zcol = 'value', tension = tension, smooth = smooth,
+ segmax = segmax, npmin = npmin, flags = 'z') != 0:
+ # GTC Hole is NULL area in a raster map
+ grass.fatal(_("Failed to fill hole %s") % cat)
+
+ # v.surf.rst sometimes fails with exit code 0
+ # related bug #1813
+ if not grass.find_file(holename + '_dem')['file']:
+ try:
+ tmp_rmaps.remove(holename)
+ tmp_rmaps.remove(holename + '_grown')
+ tmp_rmaps.remove(holename + '_edges')
+ tmp_rmaps.remove(holename + '_dem')
+ tmp_vmaps.remove(holename)
+ except:
+ pass
+ grass.warning(_("Filling has failed silently. Leaving temporary maps with prefix <%s> for debugging.") % holename)
+ failed_list.append(holename)
+ continue
+
+ # append hole result to interpolated version later used to patch into original DEM
+ if first:
+ tmp_rmaps.append(filling)
+ grass.run_command('g.region', align = input, rast = holename + '_dem', quiet = quiet)
+ grass.mapcalc("$out = if(isnull($inp), null(), $dem)",
+ out = filling, inp = holename, dem = holename + '_dem')
+ first = False
+ else:
+ tmp_rmaps.append(filling + '_tmp')
+ grass.run_command('g.region', align = input, rast = (filling, holename + '_dem'), quiet = quiet)
+ grass.mapcalc("$out = if(isnull($inp), if(isnull($fill), null(), $fill), $dem)",
+ out = filling + '_tmp', inp = holename, dem = holename + '_dem', fill = filling)
+ if grass.run_command('g.rename', rast = (filling + '_tmp', filling),
+ overwrite = True, quiet = quiet) != 0:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+ tmp_rmaps.remove(filling + '_tmp') # this map has been removed. No need for later cleanup.
+
+ # remove temporary maps to not overfill disk
+ try:
+ tmp_rmaps.remove(holename)
+ tmp_rmaps.remove(holename + '_grown')
+ tmp_rmaps.remove(holename + '_edges')
+ tmp_rmaps.remove(holename + '_dem')
+ except:
+ pass
+ if grass.run_command('g.remove', quiet = quiet, flags = 'f', rast =
+ (holename, holename + '_grown', holename + '_edges', holename + '_dem')) != 0:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+ try:
+ tmp_vmaps.remove(holename)
+ except:
+ pass
+ if grass.run_command('g.remove', quiet = quiet, flags = 'f', vect = holename) != 0:
+ grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
+
#check if method is different from rst to use r.resamp.bspline
if method != 'rst':
- grass.message(_("Using %s bspline interpolation") % method)
+ grass.message(_("Using %s bspline interpolation") % method)
# clone current region
grass.use_temp_region()
grass.run_command('g.region', align = input)
- reg = grass.region()
- # launch r.resamp.bspline
- if grass.find_file(usermask, mapset = mapset)['file']:
- grass.run_command('r.resamp.bspline', input = input, mask = usermask,
- output = tmp1 + '_filled', method = method,
- se = 3 * reg['ewres'], sn = 3 * reg['nsres'],
- _lambda = 0.01, flags = 'n')
- else:
- grass.run_command('r.resamp.bspline', input = input,
- output = tmp1 + '_filled', method = method,
- se = 3 * reg['ewres'], sn = 3 * reg['nsres'],
- _lambda = 0.01, flags = 'n')
+ reg = grass.region()
+ # launch r.resamp.bspline
+ tmp_rmaps.append(prefix + 'filled')
+ if usermask:
+ grass.run_command('r.resamp.bspline', input = input, mask = usermask,
+ output = prefix + 'filled', method = method,
+ se = 3 * reg['ewres'], sn = 3 * reg['nsres'],
+ _lambda = 0.01, flags = 'n')
+ else:
+ grass.run_command('r.resamp.bspline', input = input,
+ output = prefix + 'filled', method = method,
+ se = 3 * reg['ewres'], sn = 3 * reg['nsres'],
+ _lambda = 0.01, flags = 'n')
# restoring user's mask, if present:
- if grass.find_file(usermask, mapset = mapset)['file']:
- grass.message(_("Restoring user mask (MASK)..."))
- grass.run_command('g.rename', quiet = True, rast = (usermask, 'MASK'))
+ if usermask:
+ grass.message(_("Restoring user mask (MASK)..."))
+ if grass.run_command('g.rename', quiet = quiet, rast = (usermask, 'MASK')) != 0:
+ grass.warning(_("Failed to restore user MASK!"))
+ usermask = None
# set region to original extents, align to input
grass.run_command('g.region', n = reg_org['n'], s = reg_org['s'],
- e = reg_org['e'], w = reg_org['w'], align = input)
+ e = reg_org['e'], w = reg_org['w'], align = input)
# patch orig and fill map
grass.message(_("Patching fill data into NULL areas..."))
# we can use --o here as g.parser already checks on startup
- grass.run_command('r.patch', input = (input,tmp1 + '_filled'), output = output, overwrite = True)
+ grass.run_command('r.patch', input = (input,prefix + 'filled'), output = output, overwrite = True)
# restore the real region
grass.del_temp_region()
@@ -227,6 +361,13 @@
# write cmd history:
grass.raster_history(output)
+
+ if len(failed_list) > 0:
+ grass.warning(_("Following holes where not filled. Temporary maps with are left in place to allow examination of unfilled holes"))
+ outlist = failed_list[0]
+ for hole in failed_list[1:]:
+ outlist = ', ' + outlist
+ grass.message(outlist)
grass.message(_("Done."))
More information about the grass-commit
mailing list