topic Re: simple Raster Calculator question in Python Questions

simple Raster Calculator question

DanielSchatt — Wed, 01 May 2013 13:04:58 GMT

hi, I'm completely new to Python and I'm just trying to run a basic Raster Calculator operation, just pulling out cells below a certain elevation. The code appears to run just fine but in my output, all the cells have a zero value, i.e. there are no cells with a value of "1" that satisfy the selection. I know this is not the case. Just trying to figure out what's wrong. Really appreciate any help...here's my code:

>>> import arcpy >>> from arcpy import env >>> from arcpy.sa import * >>> arcpy.CheckOutExtension("Spatial") >>> env.workspace = "c:/stuff/SeaLevelRiseAdaptation/10YearStepRasters/Hampton" >>> env.cellSize = 3.28 >>> env.extent = "MyElevation" >>> env.mask = "MyElevation" >>> arcpy.gp.RasterCalculator_sa("MyElevation" < 2.22,"C:/stuff/SeaLevelRiseAdaptation/10YearStepRasters/Hampton/Highest/2052_zone1") <geoprocessing server result object object at 0x257BFF20>

Thanks much!

Re: simple Raster Calculator question

ChrisSnyder — Wed, 01 May 2013 15:38:11 GMT

How about:

import arcpy arcpy.CheckOutExtension("Spatial") arcpy.env.workspace = "c:/stuff/SeaLevelRiseAdaptation/10YearStepRasters/Hampton" arcpy.env.cellSize = 3.28 arcpy.env.extent = "myelevation" #the raster 'myelevation' exists in arcpy.env.workspace, correct? arcpy.env.mask = "myelevation" input1 = arcpy.Raster("myelevation") #you must define your input rasters as a "raster object" in order to use them in the new raster algebra (next line) result1 =  input1 < 2.22 #once they are raster objects, Python "just knows" what to do... result1.save("my_result") #run a .save() to permanently save the result to disk, otherwise the 'result1' raster object gets a funny scratch name and poofs away when you close python

Re: simple Raster Calculator question

DanielSchatt — Thu, 02 May 2013 12:58:19 GMT

thanks Chris!

Re: simple Raster Calculator question

ChrisSnyder — Sat, 11 Dec 2021 06:44:32 GMT

It took me a while to get the syntax of the new map algebra figured out too.... Here's some examples of building some more complex expressions (see the last half):

#Process: Set some envr variables - Note: ExtractByMask tool takes care of this for us
inDem = root + "\\dem"
arcpy.env.extent = areaIdPolyFC
arcpy.env.cellSize = inDem
cellSize = float(arcpy.env.cellSize)

#Process: Extract the input DEM
demGrd = areaIdFolderPath + "\\dem"
demTmp = arcpy.sa.ExtractByMask(inDem, areaIdPolyFC); showGpMessage()
demTmp.save(demGrd)

#Process: Run a Fill
fillTmp = arcpy.sa.Fill(demTmp, ""); showGpMessage()
fillGrd = areaIdFolderPath + "\\fill"
fillTmp.save(fillGrd)

#Process: FlowDir
flowDirTmp = arcpy.sa.FlowDirection(fillTmp, "", ""); showGpMessage()
flowDirGrd = areaIdFolderPath + "\\flowdir"
flowDirTmp.save(flowDirGrd)

#Process: Delete fillTmp since we don't need it anymore    
arcpy.Delete_management(fillTmp); showGpMessage()

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

#Process: Create a streamlink
acreThreshold = .25
acreThresholdInPixels = int(43560 * acreThreshold / cellSize ** 2)
streamLinkTmp = arcpy.sa.StreamLink(arcpy.sa.Con(flowAccTmp > acreThresholdInPixels, 1), flowDirTmp); showGpMessage()
streamLinkGrd = areaIdFolderPath + "\\streamlink"
streamLinkTmp.save(streamLinkGrd)
if streamLinkTmp.hasRAT == False:
    arcpy.BuildRasterAttributeTable_management(streamLinkTmp, ""); showGpMessage()   

#Process: Zonal stats for elev and delete elev since we don't need it anymore 
elevationZoneStatTbl = fgdbPath + "\\elevation_zn_stats"
arcpy.sa.ZonalStatisticsAsTable(streamLinkTmp, "VALUE", demTmp, elevationZoneStatTbl, "DATA", "MIN_MAX_MEAN"); showGpMessage()

#Process: Curvature and zonal stats
curvePlanGrd = areaIdFolderPath + "\\curve_plan"
curveTmp = arcpy.sa.Curvature(demTmp, "", "", curvePlanGrd)
arcpy.Delete_management(curveTmp); showGpMessage()
curvePlanTmp = arcpy.sa.Raster(curvePlanGrd)
curvePlan3Tmp = arcpy.sa.FocalStatistics(curvePlanTmp, arcpy.sa.NbrRectangle(3, 3, "CELL"), "MEAN", "DATA"); showGpMessage()
arcpy.Delete_management(curvePlanTmp); showGpMessage()
curvePlan3Grd = areaIdFolderPath + "\\curve3_plan"
curvePlan3Tmp.save(curvePlan3Grd)
curvePlan3ZoneStatTbl = fgdbPath + "\\curve_plan3_zn_stats"
arcpy.sa.ZonalStatisticsAsTable(streamLinkTmp, "VALUE", curvePlan3Tmp, curvePlan3ZoneStatTbl, "DATA", "MEAN"); showGpMessage()
arcpy.Delete_management(curvePlan3Tmp); showGpMessage()