I'm trying to figure out an automated way to identify mountainslopes from where a rockfall could hit lakes using the hydrological tools of ArcGIS 10.7. As rockfalls do not behave like runoff water, some improvisation is needed to get only those slopes from where rocks would hit the lake directly.

(I've asked a previous question about this, see Create line feature from aspect value in ArcGIS - but while the answer is correct, the computation time is enormous and not suitable for larger applications.)

Let me try to clarify my question with the picture below. The slopes with the blue runoff arrows drain into the lake, but rockfalls from there would not necessarily reach the lake. From the slopes with the yellow runoff arrows, however, they would hit the lake.

enter image description here

So my new idea would be to identify the main inlet of the lake, make it the pour point and exclude all slopes that drain into this pixel (blue arrows) while keeping all slopes that drain into the outlet of the lake (yellow arrows). This way, I would get more or less only the slopes directly adjacent to the lake.

My question now is: How do I identify the "main inlet" of a lake?

I've tried masking the flow accumulation raster with the lake outline line feature and then take the pixel with the second highest value. Unfortunately, the runoff inside the lake area sometimes overlaps with the lake border, so I don't get the inlet but rather a point further downstream.

  • 1
    It's no big deal, but not going to work either. You need a model similar to avalanche run-out distance. Have a look at TauDem.
    – FelixIP
    Feb 27, 2020 at 19:19
  • I've tried for some time now to get TauDem running - no success. With "no big deal" do you mean: Finding the main inlet is no big deal? If so - how could I do it? Feb 28, 2020 at 10:04
  • 1
    See workflow below. However it's not the solution for landslide. I'd use binary slope, 0 - steep, 1 - flat as a weight raster in flow accumulation. Find stop points, where accumulation exceeds certain threshold, exclude catchments above it. Repeat until no stop points found. Remaining area is the one you need. Extremely simplistic, but way better than what you have suggested.
    – FelixIP
    Feb 28, 2020 at 22:53

1 Answer 1


I assume that you have elevation model (DEM), you filled it (Filled), computed flow accumulation (FACC) and derived stream raster (Streams) containing value of 1, where accumulation exceeds certain threshold.

Let's identify lake(s):

arcpy.gp.RasterCalculator_sa('Con("Filled" > "DEM",1)', "C:/SCRATCH/DEPRSNS")
arcpy.gp.RegionGroup_sa("DEPRSNS", "C:/SCRATCH/RGP", "FOUR", "WITHIN", "NO_LINK", "")
arcpy.gp.ExtractByAttributes_sa("RGP", '"COUNT" >=354983', "C:/SCRATCH/LAKE")

enter image description here

Find inlets/outlet points:

arcpy.gp.Expand_sa("LAKE", "C:/SCRATCH/EXPAND", "1", "7989")
arcpy.gp.RasterCalculator_sa('Con(IsNull("LAKE"),"STREAMS" * "FACC" * "EXPAND" / "EXPAND")', "C:/SCRATCH/IOS")
arcpy.RasterToPoint_conversion(in_raster="IOS", out_point_features="C:/SCRATCH/IOS_SHP.shp", raster_field="Value")

enter image description here

Find nearest high flow point inside lake:

arcpy.gp.RasterCalculator_sa('"LAKE" / "LAKE" * "STREAMS" * "FACC"', "C:/SCRATCH/HIGH_FLOW")
arcpy.RasterToPoint_conversion(in_raster="HIGH_FLOW", out_point_features="C:/SCRATCH/INSIDE.shp", raster_field="Value")
arcpy.Near_analysis(in_features="IOS_SHP", near_features="INSIDE", search_radius="", location="NO_LOCATION", angle="NO_ANGLE", method="PLANAR")

As expected, there is only one point (outlet) where flow accumulation at the point is greater than one at the nearest high flow point inside lake:

enter image description here

Yes, you can rank remaining points/inlets (grid_code in above table) and pick one with largest catchment, but this will not solve your bigger problem.

Note: it is save to assume that IO point with greatest flow accumulation is outlet, thus no need to estimate high flow points inside lake.

UPDATE: I tested idea from my comments, using weight raster and stop points. Picture below shows result using following criteria. Weight cell = 1 if it is part of stream and has slope < 20 degrees, 0 otherwise. Consider this a friction surface. 2nd assumption - 100 such points (with cell size of 2 it is equivalent of approx. 200m) will stop rockfall happened upstream.

It is enough to erase drainage areas of most downstream stop points from lakes' catchment. Workflow below calculates such points.

a)Find connected groups of streams:

    arcpy.gp.RasterCalculator_sa('Con("facc" >= 101,1)', "C:/SCRATCH/ONE")
    arcpy.gp.RegionGroup_sa("ONE", "C:/SCRATCH/RGP", "EIGHT", "WITHIN", "NO_LINK", "")
    arcpy.gp.ZonalStatisticsAsTable_sa("RGP", "VALUE", "facc", "C:/SCRATCH/SCRATCH.gdb/ZS_MAX", "DATA", "MAXIMUM")

b) Compute flow accumulation at last point

arcpy.AddField_management(in_table="RGP", field_name="MAX", field_type="LONG", field_precision="", field_scale="", field_length="", field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")
arcpy.CalculateField_management(in_table="ZS_MAX", field="MAX", expression="!MAX!//100*100", expression_type="PYTHON_9.3", code_block="")

and transfer it to output of region group tool: enter image description here

Derive upper ("safe") reaches of streams:

arcpy.Lookup_3d(in_raster="RGP", lookup_field="MAX", out_raster="C:/SCRATCH/REDUCED")
arcpy.gp.RasterCalculator_sa('Con("facc" <= "REDUCED",1)', "C:/SCRATCH/streams.tif")

and compute their catchments:

enter image description here

  • Wow, @FelixIP, thanks so much for your suggestions!! I'm trying to follow your updated approch, but could you make it a little bit clearer for me? I have created the accumulation raster with the binary raster as weight input. Could you explain from where I would get the threshold and how I can calculate the stop points? Or should I put this in another question? Feb 29, 2020 at 12:34
  • There are 2 completely subjective parameters, flat slope and friction length, 20 deg and 200 m. If you take a smaller length, the area around the lake will reduce. Reclass your flow acc. And convert to streams. All dangle points are you first set of stop points.
    – FelixIP
    Feb 29, 2020 at 19:08
  • Parameters of such models need calibration by using data of known rockfalls. Chances are that this overly simple model might not work at all.
    – FelixIP
    Feb 29, 2020 at 19:22
  • Thanks, man! I know that I don't do hard science here, but I need to approximate somehow and you helped a bunch! Feb 29, 2020 at 21:34
  • Good luck, mate. Perhaps search for a script capable of creating same size subcatchments using flow accumulation raster. Doing serie by serie might take very long time.
    – FelixIP
    Feb 29, 2020 at 23:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.