[GRASS-SVN] r38144 - grass-addons/vector/v.autokrige
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jul 1 10:45:33 EDT 2009
Author: mathieug
Date: 2009-07-01 10:45:33 -0400 (Wed, 01 Jul 2009)
New Revision: 38144
Modified:
grass-addons/vector/v.autokrige/v_autokrige.py
Log:
Output handling improvements
Modified: grass-addons/vector/v.autokrige/v_autokrige.py
===================================================================
--- grass-addons/vector/v.autokrige/v_autokrige.py 2009-07-01 14:39:07 UTC (rev 38143)
+++ grass-addons/vector/v.autokrige/v_autokrige.py 2009-07-01 14:45:33 UTC (rev 38144)
@@ -22,27 +22,27 @@
#############################################################################
#%Module
-#% description: automatic kriging interpolation from vector point data
-#% keywords: kriging, autokrige, RGrass
+#% description: automatic ordinary kriging interpolation from vector point data with R automap package.
+#% keywords: kriging, automap, R
#%End
#%option
#% key: input
#% type: string
#% gisprompt: old,vector,vector
-#% description: Name of the vector containing sites data (values to interpolate) and geometry. Must not contain any null values.
+#% description: Name of the vector containing values to interpolate (must not contain any null values)
#% required : yes
#%end
#%option
#% key: column
#% type: string
-#% description: Attribute column to interpolate (must be numeric).
+#% description: Attribute column to interpolate (must be numeric)
#% required : yes
#%end
#%option
#% key: nbcell
#% type: integer
#% answer: 100
-#% description: number of cell of the created grid. If set, GRASS region is affected.
+#% description: number of cell of the created grid. If set, GRASS region is affected
#% required : no
#%end
#%option
@@ -51,13 +51,13 @@
#% options: Exp,Sph,Gau,Mat
#% multiple : yes
#% answer:
-#% description: List of variogram models to be automatically tested. Select only one to fix the model.
+#% description: List of variogram models to be automatically tested. Select only one to fix the model
#% required : no
#%end
#%option
#% key: range
#% type: integer
-#% description: Range value. Automatically fixed if not set.
+#% description: Range value. Automatically fixed if not set
#% required : no
#%end
#%option
@@ -75,13 +75,13 @@
#%option
#% key: nmax
#% type: double
-#% description: Local kriging. Number of nearest observations that should be used for a kriging prediction. Default : all observations are used.
+#% description: Local kriging. Number of nearest observations that should be used for a kriging prediction (default : all observations are used)
#% required : no
#%end
#%option
#% key: maxdist
#% type: double
-#% description: Local kriging. Maximum distance for observations from the prediction location to be included for prediction. Default : no limitation.
+#% description: Local kriging. Maximum distance for observations from the prediction location to be included for prediction (default : no limitation)
#% required : no
#%end
#%option
@@ -104,17 +104,17 @@
#%end
#%flag
#% key: r
-#% description: Don't set region from sites extent. Region would have been previously set in grass.
+#% description: Don't set region from sites extent. Region would have been previously set in grass
#%end
#%flag
#% key: l
-#% description: Log R output to v_autokrige.py.log.
+#% description: Log R output to v_autokrige.py.log
#%end
import sys
import os
import re
-from subprocess import Popen, PIPE
+from subprocess import Popen, call, PIPE
import traceback
##see http://trac.osgeo.org/grass/browser/grass/trunk/lib/python
@@ -145,61 +145,10 @@
self.RscriptFile = None
logfilename = 'v_autokrige.py.log'
self.logfile = os.path.join(os.getenv('LOGDIR'),logfilename) if os.getenv('LOGDIR') else logfilename
-
- def __checkLayers(self, input, output):
- """
- Preliminary checks before starting kriging.
- Note : for this to work with grass6.3, in find_file function from core.py,
- command should be (n flag removed because 6.4 specific):
- s = read_command("g.findfile", element = element, file = name, mapset = mapset)
- """
- ##Test for input vector map.
- testInput = grass.find_file(input, element = 'vector')
- if testInput['fullname'] == '':
- raise AutoKrigeError("vector map " + input + " not found in mapset search path.")
- ##Test if output raster map already exists.
- testOutput = grass.find_file(output, element = 'cell')
- if testOutput['fullname'] != '':
- if not grass.overwrite() is True:
- raise AutoKrigeError("raster map " + output + " already exists in mapset search path. \n \
-Use the --o flag to overwrite.")
- else:
- print "Warning: raster map " + output + " will be overwritten."
- ##Test also variance raster.
- if self.varianceFlag is True:
- testVarOutput = grass.find_file(output + '_var', element = 'cell')
- if testVarOutput['fullname'] != '':
- if not grass.overwrite() is True:
- raise AutoKrigeError("raster map " + output + " already exists in mapset search path. \n \
-Use the --o flag to overwrite.")
- else:
- print "Warning: raster map " + output + '_var' + " will be overwritten."
-
-
- def __getGridCellSize(self, input, nbcell):
- """Define kriged grid cell size. Raster resolution but also computation time depends on it.
- We take region resolution as cell size but we can fix this resolution with the nbcell parameter.
- Only one value is needed because the R script use square cells."""
- if self.regionFlag is not True:
- grass.run_command("g.region", vect = input)
- if nbcell is not None:
- self.__fixRegionResFromNumberOfCells(nbcell)
- regionParams = grass.region()
- cellsize = float(regionParams['nsres']) if float(regionParams['nsres']) > float(regionParams['ewres']) \
- else float(regionParams['ewres'])
- return cellsize
-
- def __fixRegionResFromNumberOfCells(self, nbcell=100):
- """Adjust one of the two region dimensions so that we have a maximum of 'nbcell' cells in both."""
- regionParams = grass.region()
- nsResForGivenNbCell = (float(regionParams['n'])- float(regionParams['s'])) / float(nbcell)
- ewResForGivenNbCell = (float(regionParams['e'])- float(regionParams['w'])) / float(nbcell)
- tmpRegionRes = nsResForGivenNbCell if nsResForGivenNbCell > ewResForGivenNbCell else ewResForGivenNbCell
- grass.use_temp_region
- grass.run_command("g.region", res = tmpRegionRes)
-
+ grass.try_remove(self.logfile)
+
def __prepareRModelsString(self, models=None):
- """Try to create the R argument as expected in the R script,
+ """!Create the R argument as expected in the R script,
to avoid string manipulations with R."""
modelsString = ''
if models is None:
@@ -215,6 +164,7 @@
return modelsString
def __writeRScript(self):
+ """!Create the R script for automap"""
RscriptFile = grass.tempfile()
fileHandle = open(RscriptFile, 'a')
script = """
@@ -303,7 +253,6 @@
attr(sitesR, "proj4string") <-CRS(G$proj4)
cat("ordinary kriging","\n")
- #[note : ajouter une option pour permettre de réaliser le krigeage universel]
kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), sitesR[column], mask_SG, model = modelslist, fix.values = c(nugget,range,sill), debug.level=-1, verbose=TRUE, nmax=nobsmax, maxdist=maxdistance)
cat("send raster to GRASS","\n")
@@ -324,12 +273,14 @@
cat("R script done","\n")
quit(status = 0)
}, interrupt = function(ex) {
- cat("An interrupt was detected while executing R script.\n");
+ cat("An interrupt was detected while executing R script.\n")
+ cat(conditionMessage(ex), "\n")
quit(status = 4)
},
error = function(ex) {
- cat("Error while executing R script.\n");
- cat(conditionMessage(ex), "\n")
+ cat("\n");
+ cat("Error while executing R script\n")
+ cat(conditionMessage(ex), "\n")
quit(status = 5)
}
) ## tryCatch()
@@ -339,52 +290,120 @@
fileHandle.close()
return RscriptFile
- def __finalize(self, output, colormap):
+ def __execRCommand(self, command, logFile = False):
+ """!R command execution method using Popen in shell mode"""
+ if logFile is not False:
+ command = command + ' >> ' + logFile
+ ##use redirection on stdout only
+ command = command + " 2>&1"
+ p = Popen(command, shell = True, stdout = PIPE)
+ retcode = p.wait()
+ com = p.communicate()
+ if retcode == 0:
+ return retcode
+ else:
+ lines = None
+ if logFile is not False:
+ logHandle = open(logFile,"r")
+ lines = logHandle.readlines()
+ else:
+ r = re.compile('\n')
+ lines = r.split(com[0])
+ out = lines[-2].strip() + ":" + lines[-1].strip()
+ raise AutoKrigeError(out)
+
+ def printMessage(self, message, type = 'info'):
+ if type == 'error':
+ grass.error(message)
+ elif type == 'warning':
+ grass.warning(message)
+ elif type == 'info' and grass.verbosity > 0:
+ grass.info(message)
+
+ def checkLayers(self, input, output, testVarianceRast = False):
"""
- Group operations non-related to kriging.
+ !Preliminary checks before starting kriging.
+ Note : for this to work with grass6.3, in find_file function from core.py,
+ command should be (n flag removed because 6.4 specific):
+ s = read_command("g.findfile", element = element, file = name, mapset = mapset)
+ """
+ ##Test for input vector map.
+ testInput = grass.find_file(input, element = 'vector')
+ if testInput['fullname'] == '':
+ raise AutoKrigeError("vector map " + input + " not found in mapset search path.")
+ ##Test if output raster map already exists.
+ testOutput = grass.find_file(output, element = 'cell')
+ if testOutput['fullname'] != '':
+ if not grass.overwrite() is True:
+ raise AutoKrigeError("raster map " + output + " already exists in mapset search path. \n \
+Use the --o flag to overwrite.")
+ else:
+ self.printMessage("raster map " + output + " will be overwritten.", type = 'warning')
+ #print "Warning: raster map " + output + " will be overwritten."
+ ##Test also variance raster.
+ if testVarianceRast is True:
+ testVarOutput = grass.find_file(output + '_var', element = 'cell')
+ if testVarOutput['fullname'] != '':
+ if not grass.overwrite() is True:
+ raise AutoKrigeError("raster map " + output + " already exists in mapset search path. \n \
+Use the --o flag to overwrite.")
+ else:
+ self.printMessage("raster map " + output + '_var' + " will be overwritten.", type = 'warning')
+ #print "Warning: raster map " + output + '_var' + " will be overwritten."
+
+
+ def getGridCellSize(self, input, nbcell, adjustRegionToSites = True):
+ """!Define kriged grid cell size.
+ Raster resolution but also computation time depends on it.
+ We take region resolution as cell size but we can fix this resolution with the nbcell parameter.
+ Only one value is needed because the R script use square cells."""
+ if adjustRegionToSites is True:
+ grass.run_command("g.region", vect = input)
+ if nbcell is not None:
+ self.fixRegionResFromNumberOfCells(nbcell)
+ regionParams = grass.region()
+ cellsize = float(regionParams['nsres']) if float(regionParams['nsres']) > float(regionParams['ewres']) \
+ else float(regionParams['ewres'])
+ return cellsize
+
+ def fixRegionResFromNumberOfCells(self, nbcell=100):
+ """!Adjust one of the two region dimensions so that we have a maximum of 'nbcell' cells in both."""
+ regionParams = grass.region()
+ nsResForGivenNbCell = (float(regionParams['n'])- float(regionParams['s'])) / float(nbcell)
+ ewResForGivenNbCell = (float(regionParams['e'])- float(regionParams['w'])) / float(nbcell)
+ tmpRegionRes = nsResForGivenNbCell if nsResForGivenNbCell > ewResForGivenNbCell else ewResForGivenNbCell
+ grass.use_temp_region
+ grass.run_command("g.region", res = tmpRegionRes)
+
+ def finalizeOutput(self):
+ """!Final operations after successful kriging.
We don't want to stop execution if an error occurs here.
"""
try:
##convert plot to png
- self.__execShellCommand('convert -alpha off Rplots.pdf autokrige_result.png')
+ call('convert -alpha off Rplots.pdf autokrige_result.png', shell = True)
grass.try_remove('Rplots.pdf')
##apply colormap
- grass.run_command("r.colors", map = output, rules = colormap, quiet = True)
+ grass.run_command("r.colors", map = self.output, rules = self.colormap, quiet = True)
if self.varianceFlag is True:
- grass.run_command("r.colors", map = output + '_var', rules = colormap, quiet = True)
+ grass.run_command("r.colors", map = self.output + '_var', rules = self.colormap, quiet = True)
except:
pass
-
- def __execShellCommand(self, command, stderrRedirection=False, logFile=False):
- """General purpose function, maybe should be added in some way to core.py """
- if logFile is not False:
- command = command + ' >> ' + logFile
- if stderrRedirection is True:
- command = command + " 2>&1"
- p = Popen(command, shell=True, stdout=PIPE)
- retcode = p.wait()
- com = p.communicate()
- if retcode == 0:
- return com[0]
- else:
- errorMessage = ''
- for elem in com:
- errorMessage += str(elem) + "\n"
- raise AutoKrigeError(errorMessage)
-
+
def runAutoKrige(self):
- """Autokrige class public method"""
+ """!Performs kriging process"""
##1)Necessary checks
- self.__checkLayers(self.input, self.output)
+ self.checkLayers(self.input, self.output, self.varianceFlag)
##2)Adjust interpolation resolution
- cellsize = self.__getGridCellSize(self.input, self.nbcell)
+ cellsize = self.getGridCellSize(self.input, self.nbcell, self.regionFlag)
##3)Execute R from GRASS
self.RscriptFile = self.__writeRScript()
##spgrass6 cause : column name in R has only the 10 first characters of the original column name
##may change with future versions of spgrass6
Rcolumnname = self.column[0:10]
writeVarRast = 'T' if self.varianceFlag is True else 'F'
- print "RGrass is working..."
+ self.printMessage("RGrass is working...", type = 'info')
+ #print "RGrass is working..."
##ordinary kriging for now
predictors = '1'
autoKrigeCommand = 'R --vanilla --slave --args ' + self.input + ' ' + Rcolumnname + ' ' + \
@@ -393,9 +412,9 @@
+ ' ' + self.nmax + ' ' + self.maxdist + ' ' \
+ writeVarRast + ' ' + predictors + ' < "' + self.RscriptFile + '"'
logFileArg = self.logfile if self.logROutput is True else False
- self.__execShellCommand(autoKrigeCommand, stderrRedirection=True, logFile=logFileArg)
- ##4)Finalize output
- self.__finalize(self.output, self.colormap)
+ self.__execRCommand(autoKrigeCommand, logFile=logFileArg)
+ ##4)Final steps
+ self.finalizeOutput()
class AutoKrigeError(Exception):
"""Errors specific to Autokrige class"""
@@ -411,15 +430,16 @@
autoKrige = AutoKrige(options, flags)
autoKrige.runAutoKrige()
except AutoKrigeError, e1:
- print >> sys.stderr, "Error \n:", e1.message
+ autoKrige.printMessage(e1.message, type = 'error')
exitStatus = 1
except:
- errorMessage = "Unexpected error \n:"
- print >> sys.stderr, errorMessage
- traceback.print_exc()
+ exceptionType, exceptionValue, exceptionTraceback = sys.exc_info()
+ errorMessage = "Unexpected error \n:" + \
+ repr(traceback.format_exception(exceptionType, exceptionValue, exceptionTraceback))
+ autoKrige.printMessage(errorMessage, type = 'error')
exitStatus = 1
else:
- print "Done"
+ autoKrige.printMessage("Done", type = 'info')
finally:
grass.try_remove(autoKrige.RscriptFile)
sys.exit(exitStatus)
More information about the grass-commit
mailing list