2

I have a DEM and I want to fill every sink with a Value of 2 meter. If a sink is less than 2 meter, it is filled full and if it is bigger, it is filled until it reach 2 meter.

The problem is that this does not work with the Spatial Analyst Tool Fill. If the sinks are deeper than the Z-Limit, they will not be filled.

After the fill, the original DEM is subtracted from the filled DEM. With this, I want to get a sink raster, which will show the depth of the water if the sinks are filled up to 2 meters.

I also tried the answer of the question asked here (Partially fill a sink - ArcGIS), but it does not really work.

This image shows a sink of my dem. The sink is roughly 10 m deep. enter image description here

Here you can see the area after the method used in (Partially fill a sink - ArcGIS). The dark blue colour has the value 2 and the white 0. That cannot be, because the outer to inner part of the sink has a difference of 10 meters.

enter image description here

This image shows how the raster roughly should look like (I am sorry, I´ve got only Paint at my disposal). You can see clearly how deep the water is at each spot and is filled up to 2 meters. enter image description here

2 Answers 2

2

After part 2 of the answer at (Partially fill a sink - ArcGIS) (Subtract the original DEM from the filled DEM) you have a raster "sinksOnly.tif".

If your original DEM is called "dem.tif" make a raster calculation with the expression

Con("sinksOnly.tif" > 2, "dem.tif" + 2, "dem.tif" + "sinksOnly.tif")

which means:

if the sink is bigger than 2 meters -> fill the original DEM by 2 Meters

else -> fill it "full"

EDIT:

After you updated the question this could be another approach:

  1. Fill the sinks smaller than 2 meter (I think this part of the sollution is appropriate) using

    Con("sinksOnly.tif" > 2, "dem.tif", "dem.tif" + "sinksOnly.tif")

  2. Remove the sinks smaller than 2 meter from your sinksOnly raster using

    Con("sinksOnly.tif" > 2, "sinksOnly.tif", "sinksOnly.tif" - "sinksOnly.tif")

  3. Run a "Segment Mean Shift" calculation over the new "sinksOnly"

  4. Convert the Segmentation result to Polygons (Raster to Polygons)

  5. Convert the DEM to Points (Raster to Points)

  6. Perform a Spatial Join with the Segmentation Polygons to the DEM points. Select "Minimum" which will be the deepest value of the respective sink.

  7. Delete the polygon without sinks.

  8. Convert the result of the spatial join back to a raster with the MIN_GRID (the deepest value of the sink) as raster value.

  9. Perform the final raster calculator:

    Con("dem.tif"-"deepestValue.tif"<2, "dem.tif" + 2 - "dem.tif"-"deepestValue.tif", "dem.tif")

enter image description here

(sorry for the german)

5
  • This has the same problem as the answer in (gis.stackexchange.com/questions/78387/…). In my DEM there are some sinks that have a depth of more than 10 meters and are bigger than 1 km². With this method, the whole sink is flooded/filled by 2 meters. But i want to get only the inner part of it.
    – Tim
    Commented Jan 14, 2019 at 13:58
  • 1
    can you post a screenshot what you mean with "I want to get only the inner part of it"? Commented Jan 14, 2019 at 14:01
  • or is it exactly the same as the "Desired Custom Fill" of the other question that you want to achieve? Commented Jan 14, 2019 at 14:04
  • 1
    I edited images into my post. The second image shows the results of your method and it is the same as in (gis.stackexchange.com/questions/78387/…). Yes it should be exactly the same, if I understand him correctly. Hope it is more clear now.
    – Tim
    Commented Jan 14, 2019 at 14:43
  • I updated the anwer. Commented Jan 14, 2019 at 15:44
1

Picture below shows depressions filled completely (blue) and filled with water up to 0.5 m deep (black):

enter image description here

WORKFLOW FOR DEPTH 0.5 m:

##    FIND DEPRESSIONS
arcpy.gp.Fill_sa("DEM", "C:/SCRATCH/filled", "")
arcpy.gp.RasterCalculator_sa('Con("filled" > "DEM",1)', "C:/SCRATCH/ONE")
##    NUMBER THEM
arcpy.gp.RegionGroup_sa("ONE", "C:/SCRATCH/DEPRESSIONS", "EIGHT", "WITHIN", "NO_LINK", "")
##    FIND SURFACE AND DEEPEST POINT LEVELS
arcpy.gp.ZonalStatistics_sa("DEPRESSIONS", "VALUE", "filled", "C:/SCRATCH/SURFACE", "MEAN", "DATA")
arcpy.gp.ZonalStatistics_sa("DEPRESSIONS", "VALUE", "DEM", "C:/SCRATCH/BOTTOM", "MINIMUM", "DATA")
##    IF DEPRESSION IS SHALLOW KEEP SURFACE LEVEL
##    OTHERWISE CHANGE A SURFACE TO BE 0.5 M ABOVE DEEPEST POINT
arcpy.gp.RasterCalculator_sa('Con(("SURFACE" - "BOTTOM")<0.5,"SURFACE","BOTTOM"+0.5)', "C:/SCRATCH/RL_05")
##    DEFINE REDUCED DEPRESSIONS OUTLINES
arcpy.gp.RasterCalculator_sa('Con("RL_05" >= "DEM",1)', "C:/SCRATCH/reduced")

Replace last calculation by:

arcpy.gp.RasterCalculator_sa('Con("RL_05" >= "DEM","RL_05" - "DEM")', "C:/SCRATCH/DEPTHS")

To get depths raster:

enter image description here

5
  • +1 for region group.. didn't know this tool exists but its perfect for this task Commented Jan 15, 2019 at 8:15
  • 1
    Yes it is a good one, e.g. gis.stackexchange.com/questions/215744/…
    – FelixIP
    Commented Jan 15, 2019 at 8:45
  • Thank you! I´ve got only one short question. Why do you use arcpy.gp and not arcpy.sa?
    – Tim
    Commented Jan 15, 2019 at 9:26
  • It's copy-paste from Results window.
    – FelixIP
    Commented Jan 15, 2019 at 18:45
  • I hope you realize that above approach fills deepest part of depression, = picture in the middle for referred post. To fill secondary etc you need erase Dem unde 1st and repeat process
    – FelixIP
    Commented Jan 15, 2019 at 19:00

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.