Introduction
I wrote the Runoff.py Python program as part of the Dumpster Project. The output of the program is a raster that shows the most direct downhill paths from given source points. This was used to visualize where the pollution from leaking dumpster goes. ArcMap did not seem to have this specific functionality readily available. I decided to use this as an opportunity to learn Python.
An interesting challenge was figuring out a way to convert the location of each source point from state plane coordinates to pixel coordinates. This was accomplished on lines 49 and 50 of the program. The pixel column coordinate is the quotient from dividing the distance the source point is from the left edge of the raster by the average width of a pixel. The pixel row coordinate is the same principle except it is the quotient after the distance from the top edge divided by average pixel height.
column = int((easting - left)//width) row = int((top - northing)//height)
The Runoff.py program uses the value from a D8 flow raster at a given position, calculates the direction of flow travel, and then moves to the adjacent pixel in that direction. The program repeats this process until it goes off the edge of the flow raster or detects an infinite loop.
Output
The output of the program will resemble a network of lines such as the ones below. The hill-shade, legend and scale bar were added separately.
Runoff.py
import arcpy import math import numpy as np project_db = "F:/dumpster/sm_dumpster.gdb" flow_file = "flow_clipped" # Name of D8 flow raster. dumpster_file = "dumpsters_simplified" # Name of source point feature class. runoff_file = "dumpster_runoff_simp" # Name of the runoff raster. arcpy.env.overwriteOutput = True arcpy.env.workspace = project_db arcpy.env.outputCoordinateSystem = flow_file flow_raster = arcpy.Raster(flow_file) flow_array = arcpy.RasterToNumPyArray(flow_raster, nodata_to_value=0) # Prevent accessing indicies that are out of bounds. def check_bounds(row, col): """Verify if (row, col) are within the boundaries of flow_raster. """ # https://www.naftaliharris.com/blog/python-multiline-if-statements/ if ( 0 <= row <= flow_raster.height - 1 and 0 <= col <= flow_raster.width - 1 ): return True else: return False class SourcePoint: @staticmethod def convert_row_col(easting, northing): """Convert state plane coordinates to corresponding row and column location for the pixel that the coordinates fall within. """ height = flow_raster.meanCellHeight width = flow_raster.meanCellWidth left = flow_raster.extent.XMin top = flow_raster.extent.YMax column = int((easting - left)//width) row = int((top - northing)//height) return (row, column) def __init__(self, obj_id, state_coords): self.id = obj_id self.cell_value = 1 self.x, self.y = state_coords self.row, self.col = SourcePoint.convert_row_col(self.x, self.y) class DownstreamTrace: """Trace runoff from a given SourcePoint. """ def __init__(self, source_point): self.runoff_array = np.full( [flow_raster.height, flow_raster.width], 0, dtype=np.uint8) self.row = source_point.row self.col = source_point.col self.cell_value = source_point.cell_value while True: # Tracing of runoff ends if (row, col) is out of bounds. This # happens when runoff flows off the edge of the DEM. if check_bounds(self.row, self.col): # Check if a cell in runoff_array has flow. This is to prevent # an infinite flow loop. A cell with value greater than zero # has flow. if self.runoff_array[self.row][self.col] == 0: self.runoff_array[self.row][self.col] = \ source_point.cell_value self.flow_value = flow_array[self.row][self.col] # There is flow when flow_value is greather than zero. # A value of zero is used to indicate no data/no flow. # Checking is necessary to prevent a mathematical error # when taking the log of flow_value. if self.flow_value > 0: # The flow_directions tuple contains tuples # representing (Δrow, Δcol). This is used for # converting the D8 flow_value to the next location. flow_directions = ( (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1,-1), (-1, 0), (-1, 1) ) self.flow_dir_index = int(math.log(self.flow_value, 2)) self.delta_row, self.delta_col = flow_directions[ self.flow_dir_index] # Update row and col to the next cell. The process will # repeat itself until flow stops, loops back on itself # or the runoff continues off the edge of the DEM. self.row = self.row + self.delta_row self.col = self.col + self.delta_col else: break else: break def __add__(self, other): self.runoff_array = self.runoff_array + other.runoff_array def save_raster(self, file): self.runoff_raster = arcpy.NumPyArrayToRaster( self.runoff_array, flow_raster.extent.lowerLeft, flow_raster.meanCellWidth, flow_raster.meanCellHeight, value_to_nodata=0) self.runoff_raster.save(file) search_c = arcpy.da.SearchCursor(dumpster_file, ["OBJECTID", "total", "SHAPE@XY"]) dumpster_list = [] for r in search_c: obj_id, total, coords = r dumpster_point = SourcePoint(obj_id, coords) dumpster_point.cell_value = total dumpster_list.append(dumpster_point) dumpster_runoff = DownstreamTrace(dumpster_list.pop()) for x in dumpster_list: dumpster_runoff + DownstreamTrace(x) dumpster_runoff.save_raster(runoff_file)