Solución:
Desde un punto de vista teórico, el llenado de depresión solo tiene una solución, aunque puede haber numerosas formas de llegar a esa solución, razón por la cual existen tantos algoritmos de llenado de depresión diferentes. Por lo tanto, en teoría, un DEM que se llena con Planchon y Darboux o Wang y Liu, o con cualquiera de los otros algoritmos de llenado de depresión, debería verse idéntico después. Sin embargo, es probable que no lo hagan y hay algunas razones por las cuales. Primero, aunque solo hay una solución para llenar una depresión, existen muchas soluciones diferentes para aplicar un gradiente en la superficie plana de la depresión llena. Es decir, generalmente no solo queremos llenar una depresión, sino que también queremos forzar el flujo sobre la superficie de la depresión llena. Eso generalmente implica agregar gradientes muy pequeños y 1) hay muchas estrategias diferentes para hacer esto (muchas de las cuales están integradas directamente en los diversos algoritmos de llenado de depresión) y 2) tratar con números tan pequeños a menudo dará como resultado pequeños errores de redondeo que pueden manifestarse en diferencias entre los DEM llenos. Mira esta imagen:
Muestra la ‘DEM de diferencia’ entre dos DEM, ambos generados a partir del DEM de origen, pero uno con depresiones rellenas con el algoritmo de Planchon y Darboux y el otro con el algoritmo de Wang y Liu. Debo decir que los algoritmos de llenado de depresión eran herramientas de Whitebox GAT y, por lo tanto, son implementaciones de los algoritmos diferentes a las que describió en su respuesta anterior. Observe que las diferencias en los DEM son todas menores de 0,008 my que están completamente contenidas dentro de las áreas de depresiones topográficas (es decir, las celdas de la cuadrícula que no están dentro de depresiones tienen exactamente las mismas elevaciones que el DEM de entrada). El pequeño valor de 8 mm refleja el pequeño valor utilizado para hacer cumplir el flujo en las superficies planas dejadas por la operación de llenado y también es probable que se vea algo afectado por la escala de errores de redondeo cuando se representan números tan pequeños con valores de coma flotante. No puede ver los dos DEM llenos que se muestran en la imagen de arriba, pero puede decir por sus entradas de leyenda que también tienen exactamente el mismo rango de valores de elevación, como era de esperar.
Entonces, ¿por qué observaría las diferencias de elevación a lo largo de los picos y otras áreas sin depresión en el DEM en su respuesta anterior? Creo que realmente solo podría reducirse a la implementación específica del algoritmo. Es probable que algo esté sucediendo dentro de la herramienta para explicar esas diferencias y no está relacionado con el algoritmo real. Esto no es tan sorprendente para mí dada la brecha entre la descripción de un algoritmo en un artículo académico y su implementación real combinada con la complejidad de cómo se manejan los datos internamente dentro de GIS.
Intentaré responder a mi propia pregunta: dun dun dun.
Usé SAGA GIS para examinar las diferencias en las cuencas hidrográficas rellenas usando su herramienta de llenado basada en Planchon y Darboux (PD) (y su herramienta de llenado basada en Wang y Liu (WL) para 6 cuencas hidrográficas diferentes (aquí solo muestro el caso de dos conjuntos de resultados: eran similares en las 6 cuencas hidrográficas) Digo “basado”, porque siempre existe la pregunta de si las diferencias se deben al algoritmo o la implementación específica del algoritmo.
Los DEM de cuencas hidrográficas se generaron recortando datos NED 30 m en mosaico utilizando archivos de formas de cuencas hidrográficas proporcionados por el USGS. Para cada DEM base, se ejecutaron las dos herramientas; solo hay una opción para cada herramienta, la pendiente mínima forzada, que se estableció en ambas herramientas en 0.01.
Una vez que se llenaron las cuencas hidrográficas, utilicé la calculadora ráster para determinar las diferencias en las cuadrículas resultantes; estas diferencias solo deberían deberse a los diferentes comportamientos de los dos algoritmos.
Las imágenes que representan las diferencias o la ausencia de diferencias (básicamente el ráster de diferencias calculadas) se presentan a continuación. La fórmula utilizada para calcular las diferencias fue: (((PD_Filled – WL_Filled) / PD_Filled) * 100): proporcione la diferencia porcentual celda por celda. Las celdas de color gris muestran ahora una diferencia, con celdas de color más rojo que indican que la elevación de PD resultante fue mayor, y celdas de color más verde que indican que la elevación de WL resultante fue mayor.
Primera cuenca: Clear Watershed, Wyoming
Aquí está la leyenda de estas imágenes:
Las diferencias solo oscilan entre -0,0915% y + 0,0910%. Las diferencias parecen estar enfocadas alrededor de picos y canales de flujo estrechos, con el algoritmo WL ligeramente más alto en los canales y PD ligeramente más alto alrededor de los picos localizados.
Cuenca clara, Wyoming, Zoom 1
Cuenca clara, Wyoming, Zoom 2
Segunda cuenca: río Winnipesaukee, NH
Aquí está la leyenda de estas imágenes:
Río Winnipesaukee, NH, Zoom 1
Las diferencias solo oscilan entre el -0,323% y el + 0,315%. Las diferencias parecen estar enfocadas alrededor de picos y canales de flujo estrechos, con (como antes) el algoritmo WL ligeramente más alto en los canales y PD ligeramente más alto alrededor de los picos localizados.
Taaaaaaaaaaaaaaaaaaaaaaan, pensamientos? Para mí, las diferencias que parecen triviales probablemente no afectarán los cálculos posteriores; ¿Alguien está de acuerdo? Lo estoy verificando completando mi flujo de trabajo para estas seis cuencas hidrográficas.
Editar: Más información. Parece que el algoritmo WL conduce a canales más anchos y menos distintos, lo que provoca valores altos de índice topográfico (mi conjunto de datos derivados finales). La imagen de la izquierda a continuación es el algoritmo PD, la imagen de la derecha es el algoritmo WL.
Estas imágenes muestran la diferencia en el índice topográfico en las mismas ubicaciones: áreas más amplias y húmedas (más canal, más rojo, TI más alto) en la imagen WL de la derecha; Canales más estrechos (área menos húmeda – menos rojo, área roja más estrecha, TI más baja en el área) en la imagen de DP a la izquierda.
Además, así es como PD manejó (izquierda) una depresión y cómo WL la manejó (derecha): ¿nota el segmento / línea naranja elevado (índice topográfico inferior) que cruza la depresión en la salida llena de WL?
Así que las diferencias, por pequeñas que sean, parecen filtrarse a través de los análisis adicionales.
Aquí está mi script de Python si alguien está interesado:
#! /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
# ----------------------------------------------------------------------
A nivel algorítmico, los dos algoritmos producirán los mismos resultados.
¿Por qué podría tener diferencias?
Representación de datos
Si uno de sus algoritmos usa float (32 bits) y otros usos double (64 bits), no debe esperar que produzcan el mismo resultado. De manera similar, algunas implementaciones representan valores de punto flotante que usan tipos de datos enteros, lo que también puede resultar en diferencias.
Aplicación de drenaje
Sin embargo, ambos algoritmos producirán áreas planas que no drenarán si se utiliza un método localizado para determinar las direcciones del flujo.
Planchon y Darboux abordan esto agregando un pequeño incremento a la altura del área plana para hacer cumplir el drenaje. Como se discutió en Barnes et al. (2014) “Una asignación eficiente de la dirección de drenaje sobre superficies planas en modelos de elevación digitales ráster”, la adición de este incremento puede causar que el drenaje fuera de un área plana se desvíe de forma no natural si el incremento es demasiado grande. Una solución es utilizar, por ejemplo, el nextafter función.
otros pensamientos
Wang y Liu (2006) es una variante del algoritmo Priority-Flood, como se discutió en mi artículo “Priority-flood: Un algoritmo óptimo de llenado de depresión y etiquetado de cuencas hidrográficas para modelos digitales de elevación”.
Priority-Flood tiene complejidades de tiempo tanto para datos enteros como de punto flotante. En mi artículo, noté que evitar colocar celdas en la cola de prioridad era una buena manera de aumentar el rendimiento del algoritmo. Otros autores como Zhou et al. (2016) y Wei et al. (2018) han utilizado esta idea para aumentar aún más la eficiencia del algoritmo. El código fuente de todos estos algoritmos está disponible aquí.
Con esto en mente, el algoritmo de Planchon y Darboux (2001) es la historia de un lugar donde la ciencia fracasó. Mientras que Priority-Flood trabaja en SOBRE) tiempo en datos enteros y O (N log N) tiempo en datos de punto flotante, P&D trabaja en O (N1,5) tiempo. Esto se traduce en una enorme diferencia de rendimiento que crece exponencialmente con el tamaño del DEM:
En 2001, Ehlschlaeger, Vincent, Soille, Beucher, Meyer y Gratin habían publicado colectivamente cinco artículos que detallaban el algoritmo Priority-Flood. Planchon y Darboux, y sus revisores, se perdieron todo esto e inventaron un algoritmo que era varios órdenes de magnitud más lento. Ahora es 2018 y todavía estamos construyendo mejores algoritmos, pero P&D todavía se está utilizando. Creo que eso es lamentable.
Aquí puedes ver las reseñas y valoraciones de los lectores
Si estás contento con lo expuesto, tienes el poder dejar un escrito acerca de qué te ha gustado de esta división.