On ArcGIS 10, Run Coral Reef Connectivity Simulation hangs on the "Writing density rasters" step

[jjr8]

Here is some example output:

Executing: CoralReefConnectivityRunSimulation C:\MGET_ArcGISToolbox_Tests_Output\Coral_Reef_Larval_Connectivity\SimDir1 C:\MGET_ArcGISToolbox_Tests_Output\Coral_Reef_Larval_Connectivity\OutputDir1 9/30/2009 30 2.4 10 10000 0.1 25 430;431;432;433;434;435;206;210;403;303;312;290;231;313 #
Start Time: Wed Sep 01 09:45:53 2010
Running script CoralReefConnectivityRunSimulation...
This tool may run slowly due to a geoprocessing performance degradation that ESRI introduced in ArcGIS 9.3. We are following up with ESRI about this problem and will try to get them to fix it. Please contact us if you have any questions. In the mean time, you can avoid this problem by running this tool on a computer with ArcGIS 9.1 or 9.2.
Reading coral reef data...
Query complete: 0:00:02 elapsed, 14 rows retrieved, 0:00:00.168142 per row.
Reading 10 ocean currents images...
Finished reading: 0:00:46 elapsed, 10 images read, 0:00:04.677800 per image.
The maximum Courant number is 1.#QNAN0, which is greater than or equal to 1. The simulation will be numerically unstable and the results will probably be inaccurate. To obtain a stable simulation, you must reduce the time step to 1.#QNAN or less, so that the maximum Courant number is less than 1.
Warning: The ocean currents images are missing data for some cells flagged as
containing reefs. This will affect the accuracy of the simulation. The
simulator will assume the ocean currents have a velocity of zero in these
cells. Larvae can only exit these cells via diffusion. These cells may retain
larvae in a manner that is not realistic. The IDs of the affected reefs are:
432, 433, 434
Simulating larval dispersal for 14 reefs...
Still simulating: 01:03 elapsed, 8 reefs processed, 07.889 seconds per reef, 6 reefs remaining, estimated completion time: 09:49:57
Simulation complete: 01:45 elapsed, 14 reefs processed, 07.530 seconds per reef.
Writing 26 edges to the edge list in the output geodatabase...
Finished inserting rows: 0:00:02 elapsed, 26 rows inserted, 0:00:00.104769 per row.
Created directory C:\MGET_ArcGISToolbox_Tests_Output\Coral_Reef_Larval_Connectivity\OutputDir1\DensityRasters.
Writing 31 density rasters to the output directory...

The tool appears to hang after the last message above. I let it run for 30 minutes with no further output.

The hang appears to occur in the ArcGIS SingleOutputMapAlgebra_sa tool. When I ran the simulation from a command-line geoprocessing script and pressed Ctrl-C when it hung, I got a traceback that suggested the hang was there:

C:\Program Files>Traceback (most recent call last):
  File "C:\PROGRA~1\GeoEco\ARCGIS~1\Scripts\CoralReefConnectivityRunSimulation.py", line 5, in <module>
    GeoEco.ArcGIS.ExecuteMethodFromCommandLineAsArcGISTool('GeoEco.Connectivity.CoralReefConnectivity', 'CoralReefConnectivity', 'RunSimulation')
  File "C:\Python26\ArcGIS10.0\lib\site-packages\GeoEco\ArcGIS.py", line 956, in ExecuteMethodFromCommandLineAsArcGISTool
    exec sourceCode in globals(), locals()
  File "<string>", line 1, in <module>
  File "C:\Python26\ArcGIS10.0\lib\site-packages\GeoEco\Connectivity\CoralReefConnectivity.py", line 643, in RunSimulation
    overwriteExisting=overwriteExisting)
  File "C:\Python26\ArcGIS10.0\lib\site-packages\GeoEco\DataManagement\ArcGISRasters.py", line 1419, in FromNumpyArray
    ArcInfoASCIIGrid.ToArcGISRaster(asciiFile, raster, isInteger=numpyArray.dtype.kind=='i', coordinateSystem=coordinateSystem, projectedCoordinateSystem=projectedCoordinateSystem, geographicTransformation=geographicTransformation, resamplingTechnique=resamplingTechnique, projectedCellSize=projectedCellSize, registrationPoint=registrationPoint, clippingRectangle=clippingRectangle, mapAlgebraExpression=mapAlgebraExpression, buildPyramids=buildPyramids, overwriteExisting=overwriteExisting)
  File "C:\Python26\ArcGIS10.0\lib\site-packages\GeoEco\DataManagement\ArcInfoASCIIGrids.py", line 145, in ToArcGISRaster
    gp.SingleOutputMapAlgebra_sa(mapAlgebraExpressionToExecute, outputRasterTemp2)
  File "<string>", line 2, in SingleOutputMapAlgebra_sa
  File "C:\Python26\ArcGIS10.0\lib\site-packages\GeoEco\ArcGIS.py", line 1357, in _InvokeMethod
    raise ArcGISError(_(u'The ArcGIS %(tool)s geoprocessing tool failed when given the parameters %(params)s and reported %(error)s: %(msg)s') % {u'tool': method.__name__, u'params': inspect.formatargvalues(_args, _varargs, _varkw, _locals), u'error': e.__class__.__name__, u'msg': unicode(e)})
GeoEco.ArcGIS.ArcGISError: The ArcGIS SingleOutputMapAlgebra_sa geoprocessing tool failed when given the parameters (*args=('setnull( C:\\MGET_ArcGISToolbox_Tests_Output\\Coral_Reef_Larval_Connectivity\\SimDir1\\ReefData\\water_mask == 0 OR  [C:\\Temp\\GeoEcoTemp_Jason\\tmpl_m212\\output]  < 10,  [C:\\Temp\\GeoEcoTemp_Jason\\tmpl_m212\\output] )', 'C:\\Temp\\GeoEcoTemp_Jason\\tmpl_m212\\output4')) and reported ExecuteError: ERROR 000582: Error occurred during execution.

When I commented out the SingleOutputMapAlgebra_sa call in the code, the tool did not hang.

This looks like a bug in ArcGIS 10. The simplest workaround is to eliminate the map algebra call and do the map algebra in memory using basic numpy operations. (The map algebra is just applying the water mask and setting pixels below a threshold to nodata.)

If I have time, I will try to create a small repro scenario for ESRI, and open a bug with them.

[jjr8]

Fixed in r558, released in MGET 0.8a12.