Map Algebra in v10

[ChrisSnyder]

Struggling with "best" Map Algebra syntax to use in scripts via PythonWin (NOT the Python window in ArcGIS, which seems to allow for different syntax - which is confussingm, but off topic)...

Anyway... Classic "make streams from flow accumulation" scenario...

In v9.3 I could do this:

flowAccGrd = "flow_acc" flowDirGrd = "flow_dir" somaExp = "streamlink(con(" + flowAccGrd + " > 500, 1, setnull(1)), " + flowDirGrd + ")" streamLinkGrd = "stream_link" gp.SingleOutputMapAlgeba_sa(somaExp, streamLinkGrd)

Solid...

Now as I try to rewrite some stuff in the new v10 syntax, the best I can com up with is this:

flowAccGrd = "flow_acc" flowDirGrd = "flow_dir" streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccGrd, 1, arcpy.sa.SetNull(1,1), "VALUE >= 500")), flowDirGrd) streamLinkGrd = "stream_link" streamLinkTmp.save(streamLinkGrd)

Is there a better/more compact way to write this expresion ini v10 sytax (again, writting for a "stand-alone script")?

Is that imbedded arcpy.sa.SetNull(1,1) neccessary - or is there a more elegant way to use the IsNull() and SetNull() expressions inside of a con() statement... Seemed easier in v9.3.

[PhilMorefield · Accepted solution]

Off the top of my head, cast your inputs as Raster objects to make things a little more compact. To wit:

flowAccGrd = arcpy.Raster("C:\\example\\flowaccgrd") ... streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccGrd > 500, 1), flowDirGrd) ...

Seems like that should do it. I don't think you even need SetNull, since the output of the Con statement will necessarily be constrained to values > 500.

[ChrisSnyder]

Thanks Phil... I was missing the whole "casting as a raster" thing... and the implied setnull - didn't know about that - very usefull!

Took me a while to realize that the temporary raster objects (for example "slopePctTmp" in the expresion: slopePctTmp = arcpy.sa.Slope(fillGrd, "PERCENT_RISE", "") can and should be used even after you ".save()" them, and low and behold, they now reference the updated .save() path! One you you use a variable to reference the .save() output path, that variable is no longer a "raster object", so just keep using the raster reference (slopePctTemp).

The big lightbulb though - and a **MAJOR ARCPY CRASH BUG ISSUE***:

What I had been doing: Seems that you CAN use the the output path variable (as opposed to the raster object variable) in SINGLE tool expreessions... for example:

#Process: Run a Fill
fillTmp = arcpy.sa.Fill(clipDem, "")
fillGrd = areaIdFolderPath + "\\fill"
fillTmp.save(fillGrd)
#Process: FlowDir
flowDirTmp = arcpy.sa.FlowDirection(fillGrd, "", "") #NOTE: I am using the fillGrd variable not the the fillTmp variable: BOTH work for "single tool expressions" 
flowDirGrd = areaIdFolderPath + "\\flowdir"
flowDirTmp.save(flowDirGrd)

But if you try to use the output path variables (and not the raster object variables) in MULTI TOOL EXPRESSIONS Python will CRASH BOOM BLOWUP... for example:

This is OK:

#Process: FlowAcc
flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage()
flowAccGrd = areaIdFolderPath + "\\flowacc"
flowAccTmp.save(flowAccGrd)

#Process: Create a streamlink
acreThreshold = 1
acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2)
flowAccGrd = arcpy.Raster(flowAccGrd)
streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccTmp > acreThresholdInPixels, 1), flowDirGrd)
streamLinkGrd = areaIdFolderPath + "\\strmlnk"
streamLinkTmp.save(streamLinkGrd)

This VERY bad and cause Python "Application Error" explosion:

#Process: FlowAcc
flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage()
flowAccGrd = areaIdFolderPath + "\\flowacc"
flowAccTmp.save(flowAccGrd)

#Process: Create a streamlink
acreThreshold = 1
acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2)
flowAccGrd = arcpy.Raster(flowAccGrd)
streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccGrd > acreThresholdInPixels, 1), flowDirGrd)
streamLinkGrd = areaIdFolderPath + "\\strmlnk"
streamLinkTmp.save(streamLinkGrd)

So the best solution is to ***ALWAYS*** use the raster objects andnot the path variables:

#Process: FlowAcc
flowAccTmp = arcpy.sa.FlowAccumulation(flowDirGrd, "", "INTEGER"); showGpMessage()
flowAccGrd = areaIdFolderPath + "\\flowacc"
flowAccTmp.save(flowAccGrd)

#Process: Create a streamlink
acreThreshold = 1
acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2)
flowAccGrd = arcpy.Raster(flowAccGrd)
streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccTmp > acreThresholdInPixels, 1), flowDirTmp)
streamLinkGrd = areaIdFolderPath + "\\strmlnk"
streamLinkTmp.save(streamLinkGrd)

So, moral of the story... Don't use the .save(pathvariable) path variables in raster expresions!!!

[ChrisSnyder]

New strugle:

I had written some code years ago that would create a 1 cell outline of the DATA portion of a raster (see attached image):

somaExp = "con(" + flowDirGrd + " > 0 & isnull(Focalmax(" + flowDirGrd + ", 'RECTANGLE', 3, 3, 'NODATA')), 1, setnull(" + flowDirGrd + "))"

....where 'flowDirGrd' is just some raster (doesn't neccesarily have to be a flow direction raster).

So basically just the outer DATA cells get coded as a 1 and everything else is NoData.

I came up with this for v10 sytax:

demRingTmp = arcpy.sa.Con(flowDirTmp > 0 & arcpy.sa.IsNull(arcpy.sa.FocalStatistics(flowDirTmp, arcpy.sa.NbrRectangle(3, 3, "CELL"), "MAXIMUM", "NODATA")), 1)

...which runs but produces the wrong results for some reason.

Anyone wanna take a stab at writting this ugly thing in v10 syntax (and get the correct result)?


[ATTACH=CONFIG]12500[/ATTACH]

[PhilMorefield]

Add two sets of parentheses, one around each of your conditions. So:

# your way; gave me a raster full of 1s
demRingTmp = arcpy.sa.Con(flowDirTmp > 0 & arcpy.sa.IsNull(arcpy.sa.FocalStatistics(flowDirTmp, arcpy.sa.NbrRectangle(3, 3, "CELL"), "MAXIMUM", "NODATA")), 1)

# my way; gave me a solid 'ring' around the data
from arcpy import sa
demRingTmp = sa.Con((flowDirTemp > 0) & (sa.IsNull(sa.FocalStatistics(flowDirTemp, sa.NbrRectangle(3, 3, 'CELL'), 'MAXIMUM', 'NODATA'))), 1)

In my experience, compound conditional statements in numpy (i.e., using 'numpy.where') require the same treatment.

Hope that helps.

[ChrisSnyder]

Right on...

I think I'm getting into the groove now. Thanks a million Phil!