|
|
@@ -131,8 +131,8 @@ def main():
|
|
|
segmax = int(options['segmax'])
|
|
|
npmin = int(options['npmin'])
|
|
|
lambda_ = float(options['lambda'])
|
|
|
- quiet = True # FIXME
|
|
|
-
|
|
|
+ quiet = True # FIXME
|
|
|
+
|
|
|
mapset = grass.gisenv()['MAPSET']
|
|
|
unique = str(os.getpid()) # Shouldn't we use temp name?
|
|
|
prefix = 'r_fillnulls_%s_' % unique
|
|
|
@@ -144,7 +144,7 @@ def main():
|
|
|
|
|
|
# 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!
|
|
|
@@ -158,16 +158,16 @@ def main():
|
|
|
# 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', 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"))
|
|
|
@@ -177,7 +177,7 @@ def main():
|
|
|
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:
|
|
|
@@ -187,7 +187,7 @@ def main():
|
|
|
except CalledModuleError:
|
|
|
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')
|
|
|
@@ -210,7 +210,7 @@ def main():
|
|
|
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
|
|
|
try:
|
|
|
grass.run_command('r.to.vect', flags='v',
|
|
|
@@ -219,7 +219,7 @@ def main():
|
|
|
except:
|
|
|
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(False)
|
|
|
grass.run_command('v.db.select', flags = 'c', map = prefix + 'holes', columns = 'cat', file = cats_file_name, quiet = quiet)
|
|
|
@@ -229,10 +229,10 @@ def main():
|
|
|
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
|
|
|
@@ -250,19 +250,20 @@ def main():
|
|
|
except CalledModuleError:
|
|
|
grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
|
|
|
tmp_vmaps.append(holename + '_pol')
|
|
|
-
|
|
|
+
|
|
|
# zoom to specific hole with a buffer of two cells around the hole to remove rest of data
|
|
|
try:
|
|
|
grass.run_command('g.region',
|
|
|
- vector=holename + '_pol', align=input,
|
|
|
+ vector=holename + '_pol', align=input,
|
|
|
w = 'w-%d' % (edge * 2 * ew_res),
|
|
|
- e = 'e+%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)
|
|
|
except CalledModuleError:
|
|
|
- grass.fatal(_("abandoned. Removing temporary maps, restoring user mask if needed:"))
|
|
|
-
|
|
|
+ grass.fatal(_("abandoned. Removing temporary maps, restoring "
|
|
|
+ "user mask if needed:"))
|
|
|
+
|
|
|
# remove temporary map to not overfill disk
|
|
|
try:
|
|
|
grass.run_command('g.remove', flags='fb', type='vector',
|
|
|
@@ -270,16 +271,16 @@ def main():
|
|
|
except CalledModuleError:
|
|
|
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
|
|
|
+
|
|
|
+ # 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')
|
|
|
try:
|
|
|
@@ -287,12 +288,12 @@ def main():
|
|
|
old=-1, out=holename + '_grown', quiet=quiet)
|
|
|
except CalledModuleError:
|
|
|
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)",
|
|
|
+ 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)
|
|
|
try:
|
|
|
@@ -301,7 +302,7 @@ def main():
|
|
|
type='point', flags='z', quiet=quiet)
|
|
|
except CalledModuleError:
|
|
|
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)
|
|
|
@@ -310,23 +311,23 @@ def main():
|
|
|
grass.verbose(_("No points to interpolate"))
|
|
|
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')
|
|
|
try:
|
|
|
grass.run_command('v.surf.rst', quiet=quiet,
|
|
|
input=holename, elev=holename + '_dem',
|
|
|
- tension=tension, smooth=smooth,
|
|
|
+ tension=tension, smooth=smooth,
|
|
|
segmax=segmax, npmin=npmin)
|
|
|
except CalledModuleError:
|
|
|
# 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']:
|
|
|
@@ -341,18 +342,18 @@ def main():
|
|
|
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, raster = holename + '_dem', quiet = quiet)
|
|
|
- grass.mapcalc("$out = if(isnull($inp), null(), $dem)",
|
|
|
+ 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, raster = (filling, holename + '_dem'), quiet = quiet)
|
|
|
- grass.mapcalc("$out = if(isnull($inp), if(isnull($fill), null(), $fill), $dem)",
|
|
|
+ grass.mapcalc("$out = if(isnull($inp), if(isnull($fill), null(), $fill), $dem)",
|
|
|
out = filling + '_tmp', inp = holename, dem = holename + '_dem', fill = filling)
|
|
|
try:
|
|
|
grass.run_command('g.rename',
|
|
|
@@ -361,7 +362,7 @@ def main():
|
|
|
except CalledModuleError:
|
|
|
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)
|
|
|
@@ -388,7 +389,7 @@ def main():
|
|
|
type='vector', name=holename)
|
|
|
except CalledModuleError:
|
|
|
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)
|
|
|
@@ -402,13 +403,13 @@ def main():
|
|
|
tmp_rmaps.append(prefix + 'filled')
|
|
|
if usermask:
|
|
|
grass.run_command('r.resamp.bspline', input = input, mask = usermask,
|
|
|
- output = prefix + 'filled', method = method,
|
|
|
- ew_step = 3 * reg['ewres'], ns_step = 3 * reg['nsres'],
|
|
|
+ output = prefix + 'filled', method = method,
|
|
|
+ ew_step = 3 * reg['ewres'], ns_step = 3 * reg['nsres'],
|
|
|
lambda_=lambda_, flags='n')
|
|
|
else:
|
|
|
grass.run_command('r.resamp.bspline', input = input,
|
|
|
- output = prefix + 'filled', method = method,
|
|
|
- ew_step = 3 * reg['ewres'], ns_step = 3 * reg['nsres'],
|
|
|
+ output = prefix + 'filled', method = method,
|
|
|
+ ew_step = 3 * reg['ewres'], ns_step = 3 * reg['nsres'],
|
|
|
lambda_=lambda_, flags='n')
|
|
|
|
|
|
# restoring user's mask, if present:
|
|
|
@@ -421,7 +422,7 @@ def main():
|
|
|
usermask = None
|
|
|
|
|
|
# set region to original extents, align to input
|
|
|
- grass.run_command('g.region', n = reg_org['n'], s = reg_org['s'],
|
|
|
+ grass.run_command('g.region', n = reg_org['n'], s = reg_org['s'],
|
|
|
e = reg_org['e'], w = reg_org['w'], align = input)
|
|
|
|
|
|
# patch orig and fill map
|
|
|
@@ -436,7 +437,7 @@ def main():
|
|
|
|
|
|
# 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]
|
Pietro Zambelli