[GRASS-SVN] r41730 - grass/branches/develbranch_6/scripts/v.krige
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 5 11:34:53 EDT 2010
Author: martinl
Date: 2010-04-05 11:34:52 -0400 (Mon, 05 Apr 2010)
New Revision: 41730
Added:
grass/branches/develbranch_6/scripts/v.krige/v.krige.py
Removed:
grass/branches/develbranch_6/scripts/v.krige/v.krige
Modified:
grass/branches/develbranch_6/scripts/v.krige/Makefile
Log:
rename script
Modified: grass/branches/develbranch_6/scripts/v.krige/Makefile
===================================================================
--- grass/branches/develbranch_6/scripts/v.krige/Makefile 2010-04-05 13:43:20 UTC (rev 41729)
+++ grass/branches/develbranch_6/scripts/v.krige/Makefile 2010-04-05 15:34:52 UTC (rev 41730)
@@ -1,7 +1,7 @@
MODULE_TOPDIR = ../..
-PGM = v.krige
+PGM = v.krige.py
include $(MODULE_TOPDIR)/include/Make/Script.make
-default: script
+default:
Deleted: grass/branches/develbranch_6/scripts/v.krige/v.krige
===================================================================
--- grass/branches/develbranch_6/scripts/v.krige/v.krige 2010-04-05 13:43:20 UTC (rev 41729)
+++ grass/branches/develbranch_6/scripts/v.krige/v.krige 2010-04-05 15:34:52 UTC (rev 41730)
@@ -1,421 +0,0 @@
-#!/usr/bin/env python
-"""
-MODULE: v.krige
-
-AUTHOR(S): Anne Ghisla <a.ghisla AT gmail.com>
-
-PURPOSE: Performs ordinary or block kriging
-
-DEPENDS: R 2.x, packages gstat, maptools and spgrass6, optional: automap
-
-COPYRIGHT: (C) 2009 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.
-"""
-
-## g.parser informations
-
-#%module
-#% description: Performs ordinary or block kriging.
-#% keywords: vector, interpolation, kriging
-#%end
-
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,vector,vector
-#% description: Name of point vector map containing sample data
-#% required: yes
-#%end
-#%option
-#% key: column
-#% type: string
-#% gisprompt: old,dbcolumn,dbcolumn
-#% description: Name of attribute column with numerical value to be interpolated
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% label: Name for output raster map
-#% description: If omitted, will be <input name>_kriging
-#% required : no
-#%end
-#%option
-#% key: package
-#% type: string
-#% options: gstat
-#% answer: gstat
-#% description: R package to use
-#% required: no
-#%end
-#%option
-#% key: model
-#% type: string
-#% options: Exp,Sph,Gau,Mat,Lin
-#% multiple: yes
-#% label: Variogram model(s)
-#% description: Leave empty to test all models (requires automap)
-#% required: no
-#%end
-#%option
-#% key: block
-#% type: integer
-#% multiple: no
-#% label: Block size (square block)
-#% description: Block size. Used by block kriging.
-#% required: no
-#%end
-#%option
-#% key: range
-#% type: integer
-#% label: Range value
-#% description: Automatically fixed if not set
-#% required : no
-#%end
-#%option
-#% key: nugget
-#% type: integer
-#% label: Nugget value
-#% description: Automatically fixed if not set
-#% required : no
-#%end
-#%option
-#% key: sill
-#% type: integer
-#% label: Sill value
-#% description: Automatically fixed if not set
-#% required : no
-#%end
-#%option
-#% key: output_var
-#% type: string
-#% gisprompt: new,cell,raster
-#% label: Name for output variance raster map
-#% description: If omitted, will be <input name>_kriging_var
-#% required : no
-#%end
-
-import os, sys
-from tempfile import gettempdir
-import time
-import thread
-
-GUIModulesPath = os.path.join(os.getenv("GISBASE"), "etc", "wxpython", "gui_modules")
-sys.path.append(GUIModulesPath)
-#GUIPath = os.path.join(os.getenv("GISBASE"), "etc", "wxpython")
-#sys.path.append(GUIPath)
-#
-#import globalvar
-#if not os.getenv("GRASS_WXBUNDLED"):
-# globalvar.CheckForWx()
-#import gselect
-#import goutput
-#import menuform
-#from preferences import globalSettings as UserSettings
-##import help
-#
-#import wx
-#import wx.lib.flatnotebook as FN
-##import wx.lib.plot as plot # for plotting the variogram.
-
-### i18N
-import gettext
-gettext.install('grasswxpy', os.path.join(os.getenv("GISBASE"), 'locale'), unicode = True)
-
-### dependencies to be checked once, as they are quite time-consuming. cfr. grass.parser.
-# GRASS binding
-try:
- import grass.script as grass
-except ImportError:
- sys.exit(_("No GRASS-python library found."))
-
-# move other checks in functions, as R?
-
-# globals
-
-Region = grass.region()
-Command = None
-InputData = None
-Variogram = None
-VariogramFunction = None
-robjects = None
-rinterface = None
-
-#classes in alphabetical order. methods in logical order :)
-
-class Controller:
- """ Executes analysis. For the moment, only with gstat functions."""
-
- def ImportMap(self, map, column):
- """ Imports GRASS map as SpatialPointsDataFrame and adds x/y columns to attribute table.
- Checks for NULL values in the provided column and exits if they are present."""
-
- #@NOTE: new way with R - as it doesn't alter original data
- Rpointmap = robjects.r.readVECT6(map, type = 'point')
- # checks if x,y columns are present in dataframe. If they do are present, but with different names,
- # they'll be duplicated.
- if "x" not in robjects.r.names(Rpointmap):
- # extract coordinates with S4 method
- coordinatesPreDF = robjects.r['as.data.frame'](robjects.r.coordinates(Rpointmap))
- coordinatesDF = robjects.r['data.frame'](x = coordinatesPreDF.r['coords.x1'][0],
- y = coordinatesPreDF.r['coords.x2'][0])
- # match coordinates with data slot of SpatialPointsDataFrame - maptools function
- # match is done on row.names
- Rpointmap = robjects.r.spCbind(Rpointmap, coordinatesDF)
-
- # GRASS checks for null values in the chosen column. R can hardly handle column as a variable,
- # looks for a hardcoded string.
- cols = grass.vector_columns(map=map, layer=1)
- nulls = int(grass.parse_command('v.univar',
- map=map,
- column=column,
- type='point',
- parse = (grass.parse_key_val,
- {'sep':': '}
- )
- )['number of NULL attributes'])
- if nulls > 0:
- grass.fatal(_("%d NULL value(s) in the selected column - unable to perform kriging.") % nulls)
- return Rpointmap
-
- def CreateGrid(self, inputdata):
- Grid = robjects.r.gmeta2grd()
-
- # addition of coordinates columns into dataframe.
- coordinatesDF = robjects.r['as.data.frame'](robjects.r.coordinates(Grid))
- data = robjects.r['data.frame'](x = coordinatesDF.r['s1'][0],
- y = coordinatesDF.r['s2'][0],
- k = robjects.r.rep(1, Region['cols']*Region['rows']))
- GridPredicted = robjects.r.SpatialGridDataFrame(Grid,
- data,
- proj4string = robjects.r.CRS(robjects.r.proj4string(inputdata)))
- return GridPredicted
-
- def ComposeFormula(self, column, isblock, inputdata):
- if isblock is True:
- predictor = 'x+y'
- else:
- predictor = 1
- Formula = robjects.r['as.formula'](robjects.r.paste(column, "~", predictor))
- #print Formula
- return Formula
-
- def FitVariogram(self, formula, inputdata, sill, nugget, range, model = ''):
- """ Fits variogram either automagically either specifying all parameters.
- Returns a list containing data and model variograms. """
-
- Variograms = {}
-
- if model is '':
- robjects.r.require('automap')
- DottedParams = {}
- #print (nugget.r_repr(), sill, range)
- DottedParams['fix.values'] = robjects.r.c(nugget, range, sill)
-
- VariogramModel = robjects.r.autofitVariogram(formula, inputdata, **DottedParams)
- #print robjects.r.warnings()
- Variograms['datavariogram'] = VariogramModel.r['exp_var'][0]
- Variograms['variogrammodel'] = VariogramModel.r['var_model'][0]
- # obtain the model name. *Too* complicated to get the string instead of level, unlike R does.
- VariogramAsDF = robjects.r['as.data.frame'](VariogramModel.r['var_model'][0]) # force conversion
- ModelDF = VariogramAsDF.r['model'][0]
- Variograms['model'] = robjects.r.levels(ModelDF).subset(ModelDF[1])[0]
- else:
- DataVariogram = robjects.r['variogram'](formula, inputdata)
- VariogramModel = robjects.r['fit.variogram'](DataVariogram,
- model = robjects.r.vgm(psill = sill,
- model = model,
- nugget = nugget,
- range = range))
- Variograms['datavariogram'] = DataVariogram
- Variograms['variogrammodel'] = VariogramModel
- Variograms['model'] = model
- return Variograms
-
- def DoKriging(self, formula, inputdata, grid, model, block):
- DottedParams = {'debug.level': -1} # let krige() print percentage status
- if block is not '': #@FIXME(anne): but it's a string!! and krige accepts it!!
- DottedParams['block'] = block
- #print DottedParams
- KrigingResult = robjects.r.krige(formula, inputdata, grid, model, **DottedParams)
- return KrigingResult
-
- def ExportMap(self, map, column, name, overwrite, command, variograms):
- # add kriging parameters to raster map history
- robjects.r.writeRAST6(map, vname = name, zcol = column, overwrite = overwrite)
- grass.run_command('r.support',
- map = name,
- title = 'Kriging output',
- history = 'Issued from command v.krige ' + command)
- if command.find('model') is -1: # if the command has no model option, add automap chosen model
- grass.run_command('r.support',
- map = name,
- history = 'Model chosen by automatic fitting: ' + variograms['model'])
-
- def Run(self, input, column, output, package, sill, nugget, range, logger, \
- overwrite, model, block, output_var, command, **kwargs):
- """ Wrapper for all functions above. """
- logger.message(_("Processing %d cells. Computing time raises exponentially with resolution." % Region['cells']))
- logger.message(_("Importing data..."))
- if globals()["InputData"] is None:
- globals()["InputData"] = self.ImportMap(input, column)
- # and from here over, InputData refers to the global variable
- #print(robjects.r.slot(InputData, 'data').names)
- logger.message(_("Data successfully imported."))
-
- GridPredicted = self.CreateGrid(InputData)
-
- logger.message(_("Fitting variogram..."))
- isblock = block is not ''
- Formula = self.ComposeFormula(column, isblock, InputData)
- if globals()["Variogram"] is None:
- globals()["Variogram"] = self.FitVariogram(Formula,
- InputData,
- model = model,
- sill = sill,
- nugget = nugget,
- range = range)
- logger.message(_("Variogram fitted."))
-
- logger.message(_("Kriging..."))
- KrigingResult = self.DoKriging(Formula, InputData, GridPredicted, Variogram['variogrammodel'], block) # using global ones
- logger.message(_("Kriging performed."))
-
- self.ExportMap(map = KrigingResult,
- column='var1.pred',
- name = output,
- overwrite = overwrite,
- command = command,
- variograms = Variogram)
- if output_var is not '':
- self.ExportMap(map = KrigingResult,
- column='var1.var',
- name = output_var,
- overwrite = overwrite,
- command = command,
- variograms = Variogram)
-
-def main(argv = None):
- """ Main. Calls either GUI or CLI, depending on arguments provided. """
- #@FIXME: solve this double ifelse. the control should not be done twice.
-
- controller = Controller()
-
- if argv is None:
- importR()
- argv = sys.argv[1:] #stripping first item, the full name of this script
- # wxGUI call.
- import v_krige_wxGUI
- import wx
-
- app = wx.App()
- KrigingFrame = v_krige_wxGUI.KrigingModule(parent = None,
- Rinstance = robjects,
- controller = controller)
- KrigingFrame.Centre()
- KrigingFrame.Show()
- app.MainLoop()
-
- else:
- #CLI
- options, flags = argv
- #@TODO: Work on verbosity. Sometimes it's too verbose (R), sometimes not enough.
- if grass.find_file(options['input'], element = 'vector')['fullname'] is '':
- grass.fatal(_("option: <input>: Vector map not found.")) #TODO cosmetics, insert real map name
-
- #@TODO: elaborate input string, if contains mapset or not.. thanks again to Bob for testing on 64bit.
-
- # create output map name, if not specified
- if options['output'] is '':
- try: # to strip mapset name from fullname. Ugh.
- options['input'] = options['input'].split("@")[0]
- except:
- pass
- options['output'] = options['input'] + '_kriging'
-
- # check for output map with same name. g.parser can't handle this, afaik.
- if grass.find_file(options['output'], element = 'cell')['fullname'] and os.getenv("GRASS_OVERWRITE") == None:
- grass.fatal(_("option: <output>: Raster map already exists."))
- if options['output_var'] is not '' and (grass.find_file(options['output_var'], element = 'cell')['fullname'] and os.getenv("GRASS_OVERWRITE") == None):
- grass.fatal(_("option: <output>: Variance raster map already exists."))
-
- importR()
- if options['model'] is '':
- try:
- robjects.r.require("automap")
- except ImportError, e:
- grass.fatal(_("R package automap is missing, no variogram autofit available."))
- else:
- if options['sill'] is '' or options['nugget'] is '' or options['range'] is '':
- grass.fatal(_("You have specified model, but forgot at least one of sill, nugget and range."))
-
- #@TODO: let GRASS remount its commandstring. Until then, keep that 4 lines below.
- #print grass.write_command(argv)
- command = ""
- notnulloptions = {}
- for k, v in options.items():
- if v is not '':
- notnulloptions[k] = v
- command = command.join("%s=%s " % (k, v) for k, v in notnulloptions.items())
- #print command
-
- # re-cast integers from strings, as parser() cast everything to string.
- for each in ("sill","nugget","range"):
- if options[each] is not '':
- options[each] = int(options[each])
- else:
- options[each] = robjects.r('''NA''')
-
- #controller = Controller()
- controller.Run(input = options['input'],
- column = options['column'],
- output = options['output'],
- overwrite = os.getenv("GRASS_OVERWRITE") == 1,
- package = options['package'],
- model = options['model'],
- block = options['block'],
- sill = options['sill'],
- nugget = options['nugget'],
- range = options['range'],
- output_var = options['output_var'],
- command = command,
- logger = grass)
-
-def importR():
- # R
- # unuseful since rpy2 will complain adequately.
- # try:
- # #@FIXME: in Windows, it launches R terminal
- # grass.find_program('R')
- # except:
- # sys.exit(_("R is not installed. Install it and re-run, or modify environment variables."))
-
- # rpy2
- global robjects
- global rinterface
- grass.message(_('Loading packages, please wait...'))
- try:
- import rpy2.robjects as robjects
- import rpy2.rinterface as rinterface #to speed up kriging? for plots.
- haveRpy2 = True
- except ImportError:
- print >> sys.stderr, _("Rpy2 not found. Please install it and re-run.") # ok for other OSes?
- haveRpy2 = False
- if not haveRpy2:
- sys.exit(1)
-
- # R packages check.
- # @FIXME: it leaves a Rtmpxxxx folder into the make tempfolder and causes make complain. [markus]
- for each in ["gstat", "spgrass6", "maptools"]:
- if not robjects.r.require(each, quietly = True)[0]:
- sys.exit(_("R package '%s' is missing. Install it and re-run v.krige.") % each)
-
-if __name__ == '__main__':
- if len(sys.argv) > 1:
- sys.exit(main(argv = grass.parser()))
- else:
- main()
Copied: grass/branches/develbranch_6/scripts/v.krige/v.krige.py (from rev 41705, grass/branches/develbranch_6/scripts/v.krige/v.krige)
===================================================================
--- grass/branches/develbranch_6/scripts/v.krige/v.krige.py (rev 0)
+++ grass/branches/develbranch_6/scripts/v.krige/v.krige.py 2010-04-05 15:34:52 UTC (rev 41730)
@@ -0,0 +1,421 @@
+#!/usr/bin/env python
+"""
+MODULE: v.krige
+
+AUTHOR(S): Anne Ghisla <a.ghisla AT gmail.com>
+
+PURPOSE: Performs ordinary or block kriging
+
+DEPENDS: R 2.x, packages gstat, maptools and spgrass6, optional: automap
+
+COPYRIGHT: (C) 2009 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.
+"""
+
+## g.parser informations
+
+#%module
+#% description: Performs ordinary or block kriging.
+#% keywords: vector, interpolation, kriging
+#%end
+
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Name of point vector map containing sample data
+#% required: yes
+#%end
+#%option
+#% key: column
+#% type: string
+#% gisprompt: old,dbcolumn,dbcolumn
+#% description: Name of attribute column with numerical value to be interpolated
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% label: Name for output raster map
+#% description: If omitted, will be <input name>_kriging
+#% required : no
+#%end
+#%option
+#% key: package
+#% type: string
+#% options: gstat
+#% answer: gstat
+#% description: R package to use
+#% required: no
+#%end
+#%option
+#% key: model
+#% type: string
+#% options: Exp,Sph,Gau,Mat,Lin
+#% multiple: yes
+#% label: Variogram model(s)
+#% description: Leave empty to test all models (requires automap)
+#% required: no
+#%end
+#%option
+#% key: block
+#% type: integer
+#% multiple: no
+#% label: Block size (square block)
+#% description: Block size. Used by block kriging.
+#% required: no
+#%end
+#%option
+#% key: range
+#% type: integer
+#% label: Range value
+#% description: Automatically fixed if not set
+#% required : no
+#%end
+#%option
+#% key: nugget
+#% type: integer
+#% label: Nugget value
+#% description: Automatically fixed if not set
+#% required : no
+#%end
+#%option
+#% key: sill
+#% type: integer
+#% label: Sill value
+#% description: Automatically fixed if not set
+#% required : no
+#%end
+#%option
+#% key: output_var
+#% type: string
+#% gisprompt: new,cell,raster
+#% label: Name for output variance raster map
+#% description: If omitted, will be <input name>_kriging_var
+#% required : no
+#%end
+
+import os, sys
+from tempfile import gettempdir
+import time
+import thread
+
+GUIModulesPath = os.path.join(os.getenv("GISBASE"), "etc", "wxpython", "gui_modules")
+sys.path.append(GUIModulesPath)
+#GUIPath = os.path.join(os.getenv("GISBASE"), "etc", "wxpython")
+#sys.path.append(GUIPath)
+#
+#import globalvar
+#if not os.getenv("GRASS_WXBUNDLED"):
+# globalvar.CheckForWx()
+#import gselect
+#import goutput
+#import menuform
+#from preferences import globalSettings as UserSettings
+##import help
+#
+#import wx
+#import wx.lib.flatnotebook as FN
+##import wx.lib.plot as plot # for plotting the variogram.
+
+### i18N
+import gettext
+gettext.install('grasswxpy', os.path.join(os.getenv("GISBASE"), 'locale'), unicode = True)
+
+### dependencies to be checked once, as they are quite time-consuming. cfr. grass.parser.
+# GRASS binding
+try:
+ import grass.script as grass
+except ImportError:
+ sys.exit(_("No GRASS-python library found."))
+
+# move other checks in functions, as R?
+
+# globals
+
+Region = grass.region()
+Command = None
+InputData = None
+Variogram = None
+VariogramFunction = None
+robjects = None
+rinterface = None
+
+#classes in alphabetical order. methods in logical order :)
+
+class Controller:
+ """ Executes analysis. For the moment, only with gstat functions."""
+
+ def ImportMap(self, map, column):
+ """ Imports GRASS map as SpatialPointsDataFrame and adds x/y columns to attribute table.
+ Checks for NULL values in the provided column and exits if they are present."""
+
+ #@NOTE: new way with R - as it doesn't alter original data
+ Rpointmap = robjects.r.readVECT6(map, type = 'point')
+ # checks if x,y columns are present in dataframe. If they do are present, but with different names,
+ # they'll be duplicated.
+ if "x" not in robjects.r.names(Rpointmap):
+ # extract coordinates with S4 method
+ coordinatesPreDF = robjects.r['as.data.frame'](robjects.r.coordinates(Rpointmap))
+ coordinatesDF = robjects.r['data.frame'](x = coordinatesPreDF.r['coords.x1'][0],
+ y = coordinatesPreDF.r['coords.x2'][0])
+ # match coordinates with data slot of SpatialPointsDataFrame - maptools function
+ # match is done on row.names
+ Rpointmap = robjects.r.spCbind(Rpointmap, coordinatesDF)
+
+ # GRASS checks for null values in the chosen column. R can hardly handle column as a variable,
+ # looks for a hardcoded string.
+ cols = grass.vector_columns(map=map, layer=1)
+ nulls = int(grass.parse_command('v.univar',
+ map=map,
+ column=column,
+ type='point',
+ parse = (grass.parse_key_val,
+ {'sep':': '}
+ )
+ )['number of NULL attributes'])
+ if nulls > 0:
+ grass.fatal(_("%d NULL value(s) in the selected column - unable to perform kriging.") % nulls)
+ return Rpointmap
+
+ def CreateGrid(self, inputdata):
+ Grid = robjects.r.gmeta2grd()
+
+ # addition of coordinates columns into dataframe.
+ coordinatesDF = robjects.r['as.data.frame'](robjects.r.coordinates(Grid))
+ data = robjects.r['data.frame'](x = coordinatesDF.r['s1'][0],
+ y = coordinatesDF.r['s2'][0],
+ k = robjects.r.rep(1, Region['cols']*Region['rows']))
+ GridPredicted = robjects.r.SpatialGridDataFrame(Grid,
+ data,
+ proj4string = robjects.r.CRS(robjects.r.proj4string(inputdata)))
+ return GridPredicted
+
+ def ComposeFormula(self, column, isblock, inputdata):
+ if isblock is True:
+ predictor = 'x+y'
+ else:
+ predictor = 1
+ Formula = robjects.r['as.formula'](robjects.r.paste(column, "~", predictor))
+ #print Formula
+ return Formula
+
+ def FitVariogram(self, formula, inputdata, sill, nugget, range, model = ''):
+ """ Fits variogram either automagically either specifying all parameters.
+ Returns a list containing data and model variograms. """
+
+ Variograms = {}
+
+ if model is '':
+ robjects.r.require('automap')
+ DottedParams = {}
+ #print (nugget.r_repr(), sill, range)
+ DottedParams['fix.values'] = robjects.r.c(nugget, range, sill)
+
+ VariogramModel = robjects.r.autofitVariogram(formula, inputdata, **DottedParams)
+ #print robjects.r.warnings()
+ Variograms['datavariogram'] = VariogramModel.r['exp_var'][0]
+ Variograms['variogrammodel'] = VariogramModel.r['var_model'][0]
+ # obtain the model name. *Too* complicated to get the string instead of level, unlike R does.
+ VariogramAsDF = robjects.r['as.data.frame'](VariogramModel.r['var_model'][0]) # force conversion
+ ModelDF = VariogramAsDF.r['model'][0]
+ Variograms['model'] = robjects.r.levels(ModelDF).subset(ModelDF[1])[0]
+ else:
+ DataVariogram = robjects.r['variogram'](formula, inputdata)
+ VariogramModel = robjects.r['fit.variogram'](DataVariogram,
+ model = robjects.r.vgm(psill = sill,
+ model = model,
+ nugget = nugget,
+ range = range))
+ Variograms['datavariogram'] = DataVariogram
+ Variograms['variogrammodel'] = VariogramModel
+ Variograms['model'] = model
+ return Variograms
+
+ def DoKriging(self, formula, inputdata, grid, model, block):
+ DottedParams = {'debug.level': -1} # let krige() print percentage status
+ if block is not '': #@FIXME(anne): but it's a string!! and krige accepts it!!
+ DottedParams['block'] = block
+ #print DottedParams
+ KrigingResult = robjects.r.krige(formula, inputdata, grid, model, **DottedParams)
+ return KrigingResult
+
+ def ExportMap(self, map, column, name, overwrite, command, variograms):
+ # add kriging parameters to raster map history
+ robjects.r.writeRAST6(map, vname = name, zcol = column, overwrite = overwrite)
+ grass.run_command('r.support',
+ map = name,
+ title = 'Kriging output',
+ history = 'Issued from command v.krige ' + command)
+ if command.find('model') is -1: # if the command has no model option, add automap chosen model
+ grass.run_command('r.support',
+ map = name,
+ history = 'Model chosen by automatic fitting: ' + variograms['model'])
+
+ def Run(self, input, column, output, package, sill, nugget, range, logger, \
+ overwrite, model, block, output_var, command, **kwargs):
+ """ Wrapper for all functions above. """
+ logger.message(_("Processing %d cells. Computing time raises exponentially with resolution." % Region['cells']))
+ logger.message(_("Importing data..."))
+ if globals()["InputData"] is None:
+ globals()["InputData"] = self.ImportMap(input, column)
+ # and from here over, InputData refers to the global variable
+ #print(robjects.r.slot(InputData, 'data').names)
+ logger.message(_("Data successfully imported."))
+
+ GridPredicted = self.CreateGrid(InputData)
+
+ logger.message(_("Fitting variogram..."))
+ isblock = block is not ''
+ Formula = self.ComposeFormula(column, isblock, InputData)
+ if globals()["Variogram"] is None:
+ globals()["Variogram"] = self.FitVariogram(Formula,
+ InputData,
+ model = model,
+ sill = sill,
+ nugget = nugget,
+ range = range)
+ logger.message(_("Variogram fitted."))
+
+ logger.message(_("Kriging..."))
+ KrigingResult = self.DoKriging(Formula, InputData, GridPredicted, Variogram['variogrammodel'], block) # using global ones
+ logger.message(_("Kriging performed."))
+
+ self.ExportMap(map = KrigingResult,
+ column='var1.pred',
+ name = output,
+ overwrite = overwrite,
+ command = command,
+ variograms = Variogram)
+ if output_var is not '':
+ self.ExportMap(map = KrigingResult,
+ column='var1.var',
+ name = output_var,
+ overwrite = overwrite,
+ command = command,
+ variograms = Variogram)
+
+def main(argv = None):
+ """ Main. Calls either GUI or CLI, depending on arguments provided. """
+ #@FIXME: solve this double ifelse. the control should not be done twice.
+
+ controller = Controller()
+
+ if argv is None:
+ importR()
+ argv = sys.argv[1:] #stripping first item, the full name of this script
+ # wxGUI call.
+ import v_krige_wxGUI
+ import wx
+
+ app = wx.App()
+ KrigingFrame = v_krige_wxGUI.KrigingModule(parent = None,
+ Rinstance = robjects,
+ controller = controller)
+ KrigingFrame.Centre()
+ KrigingFrame.Show()
+ app.MainLoop()
+
+ else:
+ #CLI
+ options, flags = argv
+ #@TODO: Work on verbosity. Sometimes it's too verbose (R), sometimes not enough.
+ if grass.find_file(options['input'], element = 'vector')['fullname'] is '':
+ grass.fatal(_("option: <input>: Vector map not found.")) #TODO cosmetics, insert real map name
+
+ #@TODO: elaborate input string, if contains mapset or not.. thanks again to Bob for testing on 64bit.
+
+ # create output map name, if not specified
+ if options['output'] is '':
+ try: # to strip mapset name from fullname. Ugh.
+ options['input'] = options['input'].split("@")[0]
+ except:
+ pass
+ options['output'] = options['input'] + '_kriging'
+
+ # check for output map with same name. g.parser can't handle this, afaik.
+ if grass.find_file(options['output'], element = 'cell')['fullname'] and os.getenv("GRASS_OVERWRITE") == None:
+ grass.fatal(_("option: <output>: Raster map already exists."))
+ if options['output_var'] is not '' and (grass.find_file(options['output_var'], element = 'cell')['fullname'] and os.getenv("GRASS_OVERWRITE") == None):
+ grass.fatal(_("option: <output>: Variance raster map already exists."))
+
+ importR()
+ if options['model'] is '':
+ try:
+ robjects.r.require("automap")
+ except ImportError, e:
+ grass.fatal(_("R package automap is missing, no variogram autofit available."))
+ else:
+ if options['sill'] is '' or options['nugget'] is '' or options['range'] is '':
+ grass.fatal(_("You have specified model, but forgot at least one of sill, nugget and range."))
+
+ #@TODO: let GRASS remount its commandstring. Until then, keep that 4 lines below.
+ #print grass.write_command(argv)
+ command = ""
+ notnulloptions = {}
+ for k, v in options.items():
+ if v is not '':
+ notnulloptions[k] = v
+ command = command.join("%s=%s " % (k, v) for k, v in notnulloptions.items())
+ #print command
+
+ # re-cast integers from strings, as parser() cast everything to string.
+ for each in ("sill","nugget","range"):
+ if options[each] is not '':
+ options[each] = int(options[each])
+ else:
+ options[each] = robjects.r('''NA''')
+
+ #controller = Controller()
+ controller.Run(input = options['input'],
+ column = options['column'],
+ output = options['output'],
+ overwrite = os.getenv("GRASS_OVERWRITE") == 1,
+ package = options['package'],
+ model = options['model'],
+ block = options['block'],
+ sill = options['sill'],
+ nugget = options['nugget'],
+ range = options['range'],
+ output_var = options['output_var'],
+ command = command,
+ logger = grass)
+
+def importR():
+ # R
+ # unuseful since rpy2 will complain adequately.
+ # try:
+ # #@FIXME: in Windows, it launches R terminal
+ # grass.find_program('R')
+ # except:
+ # sys.exit(_("R is not installed. Install it and re-run, or modify environment variables."))
+
+ # rpy2
+ global robjects
+ global rinterface
+ grass.message(_('Loading packages, please wait...'))
+ try:
+ import rpy2.robjects as robjects
+ import rpy2.rinterface as rinterface #to speed up kriging? for plots.
+ haveRpy2 = True
+ except ImportError:
+ print >> sys.stderr, _("Rpy2 not found. Please install it and re-run.") # ok for other OSes?
+ haveRpy2 = False
+ if not haveRpy2:
+ sys.exit(1)
+
+ # R packages check.
+ # @FIXME: it leaves a Rtmpxxxx folder into the make tempfolder and causes make complain. [markus]
+ for each in ["gstat", "spgrass6", "maptools"]:
+ if not robjects.r.require(each, quietly = True)[0]:
+ sys.exit(_("R package '%s' is missing. Install it and re-run v.krige.") % each)
+
+if __name__ == '__main__':
+ if len(sys.argv) > 1:
+ sys.exit(main(argv = grass.parser()))
+ else:
+ main()
More information about the grass-commit
mailing list