My problem makes use of ArcGIS tools generally and Spatial Analyst in particular, but my question has more to do with hydrological modelling and working with DEMs. My goal is to obtain flow distance rasters for a large number of sites based on a national DEM and sub-catchment dataset, but in my scripted calculations many of the flow distance rasters contain cells with low distances along the boundary or edge of the site catchment. I would like to use flow distance rasters to weigh the impact of agricultural land uses on stream ecosystems.

I have developed an ArcPy workflow for obtaining flow distance rasters for hundreds of sites based on three datasets:

  • Point feature class of sampling sites located along streams (~1,200 and expected to grow)
  • DEM (2m resolution, but re-sampled at 10m or 25m for reasonable computation times)
  • Polygon feature class of sub-catchments (obtained from government source)

To prepare the DEM for flow distance calculations, I first run the Fill and Flow Direction tool. Based on the location of the sites, I query and select the sub-catchments (based on the sub-catchment containing the site) to extract the flow direction raster based on the selected sub-catchments. I then use the Flow Length tool to calculate the flow distance from each cell to the catchment outlet. To correct for the distance between the site and the outlet (because the site maybe a short distance upstream from the outlet of the catchment polygons), I subtract the flow distance at the site from every cell in the catchment, and set any negative values to Null to obtain a raster representing the flow distances within the catchment to the site. Due to the large number of sites and my research needs for the flow distance data, I have re-sampled the DEM to 25m to obtain reasonable computation times at a broad spatial scale.

The problem with computing flow distances en-masse is that some of the flow distance rasters have cells at the catchment edge/boundary with distances that are clearly too small; I suspect that re-sampling the DEM to a coarser resolution causes the flow paths along the catchment boundaries to be routed to a false "outlet" away the correct catchment outlet (i.e., the site). My first attempt to resolve this issue has been to calculate the Euclidean distances for the extracted site catchments to the site and to SetNull any flow distance cells that are less than the Euclidean distance at the same cell. In many cases this correction works nicely, this doesn't really resolve the problem of why such low flow distances occur, nor does it always remove these low-value cells properly. A few catchments are left with large numbers of cells with Null flow distance.

I have attached a script that performs the workflow I described above for all the sites to document the exact steps I've explained above, and a screenshot illustrating an example of a catchment with low flow distance values along the edge and elsewhere. I hope to hear your thoughts on any flaws in my approach or how I might improve the existing workflow.

Screenshot: enter image description here


import arcpy, os, string, csv, sys, time
from arcpy import env


from arcpy.sa import *

# Local machine paths
root = r"D:\Workspace\DEM_analysis" # define top-level dir
gdb = r"D:\Workspace\DEM_analysis\EcologicalData.gdb" # input data
results = r"D:\Workspace\DEM_analysis\Results.gdb" # output results
scratch = r"D:\Workspace\DEM_analysis\Scratch.gdb" # temporary data
##output = r"D:\Workspace\DEM_analysis\outputs"

env.workspace = gdb
env.scratchWorkspace = scratch
env.overwriteOutput = True
env.parallelProcessingFactor = "100%"

subcatch = gdb + "\\" + "Catchments\GAB_EZG_CH"
sites = gdb + "\\" + "Invertebrates\DistinctSites_select"

##dem = root + "\\" + "DEM.gdb\dem_fill_20m"
fd = root + "\\" + "DEM.gdb\dem_fd_20m"
### Binary (0 or 1) raster indicating presence of streams
##riverRaster = gdb + "\\" + "BinaryRiver_2"

##dem = root + "\\" + "DEM25.gdb\dem25_fill"
##fd = root + "\\" + "DEM25.gdb\dem25_fd"
### Binary (0 or 1) raster indicating presence of streams
##riverRaster = gdb + "\\" + "BinRivers_25"

# Make feature layers
subcatch = arcpy.MakeFeatureLayer_management(subcatch, "sc_lyr")
sites = arcpy.MakeFeatureLayer_management(sites, "site_lyr")

# Dictionaries to associate sites, subcatchments, and hierarchy values
site_catch = {} # siteId : EZGNR pairs to associate sites with subcatchments
catch_hier = {} # EZGNR : H1, H2 values to associate subcatchments with catchment hierarchy
site_val = {} # siteId : flow distance to outlet

# Load list of all distinct site IDs
siteID = []

cursor = arcpy.SearchCursor(sites)
for row in cursor:
del row, cursor

# Query each site and load siteId : EZGNR pairs
# Intensive approach since multiple sites may belong to a single subcatchment
for ID in siteID:
    arcpy.SelectLayerByAttribute_management(sites, "NEW_SELECTION", "\"SiteId\" = \'" + str(ID) + "\'")
    arcpy.SelectLayerByLocation_management(subcatch, "CONTAINS", sites)

    cursor = arcpy.SearchCursor(subcatch)
    for row in cursor:
        site_catch[ID] = row.EZG_NR
        del row, cursor
        del cursor

    arcpy.SelectLayerByAttribute_management(subcatch, "CLEAR_SELECTION")
#    print str(site_catch[ID])

# Counters for print output statements
siteCounter = 0
siteFail = []
print "Script setup... done"

write = sys.stdout.write

##siteID = ['CSCF_644203']

for ID in siteID:
    start = time.time()
    write('Site: ' + str(ID) + ' | ')
    siteCounter = siteCounter + 1
    extractRaster = scratch + "\\" + "ext_" + str(ID)
##    nullRaster = scratch + "\\" + "null_" + str(ID)
##    riverSave = scratch + "\\" + "FD_river_" + str(ID)
    finalSave = results + "\\" + "final_" + str(ID)

    outputED = scratch + "\\" + "euclideanDistance_" + str(ID)
    extractSave = scratch + "\\" + "extract_" + str(ID)
    fdSave = scratch + "\\" + "fd_" + str(ID)
    minusSave = scratch + "\\" + "minus_" + str(ID)
    table = scratch + "\\" + "site_" + str(ID)
    outputRaster = scratch + "\\" + "temp_" + str(ID)

    # Select sub-catchments upstream of a sample
    arcpy.SelectLayerByAttribute_management(sites, "NEW_SELECTION", "\"SiteId\" = \'" + str(ID) + "\'")
    arcpy.SelectLayerByLocation_management(subcatch, "CONTAINS", sites)

    cursor = arcpy.SearchCursor(subcatch)
    for row in cursor:
        catch_hier[row.EZG_NR] = (row.H1, row.H2)
        del row, cursor
        del cursor

    h1 = catch_hier[site_catch[ID]][0]
    h2 = catch_hier[site_catch[ID]][1]

    arcpy.SelectLayerByAttribute_management(subcatch, "ADD_TO_SELECTION", '"H1" >= ' + str(h1) + ' AND "H1" < ' + str(h2))
    write('Query -> | ')

##    try:
    extractRaster = arcpy.gp.ExtractByMask_sa(fd, subcatch)

    write('Extract -> | ')

    # Flow distance to stream
##        nullRaster = arcpy.gp.SetNull_sa(riverRaster, extractRaster)
##        river_FD = arcpy.gp.FlowLength_sa(nullRaster, "DOWNSTREAM", "")
##        river_FD = arcpy.Raster(river_FD)
##        river_FD.save(riverSave)
##        write('River FD -> | ')

    # Flow distance to catchment outlet
    initialFD = arcpy.gp.FlowLength_sa(extractRaster, "DOWNSTREAM", "")
    initialFD = arcpy.Raster(initialFD)
##    initialFD.save(fdSave)

    write('Flow -> | ')

    # Correct flow distance to site, not to outlet
    arcpy.ExtractValuesToTable_ga(sites, initialFD, table)
    cursor = arcpy.da.SearchCursor(table, ("Value"))
    for row in cursor:
        site_val[ID] = row[0]
    del cursor

    x = site_val[ID]

    outMinus = Minus(initialFD, x)
    outSetNull = SetNull(outMinus, outMinus, "VALUE < 1")
##    n = scratch + "\\" + "null_" + str(ID)
##    outSetNull.save(n)

    # Set flow lengths less than Euclidean distance to Null
    arcpy.env.mask = initialFD
    arcpy.env.extent = initialFD
    arcpy.env.snapRaster = initialFD

    EucDist = arcpy.gp.EucDistance_sa(sites, outputED, "", "20", "")
    EucDist = Raster(EucDist)
##    e = scratch + "\\" + "ed_" + str(ID)
##    EucDist.save(e)

    outputRaster = Con((outSetNull < EucDist), -9999, outSetNull)
    outFinal = SetNull(outputRaster, outputRaster, "VALUE = -9999")

    write('Clean & Save || ')

    # Clean up and print loop time
    arcpy.SelectLayerByAttribute_management(subcatch, "CLEAR_SELECTION")
    arcpy.SelectLayerByAttribute_management(sites, "CLEAR_SELECTION")

    fdRaster = gdb + "\\" + "DOWNSTREAM"
##    arcpy.Delete_management(nullRaster)

#    arcpy.Delete_management(EucDist)
##        arcpy.Delete_management(outSetNull)

    stop = time.time()
    elapsed = (stop - start)/60
    elapsed = round(elapsed, 1)
    write('Loop ' + str(siteCounter) + ': ' + str(elapsed) + ' minutes\n')

##    except:
##        print "Site FAILED: " + str(ID)
##        siteFail.append(str(ID))


closed as unclear what you're asking by PolyGeo Jun 17 '17 at 7:26

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • I note that you have not yet taken the Tour which is intended to introduce every user of this site to its protocols. Here you seem to providing us with your code and its symptoms rather than extracting a code snippet that illustrates just the component where you are stuck. You may be lucky and find a potential answerer willing to delve into your code but I think you will get much clearer and quicker answers by extracting a code snippet to present. – PolyGeo Jan 17 '17 at 9:53
  • Thank you, but my question is not really about coding. Rather it is about workflow/methods to correct problems in the outputs of an existing ArcGIS tool. The code I did include is just to document the approach I've used and is fully explained in my text, in case anyone would like to check or compare between my explanation of what I've done and the code. – mbcaradima Jan 17 '17 at 10:07
  • If you have a problem with a single tool then this should be tagged arcgis-desktop rather than arcpy, and any ArcPy code included should just generate the data to be input to the tool, so that with known inputs anyone running the tool would see the same issue. – PolyGeo Jan 17 '17 at 10:09
  • 2
    Thinking out aloud, the discrepancy is probably caused because the sub-catchment layer (you sourced from the government) is probably built from a different DEM, or at least if the same DEM a non-resampled one. So the boundary does not align nicely with your resampled DEM and you are "clipping" flow that is essentially flowing out of the catchment into the adjacent catchment? – Hornbydd Jan 17 '17 at 13:05
  • Hi Hornbydd, I think that explanation makes sense. In that case, one solution would be to derive a new sub-catchment dataset based on a coarser DEM for flow path calculations. Another would be to derive the watersheds straight from a resampled DEM, however I have many sites that do not produce correct watersheds according to known river network data. I'll have to sit on these potential solutions some more, because any solution would have to work for around ~1,200 sites in a topographically varied country (alpine and flatland areas). – mbcaradima Jan 17 '17 at 13:38