try: import arcpy, os, sys, traceback from arcpy.sa import * arcpy.CheckOutExtension("Spatial") degreesSlope = Raster(r"C:\gtemp\FS_test\slope_degree") totalSoilThickness = Raster(r"C:\gtemp\FS_test\soildepth") saturatedSoilThickness = Raster(r"C:\gtemp\FS_test\dw") treeRootStrength = Raster(r"C:\gtemp\FS_test\rootcohesion") treeSurcharge = Raster(r"C:\gtemp\FS_test\surcharge") moistSoilUnitWeight = Raster(r"C:\gtemp\FS_test\moist_uw") numerator = "numerator" denominator = "denominator" #Constants... soilCohesion = 10.0 effectiveInternalAngleofFriction = 10.0 drySoilUnitWeight = 10.0 saturatedSoilUnitWeight = 10.0 waterUnitWeight = 10.0 rastercalc = r"C:\gtemp\FS_test\dwrastercalc" arcpy.env.cellSize = "MAXOF" arcpy.env.extent = "MAXOF" arcpy.gp.RasterCalculator_sa((treeRootStrength + soilCohesion + Square(Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness)+ (saturatedSoilUnitWeight-waterUnitWeight)*saturatedSoilThickness)*Tan(effectiveInternalAngleofFriction)), numerator) arcpy.gp.RasterCalculator_sa((Sin(degreesSlope)*Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness) + drySoilUnitWeight * saturatedSoilThickness), denominator) arcpy.gp.RasterCalculator_sa(numerator/denominator, rastercalc ) arcpy.AddMessage("Done!") except arcpy.ExecuteError: msgs = arcpy.GetMessages(2) arcpy.AddError(msgs) print msgs except: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1]) msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n" arcpy.AddError(pymsg) arcpy.AddError(msgs) print pymsg + "\n" print msgs
try: import arcpy, os, sys, traceback from arcpy.sa import * from arcpy import env arcpy.CheckOutExtension("Spatial") arcpy.AddMessage("\n\nCalculate Slope Factor of Safety") #Grids degreesSlope = Raster(r"C:\gtemp\FS_test\slope_degree") totalSoilThickness = Raster(r"C:\gtemp\FS_test\soildepth") saturatedSoilThickness = Raster(r"C:\gtemp\FS_test\dw") treeRootStrength = Raster(r"C:\gtemp\FS_test\rootcohesion") treeSurcharge = Raster(r"C:\gtemp\FS_test\surcharge") moistSoilUnitWeight = Raster(r"C:\gtemp\FS_test\moist_uw") numerator = "numerator" denominator = "denominator" result = "result" #Constants... soilCohesion = 10.0 effectiveInternalAngleofFriction = 10.0 drySoilUnitWeight = 10.0 saturatedSoilUnitWeight = 10.0 waterUnitWeight = 10.0 rastercalc = r"C:\gtemp\FS_test\dwrastercalc" arcpy.env.cellSize = "MAXOF" arcpy.env.extent = "MAXOF" numerator =((treeRootStrength + soilCohesion + Square(Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness)+ (saturatedSoilUnitWeight-waterUnitWeight)*saturatedSoilThickness)*Tan(effectiveInternalAngleofFriction))) denominator =((Sin(degreesSlope)*Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness) + drySoilUnitWeight * saturatedSoilThickness), denominator) result =(numerator/denominator) result.save(rastercalc) # Process: Raster Calculator #arcpy.gp.RasterCalculator_sa("treeRootStrength + soilCohesion + Square(Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness)+ (saturatedSoilUnitWeigh-waterUnitWeight)*saturatedSoilThickness)*Tan(effectiveInternalAngleofFriction))/((Sin(degreesSlope)*Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness) + drySoilUnitWeight * saturatedSoilThickness))", rastercalc) arcpy.AddMessage("Done!") except arcpy.ExecuteError: msgs = arcpy.GetMessages(2) arcpy.AddError(msgs) print msgs except: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1]) msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n" arcpy.AddError(pymsg) arcpy.AddError(msgs) print pymsg + "\n" print msgs
OK, I removed the raster calculator functions from the script. According to the literature I can use raster calculator like expression (+,-,/...) as long as I cast the raster as 'rasters'. I can get something like the below to work without error, but I cant get passed 'numerator' in the code below without getting the cell size and extent not set error.
x1 =treeRootStrength + soilCohesion x2 = Square(Cos(degreesSlope)) * treeSurcharge + moistSoilUnitWeight * (totalSoilThickness - saturatedSoilThickness) ...
try: import arcpy, os, sys, traceback from arcpy.sa import * from arcpy import env arcpy.SetProduct("ArcInfo") arcpy.CheckOutExtension("Spatial") arcpy.AddMessage("\n\nCalculate Slope Factor of Safety") arcpy.AddMessage("beta by Gerry Gabrisch, 5-4-2012") #Get input parameters using ArcToolbox tool degreesSlope = arcpy.GetParameterAsText(0) totalSoilThickness = arcpy.GetParameterAsText(1) saturatedSoilThickness = arcpy.GetParameterAsText(2) treeRootStrength = arcpy.GetParameterAsText(3) treeSurcharge = arcpy.GetParameterAsText(4) soilCohesion = float(arcpy.GetParameterAsText(5)) effectiveInternalAngleofFriction = float(arcpy.GetParameterAsText(6)) drySoilUnitWeight = float(arcpy.GetParameterAsText(7)) moistSoilUnitWeight = arcpy.GetParameterAsText(8) saturatedSoilUnitWeight = float(arcpy.GetParameterAsText(9)) waterUnitWeight = float(arcpy.GetParameterAsText(10)) rastercalc = arcpy.GetParameterAsText(11) ##Use For hardcoding ##Grids #degreesSlope = r"C:\gtemp\FS_test\slope_degree" #totalSoilThickness = r"C:\gtemp\FS_test\soildepth" #saturatedSoilThickness = r"C:\gtemp\FS_test\dw" #treeRootStrength = r"C:\gtemp\FS_test\rootcohesion" #treeSurcharge = r"C:\gtemp\FS_test\surcharge" #moistSoilUnitWeight = r"C:\gtemp\FS_test\moist_uw" ##Constants... #soilCohesion = 10.0 #effectiveInternalAngleofFriction = 10.0 #drySoilUnitWeight = 10.0 #saturatedSoilUnitWeight = 10.0 #waterUnitWeight = 10.0 #Set the output cell size, mask, and cell registration... arcpy.env.cellSize = degreesSlope arcpy.env.extent = degreesSlope arcpy.env.snapeRaster = degreesSlope #Cast input rasters to in-memory raster layers... degreesSlope = Raster(degreesSlope) totalSoilThickness = Raster(totalSoilThickness) saturatedSoilThickness = Raster(saturatedSoilThickness) treeRootStrength = Raster(treeRootStrength) treeSurcharge = Raster(treeSurcharge) moistSoilUnitWeight = Raster(moistSoilUnitWeight) numerator = "numerator" denominator = "denominator" result = "result" numerator =(treeRootStrength + soilCohesion + Square(Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness)+ (saturatedSoilUnitWeight-waterUnitWeight)*saturatedSoilThickness)*Tan(effectiveInternalAngleofFriction)) denominator =((Sin(degreesSlope)*Cos(degreesSlope))*(treeSurcharge + moistSoilUnitWeight*(totalSoilThickness - saturatedSoilThickness) + drySoilUnitWeight * saturatedSoilThickness)) #Do the last division.... result = numerator/denominator #Save the in-memory results to disk... result.save(rastercalc) arcpy.AddMessage("Done!") #Catch ESRI errors here... except arcpy.ExecuteError: msgs = arcpy.GetMessages(2) arcpy.AddError(msgs) print msgs #Catch Python errors here... except: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1]) msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n" arcpy.AddError(pymsg) arcpy.AddError(msgs) print pymsg + "\n" print msgs