Trouble migrating SOMA expression into ArcGIS 10 using arcpy?

[JamesSnider]

Hi,

In 9.3, I had successfully applied the following nested Con statement:

    while row:
        value = row.getvalue("ID")
        if value == ID99:
            print "Value99..."
            value99 = row.getvalue("VALUE")
            print value99
        elif value == ID95:
            print "Value95..."
            value95 = row.getvalue("VALUE")
            print value95
        elif value == ID90:
            print "Value90..."
            value90 = row.getvalue("VALUE")
            print value90
        elif value == ID85:
            print "Value85..."
            value85 = row.getvalue("VALUE")
            print value85
        elif value == ID80:
            print "Value80..."
            value80 = row.getvalue("VALUE")
            print value80  
        row = rows.Next()

          
    gp.SingleOutputMapAlgebra_sa("CON ( " + str(raster) + " >= " + str(value99) + ", 99 , CON ( " + str(raster) + " >= " + str(value95) + ", 95, CON ( " + str(raster) + " >= " + str(value90) + ", 90, CON ( " + str(raster) + " >= " + str(value85) + ", 85, CON ( " + str(raster) + " >= " + str(value80) + " , 80, 0 ) ) ) ) )", Output)
    print "Output percentile raster created"

However, I've had some issues implementing the same processes using arcpy.
I've tried using both the Con and Reclassify approaches, but both raise errors.

Using the Con tool, the expression I've used is:

for row in rows:
        value99 = row.getValue("VAL_99")
        value95 = row.getValue("VAL_95")
        value90 = row.getValue("VAL_90")
        value85 = row.getValue("VAL_85")
        value80 = row.getValue("VAL_80")

try:
        outRas = Con(raster >= value99, 99, Con(raster >= value95, 95, Con(raster >= value90, 90, Con(raster >= value85, 85, Con(raster >= value80, 80, 0)))))
        print "Con statement completed"
except Exception as e:
        print e.message

    try:
        outRas.save(r"F:\outRas")
        print "output raster saved"
        
    except Exception as e:
        print e.message

Which raises the following error:

ERROR 010240: Could not save raster dataset to F:\outRas with output format GRID.

Similarly, when using the Reclassify approach:

try:
        in_raster = CAVM_het_output
        reclass_field = "VALUE"
        
        valueMAX = arcpy.GetRasterProperties_management(CAVM_het_output, "MAXIMUM")
        print "max = ", valueMAX
        
        valueMIN = arcpy.GetRasterProperties_management(CAVM_het_output, "MINIMUM")
        print "min = ", valueMIN
        
        remap = ([value99, valueMAX, 99], [value95, value99, 95], [value90, value95, 90], [value85, value90, 85],[value80, value85, 80], [valueMIN, value80, 0])
        print "remap = ", remap        

        outRas = Reclassify(in_raster, reclass_field, remap, "NODATA")
        print "outRas created"
        
        print "Reclassify statement completed"
        
    except Exception as e:
        print e.message

    try:
        outRas.save(r"F:/outRas")
        print "output raster saved"
        

This approach raises the following error:

Object: Error in executing tool

In both cases, my thinking is the trouble is attributed to having variables in the expressions (e.g. value99 or valueMAX), but perhaps I'm wrong.

E.g.

remap =  ([565.0, <Result '585'>, 99], [546.0, 565.0, 95], [525.0, 546.0, 90], [504.0, 525.0, 85], [483.0, 504.0, 80], [<Result '0'>, 483.0, 0])

Does anyone have any ideas here? I'm guessing this is an easy fix, but I'm hitting the wall...

Many thanks,
James

[JamesSnider]

The Reclassify approach seems to be most promising; no success yet, but here's a slight change to that code:

try:
        in_raster = CAVM_het_output
        reclass_field = "VALUE"
        
        valueMAX = arcpy.GetRasterProperties_management(CAVM_het_output, "MAXIMUM")
        print "max = ", valueMAX
        
        valueMIN = arcpy.GetRasterProperties_management(CAVM_het_output, "MINIMUM")
        print "min = ", valueMIN
        
        remap = RemapRange([[value99, 585, 99], [value95, value99, 95], [value90, value95, 90], [value85, value90, 85],[value80, value85, 80], [0, value80, 0]])
        print "remap = ", remap        

        outRas = Reclassify(in_raster, reclass_field, remap, "NODATA")
        print "outRas created"
        
        print "Reclassify statement completed"
        
    except Exception as e:
        print e.message

Which produces the remap statement:

remap =  565.0 585 99; 546.0 565.0 95; 525.0 546.0 90; 504.0 525.0 85; 483.0 504.0 80; 0 483.0 0

Problem is, now Python crashes when it attempts the Reclassify processing step!

[ChrisSnyder]

I feel your pain... I pulled my hair for a few hours and then gave up.... Try converting something that uses .ITEM notation (i.e. using fields in the .vat or RAT or whatever they call it now). Never figured out if/how you can do it... Never found anything in the help... Anyone know how to convert something like this?

somaExp = "con(" + grid + ".FIELDNAME == 1, 1, setnull(1))"

FYI: You can still use the good ole' SOMA and MOMA tools in v10. You just need to also instantiate:

import arcgisscripting, arcpy
gp = arcgisscripting.create(9.3)
arcpy.this(blah, blah)
gp.SingleOutputMapAlgebra(blah, blah)

[JamesSnider1]

Thanks, Chris, that is very helpful to know. I'll try invoking the 9.3 processor to use some of my old SOMA expressions using arcpy. As you noted, there is little help documentation available on how exactly to make these migrations, so any tidbits like the above are greatly appreciated.

In terms of this specific problem, I did find a workaround using the reclassify approach, for those who are interested:

import sys, string, os, arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True
Study_Unit_List = range(400, 420)
print "Study_Unit_List = ", Study_Unit_List

for x in Study_Unit_List:
    print "Study unit number = ", x
    
    #Specify input files:
    raster = r"C:\Inputs\ras" + str(x) 
    print "input raster = ", raster

    Stats_Table = r"C:\Inputs\ras" + str(x) + "_stats.dbf"
    print "Stats Table = ", Stats_Table
    
    # Get distribution percentile values from Statistics Table, dbf:

    rows = arcpy.SearchCursor(Stats_Table)
    for row in rows:
        value99 = row.getValue("VAL_99")
        value95 = row.getValue("VAL_95")
        value90 = row.getValue("VAL_90")
        value85 = row.getValue("VAL_85")
        value80 = row.getValue("VAL_80")

    print "threshold values calculated"

    #Reclassify input raster based on distribution percentiles:

    try:
        valueMAX = arcpy.GetRasterProperties_management(raster, "MAXIMUM")
        print "max = ", valueMAX
        
        valueMIN = arcpy.GetRasterProperties_management(raster, "MINIMUM")
        print "min = ", valueMIN
        
        myremap = RemapRange([[value99, 585, 99], [value95, value99, 95], [value90, value95, 90], [value85, value90, 85], [value80, value85, 80], [0, value80, 0]])      

        outRas = Reclassify(raster, "VALUE", myremap)
        
    except Exception as e:
        print e.message

    # Save output raster:
    try:
        outRas.save(raster + "_pr")
        print "output raster saved"
            
    except Exception as e:
        print e.message