[GIS] Workaround for Raster data input in ArcGIS Server Geoprocessing Service

arcgis-flex-apiarcgis-serverarcpygeoprocessingraster

UPDATE

I have accomplished the Python script, but only tested in desktop environment, i.e. ArcMap.

#Import ArcPy site package module
import arcpy
from arcpy import env
from arcpy import Extent
from arcpy import Raster
from arcpy.sa import ZonalStatistics
from arcpy.sa import CreateConstantRaster

# Check out ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set Environment Workspace
ws = env.workspace = "WORKSPACE_PATH\\Data.gdb"
# Set Environment Raster Cell Size
cellSize = env.cellSize = 10
# Set Environment Processing Extent
processingExtent = Extent(602560.330672585, 3827530.65305605, 746490.330672585, 3966810.65305605)
environProcessingExtent = env.extent = processingExtent
# Constant value for Constant Raster
constantValue = 0

# Get Parameters as Text
# Strings representing rasters in the workspace
rasters = arcpy.GetParameterAsText(0).split(";")
# The dataset that defines the zones 
inZoneData = arcpy.GetParameterAsText(1)
# zoneField variable will be determined on a conditional statement, but the logic has not been implemented yet.
# Field that holds the values that define each zone
zoneField = "TRACTCE10"
# Output location 
out = arcpy.GetParameterAsText(2)
outputLocation = ws + "\\" + out

# Create a constant value raster, where value = 0
addedRaster = CreateConstantRaster(constantValue, "FLOAT", cellSize, processingExtent)

# Process each raster string
for raster in rasters:
    addedRaster = addedRaster + Raster(ws + "\\" + raster)

# Execute Zonal Statistics
zoneStats = ZonalStatistics(inZoneData, zoneField, addedRaster, "MAXIMUM", "DATA")

# Save the output
zoneStats.save(outputLocation)

ORIGINAL QUESTION

A number of my questions on GIS.SE have been related to accomplishing a single, overarching problem that I have run into. In a nutshell, I want to enable a user of a Flex client application (using the ArcGIS API for Flex) to be able to select (via checkboxes) raster data to be used in a geoprocessing task that does the following: generates a Map Algebra expression to be used in a Raster Calculator operation and then performs Zonal Statistics (the user can also select the zone data attribute), then return a result map service with predefined symbology.

The two primary challenges to this task have been twofold:
1) Creating a dynamic map algebra expression that is dependent on the rasters selected by the user. Additionally, the documentation states that the Raster Calculator tool is not to be used outside of ArcMap and is not available in ArcPy. I have found that this similar task can be accomplished with the Spatial Analyst module and using the Addition Operator, but to accomplish what I have in mind would seem to require more processing time , i.e. 

rasOne + rasTwo = rasThree ... rasThree + rasFour = rasFive ... rasFive + ras6 = ras7 ... etc

2) And the primary and seemingly insurmountable hurtle: using raster data as the input. I have tried various approaches [1], [2], and [3].

The ESRI documentation states that raster data cannot be used as input to a geoprocessing service/task. However, I have found the company Azavea has accomplished this via their Decision Tree product [4], [5]. I also found this whitepaper that states the following about possible workarounds for using raster data as inputs:

Raster Workaround Options

What I am looking for in an answer is a concrete example of one of these workarounds for raster input implemented in Flex and a Python example for the dynamic map expression. Please let me know if the relevant Flex application code I have developed, the structure of our IT architecture, or an elaboration on the models we are trying to implement will be useful.

Best Answer

From a technical standpoint, I don't think this sounds too difficult.

One issue you will face is that Python is slow and Raster I/O with geoprocessing services is also slow. If the analysis is being done synchronously (and not emailed to the user an hour later), I would really recommend checking out server object extensions (C#/Java).

I could imagine the basic workflow as follows:

  1. Load all raster that users will be able to select into one standard map service with your predefined symbology and proivde a look up table (XML / Parse JSON from REST endpoint etc.) to the Flex client so it knows the map service index associated with each of your checkboxes.

  2. Create REST-based SOE which takes raster map service index, map algebra expression, and geometry (json) digitized by the user. One major advantage you will have is that you can load rasters into memory before hand which will really help with performance. Map algebra is certainly possible in the SOE, but I would avoid a ton of it. Really make an effort to pre-crunch as much as you can and simply run the zonal statistics operation on the canned rasters. I think your user will thank you.

I absolutely love Python, but my feeling is that geoprocessing services with rasters / arcpy are slow. I just think there is too much overhead in instantiating the gp object, and loading rasters into memory.

I wrote up a quick blog post to help get you going on the SOE if you're interested.

Hope you're doing well.

Related Question