3.5.1 Bug using Slope tool in Python, works in GUI but fails in Notebook.

Question

In ArcGIS Pro 3.5.1, the arcpy.sa.Slope tool is exhibiting a bug where it fails silently when executed from a Python script, producing no output raster and no error messages via arcpy.GetMessages(), despite running perfectly with the exact same parameters (e.g., input DEM, 'DEGREES' measurement, 'GEODESIC' method, Z-factor 1.0) when initiated directly from the ArcGIS Pro GUI. This behavior prevents full automation of workflows requiring slope calculation, necessitating a workaround where users must manually compute the slope in the GUI and then load that pre-calculated raster into their script, even though other arcpy.sa tools like Fill, FlowDirection, and FlowAccumulation function as expected.

u108535adunigisatMorkisch · Answer

Hi Justin

I wonder, why you don't see it - By formatting, Dan means the indentations of your code is lost. Python's interpreter is indentation sensitive. Code won't run.

It'd be much easier, if you linked a pastebin or github link here (it's free). Pasting large code chunks ist kind of an anachronistic agressive anti-pattern (at least for me, I'm a Dev). 🙂

Sincerely,
Sebastian

JustinRego · Answer

Not sure what you mean by formatting is lost? I have print statements I am pasting the code here instead of in the code editor:

import arcpy
from arcpy import env
from arcpy.sa import * # Import Spatial Analyst tools
import os # For os.path.join

import math # Import the standard math module for pi

def calculate_twi(input_dem, output_twi_raster):
"""
Calculates the Topographic Wetness Index (TWI) from a Digital Elevation Model (DEM).

Args:
input_dem (str): Path to the input Digital Elevation Model (DEM) raster.
output_twi_raster (str): Path for the output TWI raster.
Ensure the output path ends with with a .tif extension
or is a geodatabase raster name.
"""
try:
# Check out the Spatial Analyst extension
if arcpy.CheckExtension("Spatial") == "Available":
arcpy.CheckOutExtension("Spatial")
print("Spatial Analyst extension checked out successfully.")
else:
print("ERROR: Spatial Analyst extension is not available. Please ensure it's licensed.")
return

# Set general environment settings
env.overwriteOutput = True
env.scratchWorkspace = arcpy.env.scratchFolder # Use scratch folder for temporary files
print(f"Processing DEM: {input_dem}")

# 1. Fill sinks in the DEM
print("Step 1: Filling DEM sinks...")
fill_dem = Fill(input_dem)
# Save the filled DEM to a temporary file for potential inspection and to ensure it's on disk
temp_fill_dem_path = os.path.join(arcpy.env.scratchFolder, "filled_dem_temp.tif")
fill_dem.save(temp_fill_dem_path)
print("Temporary filled DEM saved.")
print(f"Fill DEM Min: {fill_dem.minimum}, Max: {fill_dem.maximum}, Mean: {fill_dem.mean}")

# --- IMPORTANT: Set Environment settings based on the filled DEM properties ---
# Get properties from the filled DEM
spatial_ref = fill_dem.spatialReference
extent = fill_dem.extent
cell_size_x = arcpy.GetRasterProperties_management(input_dem, "CELLSIZEX").getOutput(0)
cell_size_y = arcpy.GetRasterProperties_management(input_dem, "CELLSIZEY").getOutput(0)
cell_size = float(cell_size_x) # Assuming square cells
print(f"Detected cell size: {cell_size} meters")

# Apply environment settings to ensure consistent processing
env.outputCoordinateSystem = spatial_ref
env.extent = extent # Process within the extent of the filled DEM
env.cellSize = cell_size # Ensure output rasters have this cell size
env.snapRaster = fill_dem # Align outputs to the input DEM's grid
print("Environment settings (Coord System, Extent, Cell Size, Snap Raster) set from filled DEM.")
# --- END Environment settings ---

# 2. Calculate Flow Direction
print("Step 2: Calculating Flow Direction...")
flow_direction = FlowDirection(fill_dem, "NORMAL")
print(f"Flow Direction created.")

# 3. Calculate Flow Accumulation
print("Step 3: Calculating Flow Accumulation...")
flow_accumulation = FlowAccumulation(flow_direction)
print(f"Flow Accumulation Min: {flow_accumulation.minimum}, Max: {flow_accumulation.maximum}, Mean: {flow_accumulation.mean}")

# 4. Calculate Slope (in Degrees) using GEODESIC method and explicit Z-factor
print("Step 4: Calculating Slope (in Degrees) using GEODESIC method and Z-factor 1.0 (reverted to arcpy.sa.Slope)...")
# This is the correct and preferred way to call Slope for Spatial Analyst
slope_degrees = Slope(fill_dem, "DEGREES", "GEODESIC", 1.0) # <<< CORRECTED CALL

print(f"Slope (Degrees) Min: {slope_degrees.minimum}, Max: {slope_degrees.maximum}, Mean: {slope_degrees.mean}")

# 5. Convert Slope from Degrees to Radians
print("Step 5: Converting Slope to Radians...")
slope_radians = Raster(slope_degrees) * (math.pi / 180.0)
print(f"Slope (Radians) Min: {slope_radians.minimum}, Max: {slope_radians.maximum}, Mean: {slope_radians.mean}")

# 6. Calculate Tangent of Slope
print("Step 6: Calculating Tangent of Slope (handling zero slope)...")
tan_slope = Tan(Con(slope_radians == 0, 0.0001, slope_radians))
print(f"Tan(Slope) Min: {tan_slope.minimum}, Max: {tan_slope.maximum}, Mean: {tan_slope.mean}")

# 7. Calculate Specific Catchment Area (a)
print("Step 7: Calculating Specific Catchment Area (a)...")
specific_catchment_area_a = (Raster(flow_accumulation) + 1) * cell_size
print(f"Specific Catchment Area (a) Min: {specific_catchment_area_a.minimum}, Max: {specific_catchment_area_a.maximum}, Mean: {specific_catchment_area_a.mean}")

# 8. Calculate TWI: TWI = ln(a / tan(b))
print("Step 8: Calculating TWI...")
twi_ratio = specific_catchment_area_a / tan_slope
twi_raster = Ln(Con(twi_ratio <= 0, 0.0001, twi_ratio))
print(f"TWI Min: {twi_raster.minimum}, Max: {twi_raster.maximum}, Mean: {twi_raster.mean}")

# Save the final TWI raster
twi_raster.save(output_twi_raster)
print(f"TWI calculation complete. Output saved to: {output_twi_raster}")

except arcpy.ExecuteError:
print("
ERROR: ArcPy Execute Error occurred.")
print("------------------------------------")
print("All Geoprocessing Messages:")
print(arcpy.GetMessages())
print("
Status Messages (0):")
print(arcpy.GetMessages(0))
print("
Warning Messages (1):")
print(arcpy.GetMessages(1))
print("
Error Messages (2):")
print(arcpy.GetMessages(2))
print("------------------------------------")
except Exception as e:
print(f"
ERROR: An unexpected general error occurred: {e}")
finally:
# Check in the Spatial Analyst extension
if arcpy.CheckExtension("Spatial") == "Available":
arcpy.CheckInExtension("Spatial")
print("Spatial Analyst extension checked in.")

# --- How to use the script ---
if __name__ == "__main__":
input_dem_path = r"N:\Projects\2022\225437 BLM AZ Class I, Modeling, & Class III (1.CLT)\Cultural\Graphics-GIS\GIS Data\Imperial_DEM_Mask_30m"
output_twi_path = r"C:\Users\jrego\Documents\ArcGIS\Projects\Topographic Wetness Index\Topographic Wetness Index.gdb\TWI_Result"

calculate_twi(input_dem_path, output_twi_path)

DanPatterson · Answer

python formatting is lost unfortunately.

I would start with showing your print statements (I  assume there were no errors)

As for your input and output file, move them to simple folders without spaces and other flotsam and the output file would be best if it was a *.tif (esri recommended raster format), so r"C:\results\output.tif" as a test

Then check the output folder to see if an output was created... if so, add it to your project (assuming your project is open at this point)