Looping through equations in text file, Python script

[DiannProsser]

Hi all,

I am working on a script that will iterate through a number of geoprocessing tasks for each of 30+ animal species.
I am able to use os.listdir to list the inputs from different directories, then set up a for loop to have it use the species code (7 character code) in the file name to call the correct input files.
For my Single Output Map Algebra step, I would like to call the equation from a text file so that all 30+ equations are in a managable place, however I'm having trouble figuring out how to do this.

Attached is my current draft code; towards the bottom (see #Process 5) is the part where I'm getting stuck (see blue text). I'm attempting to open the text file, read the species code, and use the equation for the correct species for that iteration. As an example, here are the contents of the text file (just for 2 species and with truncated equations for testing)- not sure if I set this up right:

Example contents of the text file:
speciesCode equation
069swgo ( ild64p + ild41p + ild42p ) > (0)
074bhgo ( ild64p + ild41p ) > (0)

I'm using Arc10.0 SP3, Python 2.6, running script using PythonWin.
Finally, here is the output from the interactive window; the error message in red:

Printing shpRange
['r069swgo.shp', 'r074bhgo.shp']
printing valListShpBr
['v069swgo_br.shp', 'v074bhgo_br.shp']
Printing speciesCode: 069swgo
Printing currentValName: v069swgo_br
Printing shapfile: r069swgo.shp
Printing shapefileName: r069swgo
Running 50km buffer....
Clipping to validation layer....
Converting Vector Range to Raster (Polygon to Raster)...
Traceback (most recent call last):
File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 312, in RunScript
exec codeObject in __main__.__dict__
File "D:\GIS_Watf\Watf_Br.py", line 95, in <module>
print "Current Species Code from Text file is:" + " " + currentSppCode #just so I can see which one it's calling
NameError: name 'currentSppCode' is not defined

 
#DJP Updated 01 February 2012
#Revised for ArcGIS 10 on GIS desktop
#tips:
#   slicing: general form is l[start:end:step]
#   to indent, use 4 spaces not a tab since tabs can be interpreted differently depending on the platform or tool being used
#   if, while, and for statements

# Import system modules
import sys, string, os, arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create()

# Set the necessary product code
gp.SetProduct("ArcInfo")

# Check out any necessary licenses
gp.CheckOutExtension("spatial")

#Python will overwrite exsisting outputfiles
gp.overwriteoutput = 1

# Load required toolboxes...(## commands below are for ArcGIS 9.*)
gp.AddToolbox("C:/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
gp.AddToolbox("C:/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Conversion Tools.tbx")
gp.AddToolbox("C:/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Analysis Tools.tbx")
#gp.AddToolbox("C:/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
#gp.AddToolbox("C:/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
#gp.AddToolbox("C:/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
#gp.AddToolbox("C:/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

gp.workspace = "D:/GIS_Watf" #(**change back to C: for laptop)
outWorkspace = "D:/GIS_Watf/Output/"
inRange = "D:/GIS_Watf/Data/Ranges/" #path to input range maps

RangeList = os.listdir(inRange) #list what's in the Ranges folder

shpRange = [] #create empty variable 
[shpRange.append(i) for i in RangeList if i[-4:] == '.shp'] #for each file in the RangeList list, append if it ends in ".shp"
print "Printing shpRange"
print shpRange

inValBr = "D:/GIS_Watf/Data/Valida/Breeding/" #create variable named inValBr that shows path to Breeding Validation data
inValWi = "D:/GIS_Watf/Data/Valida/Winter/" #create variable named inValBr that shows path to Breeding Validation data
valListFull = os.listdir(inValBr) #list all files in Valida/Breeding folder
valListShpBr = []
[valListShpBr.append(i) for i in valListFull if i[-4:] == '.shp'] #append to list only if they are shapefiles
print "printing valListShpBr"
print valListShpBr

for shapefile in shpRange: #starting a loop with new variable "shapefile" that loops thru shpRange
    shapefileName = shapefile[:-4] #new variable shapefileName that holds the shapefile name minus the .shp 
    speciesCode = shapefileName[1:8] #new variable speciesCode that holds the 2nd thru 4th character from the file name (this is the species number)
    print "Printing speciesCode:" + " " + speciesCode #print message and the species number
   
    for valShapefile in valListShpBr: #starts new loop with variable valShapefile and loops thru the valListShp list
        if speciesCode in valShapefile: #chooses validation file having the speciesCode (from Range file) 
            currentValFile = valShapefile
    currentValName = currentValFile[:-4] #currentValName is the validation shapefile (minus the ".shp")
    print "Printing currentValName:" + " " + currentValName
    print "Printing shapfile:" + " " + shapefile
    print "Printing shapefileName:" + " " + shapefileName

    # Process 1: BHGO Summer Range (Select)...
    selectOutput = outWorkspace + shapefileName + "_br.shp" #** change this to "_br.shp"
    selectInput = inRange + shapefile
    gp.Select_analysis(selectInput, selectOutput, "season = 'summer'")
    
    # Process 2: (Buffer) 50 km...
    print "Running 50km buffer...."
    bufferOutput = outWorkspace + shapefileName + "_buff50.shp"
    gp.Buffer_analysis(selectOutput, bufferOutput, "50 Kilometers", "FULL", "FLAT", "NONE", "")    
   
    # Process 3: Clip Winter Validation Points by Range...
    print "Clipping to validation layer...."
    valClipOutput = outWorkspace + currentValName + "_clip.shp"
    valInput = inValBr + currentValFile
    gp.Clip_analysis(valInput, bufferOutput, valClipOutput, "")

    # Process 4: Convert Vector Range to Raster  (Polygon to Raster)...'
    print "Converting Vector Range to Raster (Polygon to Raster)..."
    convertRangeOut = outWorkspace + shapefileName + "_br"
    gp.PolygonToRaster_conversion(bufferOutput, "season", convertRangeOut, "CELL_CENTER", "", "1000")

    # Process 5: Create Breeding Distribution Model for each species (Map Algebra)...
    # This process should find the correct species equation in a text file and perform the SOMA on it    SppEquaFile = "D:/GIS_Watf/Command/SpeciesEquations.txt"
    f = open(SppEquaFile)
    f.read()
    for speciesCode1 in f:
        if speciesCode1 in f:
            currentSppCode = speciesCode1 #Here I'm trying to get it to read the species code, move to the correct one, then use the equation that follows it; I realize I probably need another step here...    print "Current Species Code from Text file is:" + " " + currentSppCode #just so I can see which one it's calling

    v1_if_True = "1"
    v0_if_False = "0"
    SOMAoutput = outWorkspace + "a" + currentSppCode + "_br" #example output name is a069swgo_br
    gp.SingleOutputMapAlgebra_sa(ImNotSure, SOMAoutput) #here see 6 lines above where I'm getting stuck; ImNotSure will change to the correct input variable...

Thank you in advance for any help you can provide. As you might guess, I'm new to Python.

[ChrisSnyder]

Why not store the equations in the python script, which is also a text file... Then you only have one file to worry about.... But...

Maybe you want to use the .readlines() method instead?

If the text file looks like this:
my
cat
is
bad

x = f.read() will produce:

'my\ncat\nis\nbad'

where x = f.readlines() will producee a list of the lines like this:

['my\n', 'cat\n', 'is\n', 'bad']


Maybe even better, you can read the text file into a Dictionary object. For example if the text file looks like this:

cat snipps+snails
dog sugar+spice

f = open(txtFile)
x = f.readlines()
f.close()
eqDict = {}
for i in x:
   species = i.split(" ")[0]
   equation = i.split(" ")[-1].replace("\n", "")
   eqDict[species] = equation

Then you have a in-memory look up table you can access like this:

>>> print eqDict["cat"]
'snipps+snails'

[DiannProsser]

Thank you - I used your Dictionary object example and included a for/if loop, and it appears to be finding the correct equation.
I am having a problem with the SOMA process - it won't recognize my input variables (see ild64p in error message below as an example). I put them in their own directory, but added that within my syntax. Can you identify where I'm going wrong?

 
 # Process 5: Create Breeding Distribution Model for each species (Map Algebra)...
    # This process should find the correct species equation in a text file and perform the SOMA on it

    SppEquaFile = "D:/GIS_Watf/Command/SpeciesEquations.txt"
    f = open(SppEquaFile)
    x = f.readlines()
    print x
    f.close()
    eqDict = {}
    for i in x:
        species = i.split(" ") [0]
        print "species is:"
        print species
        equation = i.split("x") [-1].replace("\n", "") #the split is denoted by "x"; also .replace replaces "\n" with nothing so we are left with the equation
        print "equation is:"
        print equation
        eqDict[species] = equation
        if species == speciesCode:
            SOMAequa = equation
        print "SOMAequa is:"
        print SOMAequa
        
    inPredictors = "D:/GIS_Watf/Data/LULC/" #identifying directory for the SOMA inputs    
    v1_if_True = "1"
    v0_if_False = "0"
    SOMAoutput = outWorkspace + "a" + speciesCode + "_br" #example output name is a069swgo_br
    gp.SingleOutputMapAlgebra_sa (SOMAequa, SOMAoutput,inPredictors)  #using the optional syntax for in_data parameter (inPredictors)

the error message:

Traceback (most recent call last):
File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 312, in RunScript
exec codeObject in __main__.__dict__
File "D:\GIS_Watf\Watf_Br2_forum.py", line 112, in <module>
gp.SingleOutputMapAlgebra_sa (SOMAequa, SOMAoutput,inPredictors) #using the optional syntax for in_data parameter (inPredictors)
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000865: Map Algebra expression: ild64p does not exist.
Failed to execute (SingleOutputMapAlgebra).

[ChrisSnyder]

Is ild64p the name of the actual raster datsaset (for example r"C:\temp\ild64p") or just a variable name?

If its an actual raster, are you setting the gp.workspace (to where the raster exists) so that ArcGIS knows where to go looking for the raster?

If its a variable name, hmm.... Basically you have to feed a text string into the SOMA tool.

inGrid = r"C:\temp\input"
outGrid =  r"C:\temp\output"
somaExp = "con(" + inGrid + " == 1, 1, setnull(1))"
gp.SingleOutputMapAlgebra_sa (somaExp, outGrid) #BTW: You are using a third parameter to the SOMA tool, which you don't need in the Python script (only in ModelBuilder)!

If the equations you are reading in contain variable names (but as text strings), well that's not doing to work... Since:

"con( + inGrid +  == 1, 1, setnull(1))"

is not the same as:

"con(" + inGrid + " == 1, 1, setnull(1))"

Right?

You could do something like this:

print "con( + inGrid +  == 1, 1, setnull(1))".replace("inGrid", inGrid)

and substitute the inGrid string for the inGrid variable (which points to a different string value as defined in your script).

[ChrisSnyder]

Whoops, I'm sorry that would be:

somaExp = "con(inGrid == 1, 1, setnull(1))".replace("inGrid", inGrid)

[DiannProsser]

Ahhh - I see the problem.
I forgot a space before the in_data parameter!
gp.SingleOutputMapAlgebra_sa (SOMAequa, SOMAoutput,inPredictors) <-- missed the space between the comma and in_data

Thank you so VERY much for your help today!!



(ild64p is one of the ten or so rasters that are included in the species equations)

[DiannProsser]

Oh no - I thought I fixed the problem but I was wrong.

So, ild64p is one of the ten or so rasters that are included in the species equations.
When I put the rasters in the gp.workspace, it runs fine.

But I wanted to have these rasters in a separate "Data" folder, not the workspace, hence I am trying to use the optional in_data parameter:
(Syntax): SingleOutputMapAlgebra_sa (expression_string, out_raster, in_data)

I identify the directory here:
inPredictors = "D:/GIS_Watf/Data/LULC/"

and include it in the syntax here:
gp.SingleOutputMapAlgebra_sa (SOMAequa, SOMAoutput, inPredictors)

Can you see where I'm going wrong?

Thank you...

File "D:\GIS_Watf\Watf_Br.py", line 112, in <module>
gp.SingleOutputMapAlgebra_sa (SOMAequa, SOMAoutput, inPredictors) #using the optional syntax for in_data parameter (inPredictors)
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000865: Map Algebra expression: ild64p does not exist.
Failed to execute (SingleOutputMapAlgebra).

[ChrisSnyder]

Unfortunatly the parameter you are supplying 'inPredictors' to is not used in the way you want... It is strictly a ModelBuilder parameter that is ignored when executed in Python code. It has nothing to do with utilizing rasters from differnt workspaces.

If your input rasters are all in differnt workspaces, you would have to do something like this:

r1 = r"C:\temp1\test1"
r2 = r"C:\temp39\test14"
r3 = r"C:\temp56\test177"
r4 = r"C:\temp563\test765"
eq = 'c1 + c2 + c3'
goodToGoEq = eq.replace("c1", r1).replace("c2", r2).replace("c3", r3).replace("c4", r4)
#Note that the text string "c4" doesn't exist in the equation, and therefore the .replace() for "c4" is ignored...

>>> print eq
'C:\\temp1\\test1 + C:\\temp39\\test14 + C:\\temp56\\test177'

Like I said earlier, it would be way easier to just store the equations in the Python script file rather than read them from a text file... The reason is that you want to interpret the equations WITH the variables (that are assigned in the script), and not just as the text strings they exits as in the text file. Howevre, if all the input raster existed in the same workspace, you could just use the raw strings as input to the SOMA tool.

There is probably a better way to do this sort of string/variable substitution thing, but...

[curtvprice]

There is probably a better way to do this sort of string/variable substitution thing, but...

r1 = r"C:\temp1\test1"
r2 = r"C:\temp39\test14"
r3 = r"C:\temp56\test177"
r4 = r"C:\temp563\test765"
eq = "{0} + {1} + {2} + {3}"
eq1 = eq.format(r1,r2,r3,r4)
print eq1

The code snippet above, if pasted in a python window, prints this:

C:\temp1\test1 + C:\temp39\test14 + C:\temp56\test177 + C:\temp563\test765

If those input paths are all the same, (that is, your "Data" workspace"), it seems to me setting the gp.workspace to that path would be less messy, then just specify an output path to write the results:

import os 
SOMAoutput = os.path.join(outWorkspace, "a" + speciesCode + "_br") 
#example output name is a069swgo_br

One more thing - although ModelBuilder exports all those AddToolbox calls, if they are standard ArcToolbox tools they are unnecessary -- the gp knows where to find them.

[DiannProsser]

The rasters are all in one directory: "D:/GIS_Watf/Data/LULC/"

Why is SOMA the only process that needs inputs to be in the current workspace?
For example, I have no trouble setting the input directory for my clip and buffer processes, as you can see from my code.
I realize this is not just a problem with Python; it was the same situation when I ran the process in Model Builder.
Is there no possibility with the SOMA process to specify a single separate input directory?

The reason I want inputs for the SOMA process in their own directory is because I have multiple processes and they each pull inputs from their own directories.
For example, the clip process clips features from the Range directory; the buffer process buffers the outputs from the clip process, etc.
Also, the reason why I want the equations in a separate text file is because they will be edited multiple times by experts in the field; best to do it there rather than have to fix the code multiple times.

Thank you both for your input throughout this thread (also for the tip about the AddToolbox calls, Curtis).