Is there a difference between the Planchon & Darboux and Wang and Liu depression filling algorithms? Other than speed?

I will attempt to answer my own question - dun dun dun.

I used SAGA GIS to examine the differences in filled watersheds using their Planchon and Darboux (PD) based filling tool ( and their Wang and Liu (WL) based filling tool for 6 different watersheds. (Here I only show case two sets of results - they were similar across all 6 watersheds) I say "based", because there is always the question as to whether differences are due to the algorithm or the specific implementation of the algorithm.

Watershed DEMs were generated by clipping mosaiced NED 30 m data using USGS provided watershed shapefiles. For each base DEM, the two tools were run; there is only one option for each tool, the minimum enforced slope, which was set in both tools to 0.01.

After the watersheds were filled, I used the raster calculator to determine the differences in the resulting grids - these differences should only be due to the different behaviors of the two algorithms.

Images representing the differences or lack of differences (basically the calculated difference raster) are presented below. The formula used in calculating differences was: (((PD_Filled - WL_Filled) / PD_Filled) * 100) - give the percent difference on a cell by cell basis. Cells grey in color show now difference, with cells redder in color indicating the resulting PD elevation was greater, and cells greener in color indicating the resulting WL elevation was greater.

1st Watershed: Clear Watershed, Wyoming

Here's the legend for these images:

Differences only range from -0.0915% to +0.0910%. Differences seem to be focused around peaks and narrow stream channels, with the WL algorithm slightly higher in the channels and PD slightly higher around the localized peaks.

Clear Watershed, Wyoming, Zoom 1

Clear Watershed, Wyoming, Zoom 2

2nd Watershed: Winnipesaukee River, NH

Here's the legend for these images:

Winnipesaukee River, NH, Zoom 1

Differences only range from -0.323% to +0.315%. Differences seem to be focused around peaks and narrow stream channels, with (as before) the WL algorithm slightly higher in the channels and PD slightly higher around the localized peaks.

Sooooooo, thoughts? To me, the differences seem trivial are probably aren't likely to affect further calculations; anyone agree? I am checking by completing my workflow for these six watersheds.

Edit: More information. It seems that the WL algorithm leads to wider less distinct channels, causing high topographic index values (my final derivative data set). The image on the left below is the PD algorithm, the image on the right is the WL algorithm.

These images show the difference in topographic index at the same locations - wider wetter areas (more channel - redder, higher TI) in the WL pic on the right; narrower channels (less wet area - less red, narrower red area, lower TI in area) in the PD pic on the left.

Additionally, here is how PD handled (left) a depression and how WL handled it (right) - notice the raised orange (lower Topographic index) segment/line crossing through the depression in the WL filled output?

So the differences, however small, do seem to trickle through the additional analyses.

Here is my Python script if anyone is interested:

#! /usr/bin/env python

# ----------------------------------------------------------------------
# Create Fill Algorithm Comparison
# Author: T. Taggart
# ----------------------------------------------------------------------

import os, sys, subprocess, time



# function definitions
def runCommand_logged (cmd, logstd, logerr):
    p = subprocess.call(cmd, stdout=logstd, stderr=logerr)
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# environmental variables/paths
if (os.name == "posix"):
    os.environ["PATH"] += os.pathsep + "/usr/local/bin"
else:
    os.environ["PATH"] += os.pathsep + "C:program files (x86)SAGA-GIS"
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# global variables

WORKDIR    = "D:TomTaggartDepressionFillingTestRan_DEMs"

# This directory is the toplevel directoru (i.e. DEM_8)
INPUTDIR   = "D:TomTaggartDepressionFillingTestRan_DEMs"

STDLOG     = WORKDIR + os.sep + "processing.log"
ERRLOG     = WORKDIR + os.sep + "processing.error.log"
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# open logfiles (append in case files are already existing)
logstd = open(STDLOG, "a")
logerr = open(ERRLOG, "a")
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# initialize
t0      = time.time()
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# loop over files, import them and calculate TWI

# this for loops walks through and identifies all the folder, sub folders, and so on.....and all the files, in the directory
# location that is passed to it - in this case the INPUTDIR
for dirname, dirnames, filenames in os.walk(INPUTDIR):
    # print path to all subdirectories first.
    #for subdirname in dirnames:
        #print os.path.join(dirname, subdirname)

    # print path to all filenames.
    for filename in filenames:
        #print os.path.join(dirname, filename)
        filename_front, fileext = os.path.splitext(filename)
        #print filename
        if filename_front == "w001001":
        #if fileext == ".adf":


            # Resetting the working directory to the current directory
            os.chdir(dirname)

            # Outputting the working directory
            print "nnCurrently in Directory: " + os.getcwd()

            # Creating new Outputs directory
            os.mkdir("Outputs")

            # Checks
            #print dirname + os.sep + filename_front
            #print dirname + os.sep + "Outputs" + os.sep + ".sgrd"

            # IMPORTING Files
            # --------------------------------------------------------------
            cmd = ['saga_cmd', '-f=q', 'io_gdal', 'GDAL: Import Raster',
                   '-FILES', filename,
                   '-GRIDS', dirname + os.sep + "Outputs" + os.sep + filename_front + ".sgrd",
                   #'-SELECT', '1',
                   '-TRANSFORM',
                   '-INTERPOL', '1'
                  ]

            print "Beginning to Import Files"

            try:
                runCommand_logged(cmd, logstd, logerr)

            except Exception, e:
                logerr.write("Exception thrown while processing file: " + filename + "n")
                logerr.write("ERROR: %sn" % e)

            print "Finished importing Files"

            # --------------------------------------------------------------


            # Resetting the working directory to the ouputs directory
            os.chdir(dirname + os.sep + "Outputs")



            # Depression Filling - Wang & Liu
            # --------------------------------------------------------------
            cmd = ['saga_cmd', '-f=q', 'ta_preprocessor', 'Fill Sinks (Wang & Liu)',
                   '-ELEV', filename_front + ".sgrd",
                   '-FILLED',  filename_front + "_WL_filled.sgrd",  # output - NOT optional grid
                   '-FDIR', filename_front + "_WL_filled_Dir.sgrd",  # output - NOT optional grid
                   '-WSHED', filename_front + "_WL_filled_Wshed.sgrd",  # output - NOT optional grid
                   '-MINSLOPE', '0.0100000', 
                               ]

            print "Beginning Depression Filling - Wang & Liu"

            try:
                runCommand_logged(cmd, logstd, logerr)

            except Exception, e:
                logerr.write("Exception thrown while processing file: " + filename + "n")
                logerr.write("ERROR: %sn" % e)

            print "Done Depression Filling - Wang & Liu"


            # Depression Filling - Planchon & Darboux
            # --------------------------------------------------------------
            cmd = ['saga_cmd', '-f=q', 'ta_preprocessor', 'Fill Sinks (Planchon/Darboux, 2001)',
                   '-DEM', filename_front + ".sgrd",
                   '-RESULT',  filename_front + "_PD_filled.sgrd",  # output - NOT optional grid
                   '-MINSLOPE', '0.0100000',
                               ]

            print "Beginning Depression Filling - Planchon & Darboux"

            try:
                runCommand_logged(cmd, logstd, logerr)

            except Exception, e:
                logerr.write("Exception thrown while processing file: " + filename + "n")
                logerr.write("ERROR: %sn" % e)

            print "Done Depression Filling - Planchon & Darboux"

            # Raster Calculator - DIff between Planchon & Darboux and Wang & Liu
            # --------------------------------------------------------------
            cmd = ['saga_cmd', '-f=q', 'grid_calculus', 'Grid Calculator',
                   '-GRIDS', filename_front + "_PD_filled.sgrd",
                   '-XGRIDS', filename_front + "_WL_filled.sgrd",
                   '-RESULT',  filename_front + "_DepFillDiff.sgrd",      # output - NOT optional grid
                   '-FORMULA', "(((g1-h1)/g1)*100)",
                   '-NAME', 'Calculation',
                   '-FNAME',
                   '-TYPE', '8',
                               ]

            print "Depression Filling - Diff Calc"

            try:
                runCommand_logged(cmd, logstd, logerr)

            except Exception, e:
                logerr.write("Exception thrown while processing file: " + filename + "n")
                logerr.write("ERROR: %sn" % e)

            print "Done Depression Filling - Diff Calc"

# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# finalize
logstd.write("nnProcessing finished in " + str(int(time.time() - t0)) + " seconds.n")
logstd.close
logerr.close

# ----------------------------------------------------------------------

At an algorithmic level, the two algorithms will produce the same results.

Why might you be getting differences?

Data Representation

If one of your algorithms uses float (32-bit) and another uses double (64-bit), you should not expect them to produce the same result. Similarly, some implementations represent floating-point values uses integer data types, which may also result in differences.

Drainage Enforcement

However, both algorithms will produce flat areas which will not drain if a localized method is used to determine flow directions.

Planchon and Darboux address this by adding a small increment to the height of the flat area to enforce drainage. As discussed in Barnes et al. (2014)'s paper "An efficient assignment of drainage direction over flat surfaces in raster digital elevation models" the addition of this increment can actually cause drainage outside of a flat area to be unnaturally rerouted if the increment is too large. A solution is to use, e.g., the nextafter function.

Other Thoughts

Wang and Liu (2006) is a variant of the Priority-Flood algorithm, as discussed in my paper "Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models".

Priority-Flood has time complexities for both integer and floating-point data. In my paper, I noted that avoiding placing cells in the priority queue was a good way to increase the algorithm's performance. Other authors such as Zhou et al. (2016) and Wei et al. (2018) have used this idea to further increase the algorithm's efficiency. Source code for all these algorithms is available here.

With this in mind, the Planchon and Darboux (2001) algorithm is a story of a place where science failed. While Priority-Flood works in O(N) time on integer data and O(N log N) time on floating-point data, P&D works in O(N^1.5) time. This translates into a huge performance difference that grows exponentially with the size of the DEM:

By 2001, Ehlschlaeger, Vincent, Soille, Beucher, Meyer, and Gratin had collectively published five papers detailing the Priority-Flood algorithm. Planchon and Darboux, and their reviewers, missed all of this and invented an algorithm which was orders of magnitude slower. It's now 2018 and we're still building better algorithms, but P&D is still being used. I think that's unfortunate.