Changeset 877

Timestamp:

01/03/12 13:25:49 (17 months ago)

Author:

jjr8

Message:

Added Project Raster To Template tool.

Location:

MGET/Branches/Jason/PythonPackage/src/GeoEco

Files:

• Connectivity/CoralReefConnectivity.py (modified) (3 diffs)
• DataManagement/ArcGISRasters.py (modified) (6 diffs)
• SpatialAnalysis/Interpolation.py (modified) (4 diffs)

MGET/Branches/Jason/PythonPackage/src/GeoEco/Connectivity/CoralReefConnectivity.py

@@ -1251 +1251 @@
 
 That situation is usually urealistic and undesirable, and when it
-happens, the simulator will issue a warning describing the problem. If
-you experience this problem, you can instruct the tool to guess
-current values using one of the following algorithms:
+happens the simulator will issue a warning describing the problem. If
+you experience this problem you can instruct the tool to guess current
+values using one of the following algorithms:
 
 * Del2a - Laplacian interpolation and linear extrapolation. We
@@ -1268 +1268 @@
 
 * Del4 - Same as Del2a but instead of the Laplace operator (also
-  called the del^2 operator) it uses the biharmoic operator (also
+  called the del^2 operator) it uses the biharmonic operator (also
   called the del^4 operator). May result in more accurate
   interpolations, at some cost in speed.
@@ -1285 +1285 @@
 estimates that are likely to be inaccurate and should only be used
 when you are willing to accept those inaccuracies because there is no
-other way to proceed."""),
+other way to proceed.
+
+Thanks to John D'Errico for providing the code that implements the
+mathematical algorithms described here (click `here
+<http://www.mathworks.com/matlabcentral/fileexchange/4551>`_ for more
+information)."""),
     arcGISDisplayName=_(u'Method for estimating missing currents values'))
 

MGET/Branches/Jason/PythonPackage/src/GeoEco/DataManagement/ArcGISRasters.py

@@ -874 +874 @@
         # Perform additional validation.
 
-        from GeoEco.Types import EnvelopeTypeMetadata
         left, bottom, right, top = EnvelopeTypeMetadata.ParseFromArcGISString(extent)
 
@@ -926 +925 @@
         # Perform additional validation.
 
-        from GeoEco.Types import EnvelopeTypeMetadata
         left, bottom, right, top = EnvelopeTypeMetadata.ParseFromArcGISString(extent)
 
@@ -1722 +1720 @@
 
     @classmethod
+    def ProjectToTemplate(cls, inputRaster, templateRaster, outputRaster, resamplingTechnique, interpolationMethod=None, maxHoleSize=None, mask=False, overwriteExisting=False):
+        try:
+            # Perform additional validation.
+
+            if interpolationMethod is not None and resamplingTechnique.lower() not in [u'bilinear', u'cubic']:
+                raise ValueError(_(u'Values cannot be interpolated for No Data cells if the NEAREST or MAJORITY resampling technique is used. Please select the BILINEAR or CUBIC resampling technique instead.'))
+
+            # Look up the coordinate system, extent, and cell size of
+            # the template raster.
+
+            gp = GeoprocessorManager.GetWrappedGeoprocessor()
+            d = gp.Describe(templateRaster)
+            coordinateSystem = gp.CreateSpatialReference(d.SpatialReference).split(';')[0]
+            extent = d.Extent
+            cellSize = d.MeanCellWidth
+            
+            # Below, we will set the Extent environment variable and
+            # then call ArcGIS's ProjectRaster_management to
+            # simultaneously project and clip the input raster to the
+            # coordinate system, extent, and cell size of the template
+            # raster. Unfortunately ProjectRaster_management will not
+            # necessarily create a raster that has exactly the desired
+            # extent. It may be one cell larger or smaller in any of
+            # the four directions. Different versions of ArcGIS seem
+            # to work differently in this respect.
+            #
+            # To deal with this annoyance, we will expand the extent
+            # in all four directions by 10 cells--guaranteeing that we
+            # have a raster that is larger than the desired
+            # extent--and then use ExtractByMask_sa to obtain a raster
+            # of the desired extent.
+            #
+            # We use 10 cells rather than 1 cell because of a
+            # different problem: if the caller requested that we
+            # interpolate values for No Data cells, the most accurate
+            # values can be obtained if each No Data region is
+            # completely surrounded by cells with data. In the event
+            # that the template raster extent bisects a No Data
+            # region, the 10 cell buffer increases the chance that we
+            # will interpolate that region using cells from both sides
+            # of the extent line: the side that is within the template
+            # extent and also the side that is outside the extent but
+            # within the buffer.
+
+            from GeoEco.DataManagement.Directories import TemporaryDirectory
+            tempDir = TemporaryDirectory()
+
+            [left, bottom, right, top] = EnvelopeTypeMetadata.ParseFromArcGISString(extent)
+            oldExtent = gp.Extent
+            gp.Extent = '%r %r %r %r' % (left - cellSize*10, bottom - cellSize*10, right + cellSize*10, top + cellSize*10)
+
+            try:
+                # If this is ArcGIS 9.3, set the SnapRaster
+                # environment variable to the template raster.
+                
+                if GeoprocessorManager.GetArcGISMajorVersion() > 9 or GeoprocessorManager.GetArcGISMajorVersion() == 9 and GeoprocessorManager.GetArcGISMinorVersion() >= 3:
+                    oldSnapRaster = gp.SnapRaster
+                    gp.SnapRaster = templateRaster
+
+                # Project the raster. If this is ArcGIS 9.3 or later,
+                # the SnapRaster and Extent environment variables will
+                # ensure the projected cells are snapped and clipped
+                # as requested by the caller. If it is 9.2, use the
+                # Registration_Point parameter instead of SnapRaster,
+                # which did not exist in 9.2 (there was something like
+                # it, but it did not work correctly).
+
+                try:
+                    projectedRaster = os.path.join(tempDir.Path, u'projected.img')
+                    if GeoprocessorManager.GetArcGISMajorVersion() > 9 or GeoprocessorManager.GetArcGISMajorVersion() == 9 and GeoprocessorManager.GetArcGISMinorVersion() >= 3:
+                        gp.ProjectRaster_management(inputRaster, projectedRaster, coordinateSystem, resamplingTechnique, cellSize)
+                    else:
+                        gp.ProjectRaster_management(inputRaster, projectedRaster, coordinateSystem, resamplingTechnique, cellSize, None, extent.rsplit(' ', 2)[0])
+
+                # If this is ArcGIS 9.3, reset the SnapRaster
+                # environment variable to what it was before.
+                
+                finally:
+                    if GeoprocessorManager.GetArcGISMajorVersion() > 9 or GeoprocessorManager.GetArcGISMajorVersion() == 9 and GeoprocessorManager.GetArcGISMinorVersion() >= 3:
+                        gp.SnapRaster = oldSnapRaster
+
+            # Reset the Extent ArcGIS environment variable to what it
+            # was before.
+            
+            finally:
+                gp.Extent = oldExtent
+
+            # If the caller requested that we interpolate values for
+            # No Data regions, do it now.
+
+            if interpolationMethod is not None:
+                from GeoEco.SpatialAnalysis.Interpolation import Interpolator
+                infilledRaster = os.path.join(tempDir.Path, u'infilled.img')
+                Interpolator.InpaintArcGISRaster(projectedRaster, infilledRaster, interpolationMethod, maxHoleSize)
+                projectedRaster = infilledRaster
+
+            # We are about to use ExtractByMask_sa to extract a raster
+            # of the desired extent from the projected raster in the
+            # temp directory. But ExtractByMask_sa also has a side
+            # effect that we might not want: it sets cells that are No
+            # Data in the mask to No Data in the output raster. If the
+            # caller did not request that to happen (using the
+            # template raster as the mask), create a constant raster
+            # that has the same coordinate system, extent, and cell
+            # size as the template raster. We'll use it instead.
+
+            if not mask:
+                oldOutputCoordinateSystem = gp.OutputCoordinateSystem
+                gp.OutputCoordinateSystem = coordinateSystem
+                try:
+                    maskRaster = os.path.join(tempDir.Path, u'constant.img')
+                    gp.CreateConstantRaster_sa(maskRaster, 0, 'INTEGER', cellSize, extent)
+                finally:
+                    gp.OutputCoordinateSystem = oldOutputCoordinateSystem
+
+                # For safety, verify that the output raster has the
+                # expected number of columns and rows.
+                
+                d2 = gp.Describe(maskRaster)
+                if d2.Width != d.Width or d2.Height != d.Height:
+                    raise RuntimeError(_(u'Programming error in this tool: the constant raster does not have the same number of rows or columns as the template raster. Please contact the author of this tool for assistance.'))
+
+            else:
+                maskRaster = templateRaster
+
+            # Extract the raster.
+
+            gp.ExtractByMask_sa(projectedRaster, maskRaster, outputRaster)
+
+        except:
+            Logger.LogExceptionAsError()
+            raise
+
+    @classmethod
     def ToPolygons(cls, inputRaster, outputFeatureClass, simplify=True, field=None, projectedCoordinateSystem=None, geographicTransformation=None, resamplingTechnique=None, projectedCellSize=None, registrationPoint=None, clippingDataset=None, clippingRectangle=None, mapAlgebraExpression=None, overwriteExisting=False):
         cls.__doc__.Obj.ValidateMethodInvocation()
@@ -2928 +3060 @@
 CopyArgumentMetadata(BinaryRaster.ToArcGISRaster, u'overwriteExisting', ArcGISRaster.ProjectClipAndOrExecuteMapAlgebra, u'overwriteExisting')
 ArcGISRaster.ProjectClipAndOrExecuteMapAlgebra.__doc__.Obj.GetArgumentByName(u'overwriteExisting').ArcGISCategory = None
+
+# Public method: ArcGISRaster.ProjectToTemplate
+
+AddMethodMetadata(ArcGISRaster.ProjectToTemplate,
+    shortDescription=_(u'Projects a raster to the coordinate system, cell size, and extent of a template raster.'),
+    isExposedToPythonCallers=True,
+    isExposedByCOM=True,
+    isExposedAsArcGISTool=True,
+    arcGISDisplayName=_(u'Project Raster to Template'),
+    arcGISToolCategory=_(u'Spatial Analysis\\Project Raster to Template'),
+    dependencies=[ArcGISDependency(9, 1)])
+
+CopyArgumentMetadata(ArcGISRaster.Copy, u'cls', ArcGISRaster.ProjectToTemplate, u'cls')
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'inputRaster',
+    typeMetadata=ArcGISRasterLayerTypeMetadata(mustExist=True),
+    description=_(
+u"""Raster to project to the template raster's coordinate system, cell
+size, and extent."""),
+    arcGISDisplayName=_(u'Raster to project'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'templateRaster',
+    typeMetadata=ArcGISRasterLayerTypeMetadata(mustExist=True, mustBeDifferentThanArguments=[u'inputRaster']),
+    description=_(
+u"""Raster defining the coordinate system, cell size, and extent of
+the output raster."""),
+    arcGISDisplayName=_(u'Template raster'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'outputRaster',
+    typeMetadata=ArcGISRasterTypeMetadata(mustBeDifferentThanArguments=[u'inputRaster', u'templateRaster'], deleteIfParameterIsTrue=u'overwriteExisting', createParentDirectories=True),
+    description=_(
+u"""Output raster.
+
+If this path refers to the file system, missing directories in the
+path will be created if they do not exist."""),
+    direction=u'Output',
+    arcGISDisplayName=_(u'Output raster'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'resamplingTechnique',
+    typeMetadata=UnicodeStringTypeMetadata(makeLowercase=True, allowedValues=[u'BILINEAR', u'CUBIC']),
+    description=_(
+u"""Resampling algorithm to be used to project the input raster to the
+template raster's coordinate system. One of:
+
+* NEAREST - `Nearest neighbor assignment
+  <http://en.wikipedia.org/wiki/Nearest-neighbor_interpolation>`_.
+
+* BILINEAR - `Bilinear interpolation
+  <http://en.wikipedia.org/wiki/Bilinear_interpolation>`_.
+
+* CUBIC - Cubic convolution, also known as `bicubic interpolation
+  <http://en.wikipedia.org/wiki/Bicubic_interpolation>`_.
+
+* MAJORITY - Majority resampling. This method requires ArcGIS 9.3 or
+  later.
+
+The NEAREST and MAJORITY algorithms should be used for categorical
+data, such as a land use classification. It is not recommended that
+NEAREST or MAJORITY be used for continuous data, such as elevation
+surfaces.
+
+The BILINEAR and CUBIC options are most appropriate for continuous
+data. Do not use BILINEAR or CUBIC with categorical data.
+
+The projection is accomplished with the ArcGIS Project Raster tool.
+Please see the documentation for that tool for more information."""),
+    arcGISDisplayName=_(u'Resampling technique'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'interpolationMethod',
+    typeMetadata=UnicodeStringTypeMetadata(makeLowercase=True, allowedValues=[u'Del2a', u'Del2b', u'Del2c', u'Del4', u'Spring'], canBeNone=True),
+    description=_(
+u"""Method to use to guess values for No Data cells.
+
+Use this option to "fill in" clusters of No Data cells with values
+obtained by interpolation and extrapolation. This option is
+appropriate for rasters representing continuous surfaces, e.g. images
+of sea surface temperature in which cloudy pixels contain No Data. It
+uses algorithms based on differential calculus that may provide more
+accurate guesses than traditional ArcGIS approaches, such as computing
+the focal mean of a 3x3 neighborhood.
+
+This option is not appropriate for rasters representing categorical
+data, such as land cover classifications. Therefore, in order to use
+this option, you must select BILINEAR or CUBIC for the Resampling
+Technique parameter.
+
+The available algorithms are:
+
+* Del2a - Laplacian interpolation and linear extrapolation.
+
+* Del2b - Same as Del2a but does not build as large a linear system of
+  equations. May be faster than Del2a at the cost of some accuracy.
+
+* Del2c - Same as Del2a but solves a direct linear system of equations
+  for the No Data values. Faster than both Del2a and Del2b but is the
+  least robust to noise on the boundaries of No Data cells and least
+  able to interpolate accurately for smooth surfaces.
+
+* Del4 - Same as Del2a but instead of the Laplace operator (also
+  called the del^2 operator) it uses the biharmonic operator (also
+  called the del^4 operator). May result in more accurate
+  interpolations, at some cost in speed.
+
+* Spring - Uses a spring metaphor. Assumes springs (with a nominal
+  length of zero) connect each cell with every neighbor (horizontally,
+  vertically and diagonally). Since each cell tries to be like its
+  neighbors, extrapolation is as a constant function where this is
+  consistent with the neighboring nodes.
+
+This option is applied after the input raster has been projected to
+the coordinate system and cell size of the template raster.
+
+Although this tool can fill No Data clusters of any size, you should
+apply common sense when using it. The larger the cluster, the less
+accurate the guessed values will be, especially for rasters that
+represent a noisy surface.
+
+Thanks to John D'Errico for providing the code that implements the
+mathematical algorithms described here (click `here
+<http://www.mathworks.com/matlabcentral/fileexchange/4551>`_ for more
+information)."""),
+    arcGISDisplayName=_(u'Method for interpolating No Data cells'),
+    arcGISCategory=_(u'Interpolation and masking options'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'maxHoleSize',
+    typeMetadata=IntegerTypeMetadata(mustBeGreaterThan=0, canBeNone=True),
+    description=_(
+u"""Maximum size, in cells of the template raster, that a region of
+4-connected No Data cells may be for it to be filled in.
+
+Use this option to prevent the filling of large No Data regions (e.g.
+large clouds in remote sensing images) when you are concerned that
+values cannot be accurately guessed for those regions. If this option
+is omitted, all regions will be filled, regardless of size."""),
+    arcGISDisplayName=_(u'Maximum size of No Data regions to interpolate'),
+    arcGISCategory=_(u'Interpolation and masking options'))
+
+AddArgumentMetadata(ArcGISRaster.ProjectToTemplate, u'mask',
+    typeMetadata=BooleanTypeMetadata(),
+    description=_(
+u"""If True, the template raster will be used to mask the input raster
+after it has been projected to the template raster's coordinate system
+and values have been interpolated for No Data cells (if your requested
+that). Cells of the template raster that are No Data will be set to No
+Data in the output raster, even if you requested that values be
+interpolated for No Data cells. This is appropriate in situations
+where the template defines the areas for which you want to retain
+data; for example, when you are analyzing the ocean and you have a
+mask in which ocean cells have data and land cells are set to No Data.
+
+If False, the template raster will only be used to define the
+coordinate system, cell size, and rectangular extent of the output
+raster, and no masking will be done."""),
+    arcGISDisplayName=_(u'Use template raster as mask'),
+    arcGISCategory=_(u'Interpolation and masking options'))
 
 # Public method: ArcGISRaster.ToPolygons
@@ -3382 +3669 @@
     findAndProcessMethodName=u'FindAndExtractByMask',
     findAndProcessMethodArcGISDisplayName=u'Find ArcGIS Rasters and Extract By Mask',
-    findAndProcessMethodShortDescription=_(u'Finds and rasters in an ArcGIS workspace and extracts the cells that correspond to the areas defined by a mask.'),
+    findAndProcessMethodShortDescription=_(u'Finds rasters in an ArcGIS workspace and extracts the cells that correspond to the areas defined by a mask.'),
     findMethod=ArcGISRaster.FindAndCreateTable,
     findOutputFieldParams=[u'rasterField'],
@@ -3446 +3733 @@
     findAndProcessMethodName=u'FindAndProjectClipAndOrExecuteMapAlgebra',
     findAndProcessMethodArcGISDisplayName=u'Find ArcGIS Rasters and Project, Clip, and/or Execute Map Algebra',
-    findAndProcessMethodShortDescription=_(u'Finds and rasters in an ArcGIS workspace and projects, clips, and/or perfoms map algebra on them. You must request at least one of these three operations. If you request multiple operations, the tool performs them in the order they are listed.'),
+    findAndProcessMethodShortDescription=_(u'Finds rasters in an ArcGIS workspace and projects, clips, and/or perfoms map algebra on them. You must request at least one of these three operations. If you request multiple operations, the tool performs them in the order they are listed.'),
     findMethod=ArcGISRaster.FindAndCreateTable,
     findOutputFieldParams=[u'rasterField'],

MGET/Branches/Jason/PythonPackage/src/GeoEco/SpatialAnalysis/Interpolation.py

@@ -907 +907 @@
 
 Use this tool to guess values for small clusters of No Data cells in
-rasters that represents continuous surfaces, e.g. images of sea
-surface temperature for which some pixels have been obscured by
-scattered clouds. This tool provides several advantages over
-traditional moving-window methods provided by ArcGIS, such as the
-Focal Statistics tool of the ArcGIS Spatial Analyst:
+rasters representing continuous surfaces, e.g. images of sea surface
+temperature in which cloudy pixels contain No Data. This tool provides
+several advantages over traditional moving-window methods provided by
+ArcGIS, such as the Focal Statistics tool of the ArcGIS Spatial
+Analyst:
 
 * It uses methods based on differential calculus that may provide more
-  accurate guesses than traditional approaches (e.g. computing the
-  focal mean of a 3x3 neighborhood).
-
-* It can fill No Data regions of any size.
-
-* It accurately handles rasters representing global images, where the
-  east and west edges are connected.
+  accurate guesses than traditional approaches, such as computing the
+  focal mean of a 3x3 neighborhood.
+
+* It accurately handles rasters with global longitudinal extent, for
+  which the east and west edges are connected.
+
+* It can fill No Data clusters of any size.
 
 Although this tool can fill No Data clusters of any size, you should
@@ -947 +947 @@
     typeMetadata=ArcGISRasterTypeMetadata(mustBeDifferentThanArguments=[u'inputRaster'], deleteIfParameterIsTrue=u'overwriteExisting', createParentDirectories=True),
     description=_(
-u"""Output raster.
-
-The output raster will be identical to the input raster where the
+u"""The output raster will be identical to the input raster where the
 input raster has data. For cells of the input raster that do not have
 data, the output raster will contain values interpolated according to
@@ -976 +974 @@
 
 * Del4 - Same as Del2a but instead of the Laplace operator (also
-  called the del^2 operator) it uses the biharmoic operator (also
+  called the del^2 operator) it uses the biharmonic operator (also
   called the del^4 operator). May result in more accurate
   interpolations, at some cost in speed.
@@ -1075 +1073 @@
     findAndProcessMethodName=u'FindAndInpaintArcGISRasters',
     findAndProcessMethodArcGISDisplayName=u'Find ArcGIS Rasters and Interpolate No Data Cells',
-    findAndProcessMethodShortDescription=_(u'Finds rasters in an ArcGIS workspace and values for the No Data cells.'),
+    findAndProcessMethodShortDescription=_(u'Finds rasters in an ArcGIS workspace and interpolates values for the No Data cells.'),
     findMethod=ArcGISRaster.FindAndCreateTable,
     findOutputFieldParams=[u'rasterField'],