Czy istnieje różnica między algorytmami wypełniania depresji Planchon i Darboux a Wang i Liu? Inna niż prędkość?

Czy ktoś może mi powiedzieć, na podstawie rzeczywistych doświadczeń analitycznych, czy istnieje różnica między tymi dwoma algorytmami wypełniania depresji, innymi niż szybkość, z jaką przetwarzają i wypełniają depresje (zatapiania) w DEM?

Szybki, prosty i wszechstronny algorytm wypełniający wgłębienia cyfrowych modeli elewacji

Olivier Planchon, Frederic Darboux

i

Wydajna metoda identyfikacji i wypełniania wgłębień powierzchniowych w cyfrowych modelach wysokości do analizy i modelowania hydrologicznego

Wang i Liu

Dzięki.

Z teoretycznego punktu widzenia wypełnienie depresji ma tylko jedno rozwiązanie, chociaż istnieje wiele sposobów na uzyskanie tego rozwiązania, dlatego istnieje tak wiele różnych algorytmów wypełniania depresji. Dlatego teoretycznie DEM wypełniony albo Planchon i Darboux, albo Wang i Liu, albo dowolnym innym algorytmem wypełniania depresji, powinien wyglądać później identycznie. Prawdopodobnie jednak tego nie zrobią i istnieje kilka powodów. Po pierwsze, chociaż istnieje tylko jedno rozwiązanie do wypełnienia wgłębienia, istnieje wiele różnych rozwiązań do zastosowania gradientu na płaskiej powierzchni wypełnionego wgłębienia. Oznacza to, że zwykle nie tylko chcemy wypełnić depresję, ale chcemy również wymusić przepływ nad powierzchnią wypełnionej depresji. Zwykle wymaga to dodania bardzo małych gradientów i 1) istnieje wiele różnych strategii (wiele z nich wbudowanych jest bezpośrednio w różne algorytmy wypełniania depresji) i 2) radzenie sobie z tak małymi liczbami często powoduje małe błędy zaokrąglania, które mogą przejawiać się w różnicach między wypełnionymi DEM. Spójrz na ten obraz:

Pokazuje „DEM różnicy” między dwoma DEM, oba wygenerowane ze źródłowej DEM, ale jedna z zagłębieniami wypełnionymi przy użyciu algorytmu Planchon i Darboux, a druga z algorytmem Wanga i Liu. Powinienem powiedzieć, że algorytmy wypełniania depresji były narzędziami z Whitebox GAT i dlatego są różnymi implementacjami algorytmów niż te, które opisałeś w swojej odpowiedzi powyżej. Zauważ, że różnice w DEM są mniejsze niż 0,008 m i że są całkowicie zawarte w obszarach wgłębień topograficznych (tj. Komórki siatki, które nie znajdują się w zagłębieniach, mają dokładnie takie same rzędne jak wejściowe DEM). Mała wartość 8 mm odzwierciedla niewielką wartość stosowaną do wymuszenia przepływu na płaskich powierzchniach pozostawionych przez operację napełniania, a także prawdopodobnie w pewnym stopniu wpływa na skalę błędów zaokrąglania przy reprezentowaniu tak małych liczb za pomocą wartości zmiennoprzecinkowych. Nie widać dwóch wypełnionych elementów DEM wyświetlanych na powyższym obrazku, ale z ich wpisów legendy można stwierdzić, że mają one również dokładnie taki sam zakres wartości wysokości, jak można się spodziewać.

Dlaczego więc w powyższej odpowiedzi zaobserwowałbyś różnice wysokości wzdłuż szczytów i innych obszarów bez depresji w DEM? Myślę, że to naprawdę może sprowadzać się tylko do konkretnej implementacji algorytmu. Prawdopodobnie w narzędziu dzieje się coś, co tłumaczy te różnice i nie ma to związku z faktycznym algorytmem. Nie jest to dla mnie zaskakujące, biorąc pod uwagę różnicę między opisem algorytmu w pracy naukowej a jego faktyczną implementacją w połączeniu ze złożonością sposobu przetwarzania danych wewnątrz GIS. W każdym razie dzięki za zadanie tego bardzo interesującego pytania.

Twoje zdrowie,

Jan

Spróbuję odpowiedzieć na własne pytanie - dun dun.

Użyłem SAGA GIS do zbadania różnic w wypełnionych zlewach za pomocą narzędzia do napełniania opartego na Planchon i Darboux (PD) (oraz ich narzędzia do napełniania na podstawie Wang i Liu (WL) dla 6 różnych zlewisk. (Tutaj pokazuję tylko dwa zestawy wyników - były podobne we wszystkich 6 przełomach). Mówię „oparty”, ponieważ zawsze pojawia się pytanie, czy różnice wynikają z algorytmu, czy z konkretnej implementacji algorytmu.

Wodopochodne DEM wygenerowano przez wycięcie mozaikowych danych NED 30 m przy użyciu USGS dostarczonych kształtowych plików kształtu. Dla każdego podstawowego DEM uruchomiono dwa narzędzia; dla każdego narzędzia jest tylko jedna opcja - minimalne wymuszone nachylenie, które w obu narzędziach zostało ustawione na 0,01.

Po wypełnieniu zlewów użyłem kalkulatora rastrowego do określenia różnic w wynikowych siatkach - różnice te powinny wynikać tylko z różnych zachowań obu algorytmów.

Obrazy przedstawiające różnice lub brak różnic (w zasadzie raster obliczonych różnic) są przedstawione poniżej. Formuła użyta do obliczenia różnic była następująca: (((PD_Filled - WL_Filled) / PD_Filled) * 100) - podaje procentową różnicę dla poszczególnych komórek. Komórki w kolorze szarym wykazują teraz różnicę, komórki w kolorze czerwonym wskazują, że wynikowe podniesienie PD było większe, a komórki bardziej zielone, wskazując, że podniesienie WL było większe.

1st Watershed: Clear Watershed, Wyoming

Oto legenda tych zdjęć:

Różnice mieszczą się w zakresie od -0,0915% do + 0,0910%. Wydaje się, że różnice koncentrują się wokół pików i wąskich kanałów strumieniowych, przy czym algorytm WL jest nieco wyższy w kanałach, a PD nieco wyższy wokół zlokalizowanych pików.

Wyczyść wododział, Wyoming, Zoom 1

Wyczyść wododział, Wyoming, Zoom 2

2nd Watershed: Winnipesaukee River, NH

Oto legenda tych zdjęć:

Rzeka Winnipesaukee, NH, Zoom 1

Różnice wynoszą tylko od -0,323% do + 0,315%. Wydaje się, że różnice koncentrują się wokół pików i wąskich kanałów strumieniowych, przy czym (jak poprzednio) algorytm WL jest nieco wyższy w kanałach, a PD nieco wyższy wokół zlokalizowanych pików.

Łał, myśli? Dla mnie różnice wydają się trywialne i prawdopodobnie nie wpłyną na dalsze obliczenia; ktoś się zgadza? Sprawdzam, kończąc przepływ pracy dla tych sześciu przełomów.

Edycja: więcej informacji. Wydaje się, że algorytm WL prowadzi do szerszych, mniej wyraźnych kanałów, powodując wysokie wartości indeksu topograficznego (mój końcowy zestaw danych pochodnych). Obraz po lewej poniżej to algorytm PD, obraz po prawej to algorytm WL.

Te obrazy pokazują różnicę w indeksie topograficznym w tych samych lokalizacjach - szersze wilgotniejsze obszary (więcej kanałów - bardziej czerwony, wyższy TI) na zdjęciu WL po prawej stronie; węższe kanały (mniej mokry obszar - mniej czerwony, węższy czerwony obszar, niższe TI w obszarze) na zdjęciu PD po lewej.

Dodatkowo, oto w jaki sposób PD radził sobie z depresją i jak WL radził sobie z tym (prawy) - czy zauważyłeś podniesiony pomarańczowy segment (dolny indeks topograficzny) przecinający depresję na wyjściu wypełnionym WL?

Tak więc różnice, choć niewielkie, wydają się ściekać przez dodatkowe analizy.

Oto mój skrypt w języku Python, jeśli ktoś jest zainteresowany:

#! /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:\TomTaggart\DepressionFillingTest\Ran_DEMs"

# This directory is the toplevel directoru (i.e. DEM_8)
INPUTDIR   = "D:\TomTaggart\DepressionFillingTest\Ran_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 "\n\nCurrently 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: %s\n" % 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: %s\n" % 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: %s\n" % 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: %s\n" % e)

            print "Done Depression Filling - Diff Calc"

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


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

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

Na poziomie algorytmicznym oba algorytmy dają takie same wyniki.

Dlaczego możesz mieć różnice?

Reprezentacja danych

Jeśli jeden z algorytmów używa float(32-bit), a inny używa double(64-bit), nie należy oczekiwać, że przyniosą ten sam wynik. Podobnie niektóre implementacje reprezentują wartości zmiennoprzecinkowe, które wykorzystują całkowite typy danych, co może również powodować różnice.

Egzekwowanie drenażu

Jednak oba algorytmy utworzą płaskie obszary, które nie będą drenować, jeśli do określenia kierunków przepływu zastosowana zostanie metoda zlokalizowana.

Planchon i Darboux rozwiązują ten problem, dodając niewielki przyrost do wysokości płaskiej powierzchni w celu wymuszenia drenażu. Jak omówiono w Barnes i in. (2014) w artykule „Skuteczne przypisanie kierunku drenażu nad płaskimi powierzchniami w cyfrowych modelach rastrowych” dodanie tego przyrostu może w rzeczywistości spowodować nienaturalne przekierowanie drenażu poza obszar płaski, jeśli przyrost jest zbyt duży. Rozwiązaniem jest użycie np nextafter. Funkcji.

inne przemyślenia

Wang i Liu (2006) to wariant algorytmu Priority-Flood, jak omówiono w moim artykule „Priority-flood: Optymalny algorytm wypełniania depresji i oznaczania zlewni dla cyfrowych modeli wysokości” .

Priority-Flood ma złożoność czasową dla danych liczb całkowitych i zmiennoprzecinkowych. W swoim artykule zauważyłem, że unikanie umieszczania komórek w kolejce priorytetowej jest dobrym sposobem na zwiększenie wydajności algorytmu. Inni autorzy, tacy jak Zhou i in. (2016) oraz Wei i in. (2018) wykorzystali ten pomysł do dalszego zwiększenia wydajności algorytmu. Kod źródłowy wszystkich tych algorytmów jest dostępny tutaj .

Mając to na uwadze, algorytm Planchon i Darboux (2001) to historia miejsca, w którym nauka zawiodła. Podczas gdy zalew priorytetowy działa w czasie O (N) na danych liczbach całkowitych i czasie O (N log N) na danych zmiennoprzecinkowych, P&D działa w czasie O (N ^1,5 ). Przekłada się to na ogromną różnicę wydajności, która rośnie wykładniczo wraz z rozmiarem DEM:

Do 2001 roku Ehlschlaeger, Vincent, Soille, Beucher, Meyer i Gratin opublikowali wspólnie pięć artykułów szczegółowo opisujących algorytm Priority-Flood. Planchon i Darboux oraz ich recenzenci przeoczyli to wszystko i wymyślili algorytm, który był o rząd wielkości wolniejszy. Jest teraz 2018 i wciąż budujemy lepsze algorytmy, ale P&D jest nadal używany. Myślę, że to niefortunne.