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("\nERROR: ArcPy Execute Error occurred.")
print("------------------------------------")
print("All Geoprocessing Messages:")
print(arcpy.GetMessages())
print("\nStatus Messages (0):")
print(arcpy.GetMessages(0))
print("\nWarning Messages (1):")
print(arcpy.GetMessages(1))
print("\nError Messages (2):")
print(arcpy.GetMessages(2))
print("------------------------------------")
except Exception as e:
print(f"\nERROR: 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)
