Changeset 400

Timestamp:

05/31/09 20:35:21 (4 years ago)

Author:

jjr8

Message:

Checkin of unfinished code so that I can rebuild the Arc 9.3 toolbox for testing purposes.

Location:

MGET/Branches/Jason/PythonPackage/src/GeoEco

Files:

• R.py (modified) (4 diffs)
• Statistics/Exploratory.py (modified) (3 diffs)
• Statistics/Modeling.py (modified) (13 diffs)

MGET/Branches/Jason/PythonPackage/src/GeoEco/R.py

@@ -48 +48 @@
 
     @classmethod
-    def Evaluate(cls, statements, variableNames=None, variableValues=None, tables=None, dataFrameNamesForTables=None, xColumnName=u'x', yColumnName=u'y', zColumnName=u'z', mColumnName=u'm', stringsAsFactors=True, timeZone=None):
+    def Evaluate(cls, statements, clearGlobalEnv=True, variableNames=None, variableValues=None, tables=None, dataFrameNamesForTables=None, xColumnName=u'x', yColumnName=u'y', zColumnName=u'z', mColumnName=u'm', stringsAsFactors=True, timeZone=None):
         cls.__doc__.Obj.ValidateMethodInvocation()
 
@@ -138 +138 @@
 
     @classmethod
-    def EvaluateFile(cls, path, variableNames=None, variableValues=None, tables=None, dataFrameNamesForTables=None, xColumnName=u'x', yColumnName=u'y', zColumnName=u'z', mColumnName=u'm', stringsAsFactors=True, timeZone=None):
+    def EvaluateFile(cls, path, clearGlobalEnv=True, variableNames=None, variableValues=None, tables=None, dataFrameNamesForTables=None, xColumnName=u'x', yColumnName=u'y', zColumnName=u'z', mColumnName=u'm', stringsAsFactors=True, timeZone=None):
         cls.__doc__.Obj.ValidateMethodInvocation()
 
@@ -1838 +1838 @@
     arcGISDisplayName=_(u'R statements'))
 
+AddArgumentMetadata(R.Evaluate, u'clearGlobalEnv',
+    typeMetadata=BooleanTypeMetadata(),
+    description=_(
+u"""If True, the default, all objects will be removed from the R
+global environment after the statements are evaluated. If False, the R
+global environment will be left unchanged.
+
+This parameter is intended as a convenience to help minimize the
+memory utilization of this tool when it is run "in process". This
+parameter has no effect when this tool is run "out of process", as
+will occur when you run this tool from an ArcGIS geoprocessing model
+or script with ArcGIS version 9.2 or earlier. If you are in that
+situation, you should ignore this parameter.
+
+When this tool is run for the first time, it instantiates the R
+interpreter in the process that invoked it. If you are using ArcGIS,
+this process is typically an instance of ArcCatalog or ArcMap. When
+you run the tool again, it reuses the same R interpreter. Thus, if you
+specified False for this parameter the first time you ran this tool,
+the R variables defined by your statements will remain defined and be
+available the second time you run this tool. That can be good or bad,
+depending on what you want.
+
+The main benefit is that your second batch of statements can access
+results obtained by the first batch. This allows you to interleave
+ArcGIS operations between batches of R statements. The main downside
+is that the objects allocated by the first batch occupy memory.
+Because ArcGIS processes can only utilize about 1.5 GB of memory
+before crashing, leaving the R objects in memory may limit your
+subsequent geoprocessing operations.
+
+Because we anticipate that most users will not want to pass results
+between subsequent executions of this tool, we set the default value
+of this parameter to True, to minimize memory utilization. If you do
+want to pass results between subsequent executions, set this parameter
+to False.
+
+When this parameter is true, it cleans out the R global environment by
+executing the following R statement::
+
+    rm(list = ls())
+"""),
+    arcGISDisplayName=_(u'Clear R global environment after evaluating statements'))
+
 AddArgumentMetadata(R.Evaluate, u'variableNames',
     typeMetadata=ListTypeMetadata(elementType=UnicodeStringTypeMetadata(minLength=1), canBeNone=True),
@@ -2137 +2181 @@
     arcGISDisplayName=_(u'Text file'))
 
+CopyArgumentMetadata(R.Evaluate, u'clearGlobalEnv', R.EvaluateFile, u'clearGlobalEnv')
 CopyArgumentMetadata(R.Evaluate, u'variableNames', R.EvaluateFile, u'variableNames')
 CopyArgumentMetadata(R.Evaluate, u'variableValues', R.EvaluateFile, u'variableValues')

MGET/Branches/Jason/PythonPackage/src/GeoEco/Statistics/Exploratory.py

@@ -239 +239 @@
     arcGISDisplayName=_(u'Scatterplot Matrix for Table'),
     arcGISToolCategory=_(u'Statistics\\Explore Data'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 1), RDependency(2, 5, 0)])
 
@@ -459 +458 @@
     arcGISDisplayName=_(u'Density Histogram for Field'),
     arcGISToolCategory=_(u'Statistics\\Explore Data'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 1), RDependency(2, 5, 0)])
 
@@ -553 +551 @@
     arcGISDisplayName=_(u'Density Histogram for Point Coordinate'),
     arcGISToolCategory=_(u'Statistics\\Explore Data'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 1), RDependency(2, 5, 0)])
 

MGET/Branches/Jason/PythonPackage/src/GeoEco/Statistics/Modeling.py

@@ -775 +775 @@
     if templateRaster is None:
         templateRaster = inputPredictorRasters[0]
-    else:
-        d = gp.Describe(templateRaster)
-        if d.SpatialReference.Name is None or d.SpatialReference.Name.lower() == u'unknown':
-            Logger.RaiseException(ValueError(_(u'The template raster %(raster)s does not have a coordinate system defined. Please define the coordinate system for this raster and try again. (Use the ArcGIS Define Projection tool.)') % {u'raster': templateRaster}))
+
+    templateDescribe = gp.Describe(templateRaster)
+    if templateDescribe.SpatialReference.Name is None or templateDescribe.SpatialReference.Name.lower() == u'unknown':
+        Logger.RaiseException(ValueError(_(u'The template raster %(raster)s does not have a coordinate system defined. Please define the coordinate system for this raster and try again. (Use the ArcGIS Define Projection tool.)') % {u'raster': templateRaster}))
+
+    coordinateSystem = gp.CreateSpatialReference(templateDescribe.SpatialReference).split(';')[0]
+    cellSize = templateDescribe.MeanCellWidth
+    extent = templateDescribe.Extent
 
     # Load the fitted model.
@@ -868 +872 @@
             # First prepare the input rasters for the prediction.
 
-            _PreparePredictorRasters(r, gp, tempDir, variableNames, inputPredictorRasters, templateRaster, resamplingTechniques)
+            _PreparePredictorRasters(r, gp, tempDir, variableNames, inputPredictorRasters, templateRaster, templateDescribe, coordinateSystem, cellSize, extent, resamplingTechniques)
 
             # Do the prediction. The prediction rasters are
@@ -945 +949 @@
         r('if (exists("fam20985982305")) rm("fam20985982305", pos=globalenv())')
 
-def _PreparePredictorRasters(r, gp, tempDir, variableNames, inputPredictorRasters, templateRaster, resamplingTechniques):
-
-    # Build a dictionary that maps predictor variable names in
-    # the model to rasters provided by the caller.
+def _PreparePredictorRasters(r, gp, tempDir, variableNames, inputPredictorRasters, templateRaster, templateDescribe, coordinateSystem, cellSize, extent, resamplingTechniques):
+
+    # Build a dictionary that maps predictor variable names in the
+    # model to rasters provided by the caller.
 
     rastersForPredictors = {}
@@ -965 +969 @@
             rastersForPredictors[variable] = inputPredictorRasters[variableNames.index(variable)]
 
-    # Generate the X and Y coordinate rasters, if needed
-    # and possible.
+    # Generate the X and Y coordinate rasters, if needed and possible.
 
     if createXRaster:
@@ -992 +995 @@
         rastersForPredictors[r['yVar']] = yRaster
 
-    # Check each predictor raster. If it has a different
-    # coordinate system, extent, or cell size than that
-    # specified by the caller for the output rasters,
+    # Check each predictor raster. If it has a different coordinate
+    # system, extent, or cell size than that of the template raster,
     # project it in the temp directory.
+    #
+    # Note that ArcGIS does not always seem to represent the cell size
+    # and extent of rasters at the fullest possible precision, so we
+    # perform comparisons at less than full precision.
 
     Logger.Info('Checking projections, extents, and cell sizes of predictor rasters and reprojecting and clipping as needed...')
+
+    from GeoEco.Types import EnvelopeTypeMetadata
+    [templateLeft, templateBottom, templateRight, templateTop] = EnvelopeTypeMetadata.ParseCoordinatesFromString(extent)
 
     vars = copy.deepcopy(rastersForPredictors.keys())           # I'm going to modify the dictionary in the loop below; not sure if deepcopy is needed, but playing it safe.
@@ -1004 +1013 @@
 
     for i in range(len(vars)):
-        needToProject = True
+        if rasters[i] == templateRaster:
+            continue
+        
+        needToProject = False
+        needToClip = False
         d = gp.Describe(rasters[i])
         cs = gp.CreateSpatialReference_management(d.SpatialReference).split(u';')[0]
-        
+
         if cs != coordinateSystem:
             Logger.Debug(_(u'%(raster)s must be projected because its coordinate system (%(cs1)s) does not match the output coordinate system (%(cs2)s.') % {u'raster': rasters[i], u'cs1': cs, u'cs2': coordinateSystem})
-        elif d.Extent != extent:
-            Logger.Debug(_(u'%(raster)s must be projected because its extent (%(e1)s) does not match the output extent (%(e2)s).') % {u'raster': rasters[i], u'e1': d.Extent, u'e2': extent})
-        elif d.MeanCellWidth != cellSize:
-            Logger.Debug(_(u'%(raster)s must be projected because its cell size (%(cs1)f) does not match the output cell size (%(cs2)f).') % {u'raster': rasters[i], u'cs1': d.MeanCellWidth, u'cs2': cellSize})
+            Logger.Info(_(u'Projecting and clipping %(raster)s because its coordinate system does not match that of the template raster.') % {u'raster': rasters[i]})
+            needToProject = True
+
         else:
-            needToProject = False
+            [rasterLeft, rasterBottom, rasterRight, rasterTop] = EnvelopeTypeMetadata.ParseCoordinatesFromString(d.Extent)
+
+            if rasterLeft > templateLeft:
+                Logger.RaiseException(_(u'The template raster is too large; its left edge falls outside the left edge of the predictor raster %(raster)s. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+            if rasterRight < templateRight:
+                Logger.RaiseException(_(u'The template raster is too large; its right edge falls outside the right edge of the predictor raster %(raster)s. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+            if rasterTop > templateTop:
+                Logger.RaiseException(_(u'The template raster is too large; its top edge falls outside the top edge of the predictor raster %(raster)s. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+            if rasterBottom < templateBottom:
+                Logger.RaiseException(_(u'The template raster is too large; its bottom edge falls outside the bottom edge of the predictor raster %(raster)s. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+                
+            elif abs(1.0 - cellSize / d.MeanCellWidth) > 0.00001:
+                Logger.Debug(_(u'%(raster)s must be projected because its cell size (%(cs1)f) does not match the output cell size (%(cs2)f).') % {u'raster': rasters[i], u'cs1': d.MeanCellWidth, u'cs2': cellSize})
+                Logger.Info(_(u'Projecting and clipping %(raster)s because its cell size does not match that of the template raster.') % {u'raster': rasters[i]})
+                needToProject = True
+
+            else:
+                deltaLeft = (float(extent.split()[0]) - float(d.Extent.split()[0])) / cellSize     # Difference, in cells, between the left edge of the template and this raster
+                deltaBottom = (float(extent.split()[1]) - float(d.Extent.split()[1])) / cellSize   # Difference, in cells, between the bottom edge of the template and this raster
+                deltaRight = (float(extent.split()[2]) - float(d.Extent.split()[2])) / cellSize    # Difference, in cells, between the right edge of the template and this raster
+                deltaTop = (float(extent.split()[3]) - float(d.Extent.split()[3])) / cellSize      # Difference, in cells, between the top edge of the template and this raster
+
+                if deltaLeft < -0.01:
+                    Logger.RaiseException(_(u'The template raster is too large. The left edge of the predictor raster %(raster)s is inside the left edge of the template raster. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+                if deltaRight > -0.01:
+                    Logger.RaiseException(_(u'The template raster is too large. The left edge of the predictor raster %(raster)s is inside the left edge of the template raster. Please reduce the extent of the template raster such that it does not exceed the extents of any of the predictor rasters.') % {u'raster': rasters[i]})
+
+                if deltaLeft < 0.01 and deltaBottom < 0.01 and deltaRight < 0.01 and deltaTop < 0.01:
+                    Logger.Debug(_(u'%(raster)s has the same coordinate system, cell size, and extent as the template raster.') % {u'raster': rasters[i]})
+
+                elif abs(deltaLeft - round(deltaLeft)) < 0.01 and abs(deltaBottom - round(deltaBottom)) < 0.01 and abs(deltaRight - round(deltaRight)) < 0.01 and abs(deltaTop - round(deltaTop)) < 0.01:
+                    Logger.Debug(_(u'%(raster)s must be clipped because it has the same coordinate system and cell size as the template raster, and is aligned on the same grid coordinates, but has a different extent.') % {u'raster': rasters[i]})
+                    Logger.Info(_(u'Clipping %(raster)s to the extent of the template raster.') % {u'raster': rasters[i]})
+                    needToClip = True
+
+                else:
+                    Logger.Debug(_(u'%(raster)s must be clipped because it has the same coordinate system and cell size as the template raster, and is aligned on the same grid coordinates, but has a different extent.') % {u'raster': rasters[i]})
+                    Logger.Info(_(u'Clipping %(raster)s to the extent of the template raster.') % {u'raster': rasters[i]})
+                    needToProject = True
+                    
+
+                if abs(offsetX) > 0.01 or abs(offsetY) > 0.01:
+                    if abs(offsetX - round(offsetX)) > 0.01 or abs(offsetX - round(offsetX)) > 0.01:
+                        Logger.Debug(_(u'%(raster)s must be projected because its cell size (%(cs1)f) does not match the output cell size (%(cs2)f).') % {u'raster': rasters[i], u'cs1': d.MeanCellWidth, u'cs2': cellSize})
+                        Logger.Info(_(u'Projecting and clipping %(raster)s because its cell size does not match that of the template raster.') % {u'raster': rasters[i]})
+                        needToProject = True
 
         if needToProject:
@@ -1142 +1199 @@
     arcGISDisplayName=_(u'Fit GLM'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Generalized Linear Models'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 1), RDependency(2, 5, 0)])
 
@@ -1567 +1623 @@
     arcGISDisplayName=_(u'Predict GLM From Rasters'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Generalized Linear Models'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 2), RDependency(2, 5, 0), RPackageDependency(u'rgdal')])
 
@@ -1822 +1877 @@
     arcGISDisplayName=_(u'Fit GAM'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Generalized Additive Models'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 1), RDependency(2, 5, 0)])
 
@@ -2005 +2059 @@
     arcGISDisplayName=_(u'Predict GAM From Rasters'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Generalized Additive Models'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 2), RDependency(2, 5, 0), RPackageDependency(u'rgdal')])
 
@@ -2038 +2091 @@
     arcGISDisplayName=_(u'Predict Bayesian Probabilities for GAM From Rasters'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Generalized Additive Models'),
-    isOutOfProcArcGISTool=True,
     dependencies=[ArcGISDependency(9, 2), RDependency(2, 5, 0), RPackageDependency(u'mgcv'), RPackageDependency(u'MASS'), RPackageDependency(u'rgdal')])
 
@@ -2157 +2209 @@
     arcGISDisplayName=_(u'Plot Performance of Binary Classification Model'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Evaluate Model Performance'),
-    isOutOfProcArcGISTool=True,
     dependencies=[RDependency(2, 5, 0), RPackageDependency(u'ROCR')])
 
@@ -2501 +2552 @@
     arcGISDisplayName=_(u'Plot ROC of Binary Classification Model'),
     arcGISToolCategory=_(u'Statistics\\Model Data\\Evaluate Model Performance'),
-    isOutOfProcArcGISTool=True,
     dependencies=[RDependency(2, 5, 0), RPackageDependency(u'ROCR')])