[GRASS-SVN] r33625 - grass/trunk/scripts/r.fillnulls
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 30 16:31:56 EDT 2008
Author: glynn
Date: 2008-09-30 16:31:56 -0400 (Tue, 30 Sep 2008)
New Revision: 33625
Added:
grass/trunk/scripts/r.fillnulls/r.fillnulls.py
Modified:
grass/trunk/scripts/r.fillnulls/r.fillnulls
Log:
Convert r.fillnulls to Python
Modified: grass/trunk/scripts/r.fillnulls/r.fillnulls
===================================================================
--- grass/trunk/scripts/r.fillnulls/r.fillnulls 2008-09-30 20:30:27 UTC (rev 33624)
+++ grass/trunk/scripts/r.fillnulls/r.fillnulls 2008-09-30 20:31:56 UTC (rev 33625)
@@ -140,7 +140,7 @@
g.message "Locating and isolating NULL areas..."
#creating 0/1 map:
-r.mapcalc $TMP1="if(isnull($GIS_OPT_INPUT),1,null())"
+r.mapcalc "$TMP1 = if(isnull($GIS_OPT_INPUT),1,null())"
#generate a ring:
# the buffer is set to three times the map resolution so you get nominally
@@ -161,7 +161,7 @@
exit 1
fi
-r.mapcalc "MASK=if($TMP1.buf==2,1,null())"
+r.mapcalc "MASK = if($TMP1.buf==2,1,null())"
# now we only see the outlines of the NULL areas if looking at INPUT.
# Use this outline (raster border) for interpolating the fill data:
Added: grass/trunk/scripts/r.fillnulls/r.fillnulls.py
===================================================================
--- grass/trunk/scripts/r.fillnulls/r.fillnulls.py (rev 0)
+++ grass/trunk/scripts/r.fillnulls/r.fillnulls.py 2008-09-30 20:31:56 UTC (rev 33625)
@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE: r.fillnulls
+# AUTHOR(S): Markus Neteler <neteler itc it>
+# Updated to GRASS 5.7 by Michael Barton
+# Updated to GRASS 6.0 by Markus Neteler
+# Ring improvements by Hamish Bowman
+# Converted to Python by Glynn Clements
+# PURPOSE: fills NULL (no data areas) in raster maps
+# The script respects a user mask (MASK) if present.
+#
+# COPYRIGHT: (C) 2001,2004-2005,2008 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.
+#
+#############################################################################
+
+
+#%Module
+#% description: Fills no-data areas in raster maps using v.surf.rst splines interpolation
+#% keywords: raster, elevation, interpolation
+#%End
+#%option
+#% key: input
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Raster map in which to fill nulls
+#% required : yes
+#%end
+#%option
+#% key: output
+#% gisprompt: new,cell,raster
+#% type: string
+#% description: Output raster map with nulls filled by interpolation from surrounding values
+#% required : yes
+#%end
+#%option
+#% key: tension
+#% type: double
+#% description: Spline tension parameter
+#% required : no
+#% answer : 40.
+#%end
+#%option
+#% key: smooth
+#% type: double
+#% description: Spline smoothing parameter
+#% required : no
+#% answer : 0.1
+#%end
+
+import sys
+import os
+import atexit
+import grass
+
+# what to do in case of user break:
+def cleanup():
+ #delete internal mask and any TMP files:
+ rasts = [tmp1 + ext for ext in ['', '.buf', '_filled']]
+ grass.run_command('g.remove', flags = 'f', rast = rasts)
+ grass.run_command('g.remove', flags = 'f', vect = vecttmp)
+ grass.run_command('g.remove', rast = 'MASK')
+ if grass.find_file(usermask, mapset = mapset)['file']:
+ grass.run_command('g.rename', rast = (usermask, 'MASK'))
+
+def main():
+ global vecttmp, tmp1, usermask, mapset
+
+ input = options['input']
+ output = options['output']
+ tension = options['tension']
+ smooth = options['smooth']
+
+ #check if input file exists
+ if not grass.find_file(input)['file']:
+ grass.fatal("<%s> does not exist." % input)
+
+ mapset = grass.gisenv()['MAPSET']
+ unique = str(os.getpid())
+
+ # 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', rast = ('MASK',usermask))
+
+ #make a mask of NULL cells
+ tmp1 = "r_fillnulls_" + unique
+
+ # 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:
+ e = "%s = if(isnull(%s),1,null())" % (tmp1, input)
+ grass.run_command('r.mapcalc', expression = e)
+
+ #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).
+ 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:
+ grass.fatal("abandoned. Removing temporary map, restoring user mask if needed:")
+
+ grass.run_command('r.mapcalc', expression = "MASK=if(%s.buf==2,1,null())" % 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...")
+ if grass.run_command('r.to.vect', input = input, output = vecttmp, feature = 'point'):
+ grass.fatal("abandoned. Removing temporary maps, restoring user mask if needed:")
+
+ # count number of points to control segmax parameter for interpolation:
+ s = grass.read_command('v.info', map = vecttmp)
+ for l in s.splitlines():
+ if "Number of points:" in l:
+ s = l.split(':')[1].strip()
+ s = s.split()[0].strip()
+ pointsnumber = int(s)
+ break
+
+ 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?")
+
+ grass.message("Note: The following warnings may be ignored.")
+
+ # remove internal MASK first -- WHY???? MN 10/2005
+ grass.run_command('g.remove', rast = 'MASK')
+
+ if grass.find_file(usermask, mapset = mapset)['file']:
+ grass.message("Using user mask while interpolating")
+ maskmap = usermask
+ else:
+ maskmap = None
+
+ segmax = 600
+ if pointsnumber > segmax:
+ grass.message("Using segmentation for interpolation...")
+ segmax = None
+ else:
+ grass.message("Using no segmentation for interpolation as not needed...")
+
+ grass.run_command('v.surf.rst', input = vecttmp, elev = tmp1 + '_filled',
+ zcol = 'value', tension = tension, smooth = smooth,
+ maskmap = maskmap, segmax = segmax)
+
+ grass.message("Note: Above warnings may be ignored.")
+
+ # 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', rast = (usermask, 'MASK'))
+
+ # 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.message("Filled raster map is: %s" % output)
+
+ # write cmd history:
+ grass.raster_history(output)
+
+ grass.message("Done.")
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ main()
Property changes on: grass/trunk/scripts/r.fillnulls/r.fillnulls.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list