【问题标题】:Path of steepest descent with a known starting and end point具有已知起点和终点的最陡下降路径
【发布时间】:2021-09-03 15:55:35
【问题描述】:

我有一个高程栅格(数字高程模型)。我想找到两个已知点之间最陡峭的下降路径。从一个单元格中,您只能移动到周围的七个单元格(减去前一个单元格),我希望根据最陡峭的下降来选择该单元格。如果所有七个周围的单元格都高于当前单元格,我想选择上升最温和的单元格并继续。所有没有数据的单元格都应该被忽略。我只设法计算了坡度栅格(我知道,这是最简单的步骤),不知道如何从那里开始。

要使用的示例数据集

fp <- sp::SpatialPoints(cbind(333350,3060410)) # starting point
lp <- sp::SpatialPoints(cbind(333100,3060600)) # end point

#raster
    dput(rast)
    new("RasterLayer", file = new(".RasterFile", name = "", datanotation = "FLT4S", 
        byteorder = "little", nodatavalue = -3.4e+38, NAchanged = FALSE, 
        nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, 
        blockrows = 16L, blockcols = 128L, driver = "", open = FALSE), 
        data = new(".SingleLayerData", values = c(NA, NA, NA, NA, 
        NA, 1208.45349121094, 1207.15441894531, 1207.19006347656, 
        1207.94274902344, 1207.89782714844, 1207.94274902344, 1209.11303710938, 
        1210.1640625, 1208.50305175781, 1209.26330566406, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.56213378906, 
        1207.99792480469, 1207.19848632812, 1207.94274902344, 1207.94274902344, 
        1208.01794433594, 1208.43347167969, 1208.52490234375, 1210.61877441406, 
        1209.56994628906, 1212.35241699219, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.55212402344, 1208.19763183594, 
        1207.61181640625, 1207.9072265625, 1208.43347167969, 1208.53259277344, 
        1208.56506347656, 1208.56506347656, 1209.91125488281, 1209.72192382812, 
        1217.07543945312, 1219.02136230469, 1219.72485351562, 1217.35620117188, 
        1214.68359375, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        1208.24816894531, 1208.21142578125, 1207.71325683594, 1207.91076660156, 
        1208.19030761719, 1208.56506347656, 1208.53442382812, 1208.56506347656, 
        1211.30285644531, 1215.82885742188, 1214.61682128906, 1221.67639160156, 
        1220.3544921875, 1219.9892578125, 1217.69470214844, 1221.14685058594, 
        1219.81823730469, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        1208.55249023438, 1208.02697753906, 1208.017578125, 1207.400390625, 
        1208.65942382812, 1208.56506347656, 1208.53625488281, 1208.76940917969, 
        1208.76940917969, 1208.76940917969, 1208.76940917969, 1209.5849609375, 
        1217.29479980469, 1222.43395996094, 1220.53955078125, 1222.69262695312, 
        1221.01965332031, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        1208.05090332031, 1207.9970703125, 1207.96997070312, 1208.05041503906, 
        1208.76940917969, 1208.76940917969, 1208.76940917969, 1208.76940917969, 
        1208.76940917969, 1209.5849609375, 1209.5849609375, 1209.5849609375, 
        1222.4541015625, 1220.85107421875, 1221.82141113281, 1216.05773925781, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.44177246094, 
        1208.52172851562, 1208.26806640625, 1208.46484375, 1208.78295898438, 
        1207.81494140625, 1208.16723632812, 1208.69140625, 1209.5849609375, 
        1209.5849609375, 1209.5849609375, 1209.720703125, 1213.89770507812, 
        1214.76220703125, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, 1209.41796875, 1208.56091308594, 1208.10803222656, 
        1208.53002929688, 1208.42626953125, 1207.88256835938, 1208.81066894531, 
        1208.83459472656, 1209.5849609375, 1209.52807617188, 1209.81579589844, 
        1214.59741210938, 1223.80065917969, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.49841308594, 
        1209.14697265625, 1208.54370117188, 1208.89050292969, 1209.81579589844, 
        1209.81579589844, 1214.78759765625, 1221.84899902344, 1223.90710449219, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, 1208.6904296875, 1208.94995117188, 1209.23315429688, 
        1209.67565917969, 1209.67565917969, 1214.1171875, 1215.07580566406, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, 1208.83178710938, 1208.76831054688, 1209.67565917969, 
        1209.67565917969, 1211.46594238281, 1214.94274902344, 1215.39904785156, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, 1209.25732421875, 1208.86218261719, 1209.048828125, 
        1209.67565917969, 1209.67565917969, 1211.82763671875, 1215.14599609375, 
        1215.67370605469, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, 1209.06079101562, 1209.17163085938, 
        1208.79455566406, 1209.67565917969, 1209.67565917969, 1213.59326171875, 
        1220.58471679688, 1217.48803710938, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1210.30834960938, 
        1209.05725097656, 1210.30346679688, 1209.67565917969, 1209.67565917969, 
        1215.83471679688, 1222.40283203125, 1218.2578125, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, 1209.61926269531, 1210.16052246094, 1210.56823730469, 
        1210.79309082031, 1210.46276855469, 1214.73291015625, 1216.96398925781, 
        1219.07653808594, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, 1209.46716308594, 1209.22875976562, 
        1210.47521972656, 1210.0283203125, 1210.0283203125, 1212.83312988281, 
        1216.43481445312, 1219.17358398438, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1210.00036621094, 
        1209.9521484375, 1209.96520996094, 1210.0283203125, 1210.0283203125, 
        1213.79821777344, 1217.90661621094, 1219.22473144531, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, 1209.54418945312, 1210.26818847656, 1210.0283203125, 
        1210.03039550781, 1210.0283203125, 1212.66833496094, 1216.18420410156, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, 1209.45166015625, 1209.70434570312, 1209.33239746094, 
        1210.0283203125, 1210.0283203125, 1211.15905761719, 1212.17797851562, 
        1214.48999023438, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
        NA, NA, NA, NA, NA, NA, NA, NA, 1209.94604492188, 1210.36206054688, 
        1210.103515625, 1210.0283203125, 1210.0283203125, 1211.44226074219, 
        1214.18041992188, NA, NA), offset = 0, gain = 1, inmemory = TRUE, 
            fromdisk = FALSE, isfactor = FALSE, attributes = list(), 
            haveminmax = TRUE, min = 1207.15441894531, max = 1223.90710449219, 
            band = 1L, unit = "", names = "nakkhu_hosp_10m"), legend = new(".RasterLegend", 
            type = character(0), values = logical(0), color = logical(0), 
            names = logical(0), colortable = logical(0)), title = character(0), 
        extent = new("Extent", xmin = 333003.9801, xmax = 333403.9801, 
            ymin = 3060402.4038, ymax = 3060602.4038), rotated = FALSE, 
        rotation = new(".Rotation", geotrans = numeric(0), transfun = function () 
        NULL), ncols = 40L, nrows = 20L, crs = new("CRS", projargs = "+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs"), 
        history = list(), z = list())

编辑:回应 nniloc 的回答

我认为它会找到距离最短的路径。我想要的不一样。当我将您的方法应用于不同的栅格时,请参见下图。看起来路径在找到较短的路线时会跳过周围的单元格。

【问题讨论】:

  • 如果您有 terra pkg,请查看 ?terra::terrain arg v = 'flowdir' 或 'TPI' 它通常是如何计算的,如详细信息中所述。
  • 谢谢。 terra::terrain 确实很有帮助。但我正在努力根据terra::terrain 的结果定义一条路径。对前进的道路有什么建议吗?
  • 另一种使用topmodel的方法如下所示。
  • 采用deldir 的进一步方法,仍有待完成,但草图就在那里。

标签: r raster terrain


【解决方案1】:

leastcostpathpackage 有一个函数create_lcp,它允许您定义origin 和destination。

编辑:添加基于this vignette 的“障碍”。在示例数据中它没有明显改变路径,所以我不能确定它是否正常工作。

library(leastcostpath)

rast_slope <- create_slope_cs(rast, cost_function = 'tobler', neighbours = 8)

# add barrier anywhere rast is equal to NA
rast_nas <- is.na(rast)
rast_nas[rast_nas == 0] <- NA
rast_barrier <- create_barrier_cs(rast, rast_nas, field = 0, background = 1)


rast_slope_barrier <- rast_slope * rast_barrier
#plot(raster(rast_slope_barrier), col = grey.colors(100))


lcp <- create_lcp(rast_slope_barrier, fp, lp, TRUE)

plot(rast)
lines(lcp)
points(fp)
points(lp)

数据

library(raster)

rast <- new("RasterLayer", file = new(".RasterFile", name = "", datanotation = "FLT4S", 
                                      byteorder = "little", nodatavalue = -3.4e+38, NAchanged = FALSE, 
                                      nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, 
                                      blockrows = 16L, blockcols = 128L, driver = "", open = FALSE), 
            data = new(".SingleLayerData", values = c(NA, NA, NA, NA, 
                                                      NA, 1208.45349121094, 1207.15441894531, 1207.19006347656, 
                                                      1207.94274902344, 1207.89782714844, 1207.94274902344, 1209.11303710938, 
                                                      1210.1640625, 1208.50305175781, 1209.26330566406, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.56213378906, 
                                                      1207.99792480469, 1207.19848632812, 1207.94274902344, 1207.94274902344, 
                                                      1208.01794433594, 1208.43347167969, 1208.52490234375, 1210.61877441406, 
                                                      1209.56994628906, 1212.35241699219, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.55212402344, 1208.19763183594, 
                                                      1207.61181640625, 1207.9072265625, 1208.43347167969, 1208.53259277344, 
                                                      1208.56506347656, 1208.56506347656, 1209.91125488281, 1209.72192382812, 
                                                      1217.07543945312, 1219.02136230469, 1219.72485351562, 1217.35620117188, 
                                                      1214.68359375, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      1208.24816894531, 1208.21142578125, 1207.71325683594, 1207.91076660156, 
                                                      1208.19030761719, 1208.56506347656, 1208.53442382812, 1208.56506347656, 
                                                      1211.30285644531, 1215.82885742188, 1214.61682128906, 1221.67639160156, 
                                                      1220.3544921875, 1219.9892578125, 1217.69470214844, 1221.14685058594, 
                                                      1219.81823730469, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      1208.55249023438, 1208.02697753906, 1208.017578125, 1207.400390625, 
                                                      1208.65942382812, 1208.56506347656, 1208.53625488281, 1208.76940917969, 
                                                      1208.76940917969, 1208.76940917969, 1208.76940917969, 1209.5849609375, 
                                                      1217.29479980469, 1222.43395996094, 1220.53955078125, 1222.69262695312, 
                                                      1221.01965332031, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      1208.05090332031, 1207.9970703125, 1207.96997070312, 1208.05041503906, 
                                                      1208.76940917969, 1208.76940917969, 1208.76940917969, 1208.76940917969, 
                                                      1208.76940917969, 1209.5849609375, 1209.5849609375, 1209.5849609375, 
                                                      1222.4541015625, 1220.85107421875, 1221.82141113281, 1216.05773925781, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.44177246094, 
                                                      1208.52172851562, 1208.26806640625, 1208.46484375, 1208.78295898438, 
                                                      1207.81494140625, 1208.16723632812, 1208.69140625, 1209.5849609375, 
                                                      1209.5849609375, 1209.5849609375, 1209.720703125, 1213.89770507812, 
                                                      1214.76220703125, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, 1209.41796875, 1208.56091308594, 1208.10803222656, 
                                                      1208.53002929688, 1208.42626953125, 1207.88256835938, 1208.81066894531, 
                                                      1208.83459472656, 1209.5849609375, 1209.52807617188, 1209.81579589844, 
                                                      1214.59741210938, 1223.80065917969, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1208.49841308594, 
                                                      1209.14697265625, 1208.54370117188, 1208.89050292969, 1209.81579589844, 
                                                      1209.81579589844, 1214.78759765625, 1221.84899902344, 1223.90710449219, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, 1208.6904296875, 1208.94995117188, 1209.23315429688, 
                                                      1209.67565917969, 1209.67565917969, 1214.1171875, 1215.07580566406, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, 1208.83178710938, 1208.76831054688, 1209.67565917969, 
                                                      1209.67565917969, 1211.46594238281, 1214.94274902344, 1215.39904785156, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, 1209.25732421875, 1208.86218261719, 1209.048828125, 
                                                      1209.67565917969, 1209.67565917969, 1211.82763671875, 1215.14599609375, 
                                                      1215.67370605469, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, 1209.06079101562, 1209.17163085938, 
                                                      1208.79455566406, 1209.67565917969, 1209.67565917969, 1213.59326171875, 
                                                      1220.58471679688, 1217.48803710938, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1210.30834960938, 
                                                      1209.05725097656, 1210.30346679688, 1209.67565917969, 1209.67565917969, 
                                                      1215.83471679688, 1222.40283203125, 1218.2578125, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, 1209.61926269531, 1210.16052246094, 1210.56823730469, 
                                                      1210.79309082031, 1210.46276855469, 1214.73291015625, 1216.96398925781, 
                                                      1219.07653808594, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, 1209.46716308594, 1209.22875976562, 
                                                      1210.47521972656, 1210.0283203125, 1210.0283203125, 1212.83312988281, 
                                                      1216.43481445312, 1219.17358398438, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1210.00036621094, 
                                                      1209.9521484375, 1209.96520996094, 1210.0283203125, 1210.0283203125, 
                                                      1213.79821777344, 1217.90661621094, 1219.22473144531, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, 1209.54418945312, 1210.26818847656, 1210.0283203125, 
                                                      1210.03039550781, 1210.0283203125, 1212.66833496094, 1216.18420410156, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, 1209.45166015625, 1209.70434570312, 1209.33239746094, 
                                                      1210.0283203125, 1210.0283203125, 1211.15905761719, 1212.17797851562, 
                                                      1214.48999023438, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
                                                      NA, NA, NA, NA, NA, NA, NA, NA, 1209.94604492188, 1210.36206054688, 
                                                      1210.103515625, 1210.0283203125, 1210.0283203125, 1211.44226074219, 
                                                      1214.18041992188, NA, NA), offset = 0, gain = 1, inmemory = TRUE, 
                       fromdisk = FALSE, isfactor = FALSE, attributes = list(), 
                       haveminmax = TRUE, min = 1207.15441894531, max = 1223.90710449219, 
                       band = 1L, unit = "", names = "nakkhu_hosp_10m"), legend = new(".RasterLegend", 
                                                                                      type = character(0), values = logical(0), color = logical(0), 
                                                                                      names = logical(0), colortable = logical(0)), title = character(0), 
            extent = new("Extent", xmin = 333003.9801, xmax = 333403.9801, 
                         ymin = 3060402.4038, ymax = 3060602.4038), rotated = FALSE, 
            rotation = new(".Rotation", geotrans = numeric(0), transfun = function () 
              NULL), ncols = 40L, nrows = 20L, crs = new("CRS", projargs = "+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs"), 
            history = list(), z = list())

fp <- sp::SpatialPoints(cbind(333350,3060410)) # starting point
lp <- sp::SpatialPoints(cbind(333100,3060600)) # end point

【讨论】:

  • 编辑添加障碍。很想知道该解决方案是否适用于您更大的数据集。它在较小的栅格上运行,但似乎不会更改 LCP。我可能错过了一步。
  • 它看起来有效。我有几个问题。如果你能解释答案,那就太好了。每个基本问题都有一个 - 与我所寻找的相反,最低成本路径是否代表最温和的斜坡路径?第二个问题 - 我不应该使用 'neighbours = 8' 而不是默认值 16,因为我一次只查看周围的 8 个单元格吗?我尝试了两者,它们给出了不同的路径。
  • 好问题。回复neighbors = 8 你是对的,至少当我阅读文档时,我得出了和你一样的结论。至于你的第二个问题,有点让人费解。我想我意识到我不知道如何在栅格上计算斜率,也不知道斜率如何转换为 D8 流。
  • 如果您尝试计算流动路径,我认为最终您最好按照here、here 和here 的建议使用 GIS 程序。您问题的棘手/非标准部分是预定义的终点。
【解决方案2】:

使用您的数据,这里称为 mamu_rast:

mamu_flow <- terra::terrain(mamu_rast, v = 'flow_dir)
plot(mamu_rast)
plot(mamu_flow, add = TRUE)

在我的情节中,白框清楚地或几乎遵循您上面的最低成本路径,我认为可以从 mamu_flow 中提取,并在 mamu_rast 上重新绘制。在“河流”的尽头似乎有一个轻微的“出口”海拔问题,这可能是原始数据接近的产物。最后一个单元格比河流高......

我没有发布我的图,因为我已经炸毁了 fontconfig,所以所有本地图都显得毫无意义。 flowdir 使用 8 个邻居并且随机断开连接。 16 位邻居可能会指出另一种打破平局的方法,但我不知道。

对于您在上面的响应中有趣的“大弯河”,也许 dput(case_2) 与“flow_dir”一起使用。

由于 'flow_dir' 的输出被认为令人失望,因此使用实际的水电包可能会有所改善,这里是 topmodel,其优点(在许多中)是它已完善并返回更稀缺的栅格:

library(raster)
library(topmodel)

mamu_mtx <- raster::as.matrix(mamu_rast)
mamu_topidx <- topidx(mamu_mtx, resolution=10)
#returns 2 lists,atb (akin to slope) area (drainage to flow by cell)
#one then fishes around for values that seem to make sense
#but it's a river so the fishing is good
manu_top_b3_a1k <- river(mamu_mtx, atb=mamu_topidx$atb,area=mamu_topidx$area, res=10, thatb=3, tharea=1000)
#make it a raster again
mamu_top_b3_a1k_rast <- raster(mamu_top_be_a1k, template=mamu_rast)
#get river cells
river_cells <- Which(mamu_top_b3_a1k_rast, cells=TRUE, na.rm=TRUE)
xyFromCell(mamu_top_b3_a1k_rast, cell=river_cells)
          x       y
[1,] 333079 3060597
[2,] 333089 3060587
[3,] 333099 3060577
[4,] 333149 3060547
[5,] 333159 3060547
[6,] 333249 3060517
[7,] 333309 3060427

从空间点构造空间线的美妙痛苦留给读者。有人可能会争辩说我的出口点不存在。此时可以追加或交换点。我不太愿意与几代水文学家争论这个模型。

另一种方法,假设有人从 JAXA 日本卫星数据门户网站追踪到ALPSMLC30_N027E085_DSM.tif,然后在终端:

gdalwarp -t_srs "+proj=utm +zone=44 +datum=WGS84 +units=m +no_defs" ALPSMLC30_N027E085_DSM.tif nep27_85.tif
# clip to area around hospital
gdalwarp -te 924530.12 3066804.08 925299.38 3067322.88 nep27_85.tif nakhu_hosp_bend30.tif
Creating output file that is 26P x 18L.
Processing nep27_85.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

返回 R:

nakhu_hosp <- raster('./N027E085/N027E085/nakhu_hosp_bend30.tif')
nakhu_hosp <- readAll(nakhu_hosp) #so we have values for dput()

dput(nakhu_hosp2)
new("RasterLayer", file = new(".RasterFile", name = "", datanotation = "INT2S", 
    byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE, 
    nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, 
    blockrows = 18L, blockcols = 26L, driver = "gdal", open = FALSE), 
    data = new(".SingleLayerData", values = c(1286L, 1283L, 1281L, 
    1282L, 1282L, 1282L, 1286L, 1288L, 1287L, 1287L, 1288L, 1287L, 
    1287L, 1295L, 1294L, 1292L, 1292L, 1290L, 1289L, 1289L, 1292L, 
    1293L, 1292L, 1294L, 1298L, 1298L, 1288L, 1284L, 1282L, 1282L, 
    1282L, 1282L, 1284L, 1287L, 1287L, 1287L, 1287L, 1287L, 1287L, 
    1292L, 1292L, 1290L, 1290L, 1288L, 1289L, 1291L, 1292L, 1293L, 
    1293L, 1294L, 1296L, 1300L, 1288L, 1285L, 1283L, 1282L, 1282L, 
    1281L, 1283L, 1285L, 1287L, 1287L, 1286L, 1286L, 1286L, 1288L, 
    1286L, 1287L, 1286L, 1286L, 1287L, 1291L, 1292L, 1292L, 1291L, 
    1291L, 1296L, 1300L, 1287L, 1285L, 1284L, 1284L, 1282L, 1282L, 
    1283L, 1284L, 1285L, 1285L, 1285L, 1286L, 1287L, 1288L, 1286L, 
    1288L, 1286L, 1286L, 1285L, 1287L, 1290L, 1292L, 1291L, 1289L, 
    1294L, 1297L, 1289L, 1286L, 1285L, 1285L, 1284L, 1283L, 1283L, 
    1283L, 1286L, 1286L, 1285L, 1285L, 1287L, 1289L, 1288L, 1288L, 
    1292L, 1288L, 1286L, 1286L, 1287L, 1288L, 1289L, 1289L, 1293L, 
    1295L, 1293L, 1289L, 1287L, 1286L, 1285L, 1286L, 1285L, 1285L, 
    1285L, 1285L, 1285L, 1285L, 1289L, 1291L, 1291L, 1291L, 1293L, 
    1288L, 1286L, 1286L, 1290L, 1290L, 1288L, 1291L, 1293L, 1295L, 
    1299L, 1294L, 1293L, 1292L, 1291L, 1291L, 1290L, 1288L, 1288L, 
    1289L, 1289L, 1290L, 1293L, 1292L, 1294L, 1293L, 1294L, 1290L, 
    1287L, 1286L, 1288L, 1289L, 1288L, 1293L, 1292L, 1293L, 1300L, 
    1299L, 1297L, 1299L, 1297L, 1298L, 1298L, 1297L, 1297L, 1300L, 
    1298L, 1299L, 1297L, 1294L, 1295L, 1295L, 1293L, 1291L, 1286L, 
    1287L, 1289L, 1290L, 1293L, 1293L, 1295L, 1294L, 1301L, 1299L, 
    1300L, 1303L, 1299L, 1301L, 1303L, 1302L, 1301L, 1303L, 1303L, 
    1301L, 1301L, 1297L, 1296L, 1296L, 1294L, 1294L, 1285L, 1287L, 
    1291L, 1292L, 1290L, 1290L, 1293L, 1291L, 1301L, 1299L, 1302L, 
    1303L, 1302L, 1304L, 1305L, 1304L, 1303L, 1303L, 1302L, 1302L, 
    1303L, 1302L, 1302L, 1298L, 1295L, 1294L, 1284L, 1287L, 1288L, 
    1291L, 1293L, 1295L, 1292L, 1292L, 1300L, 1300L, 1302L, 1303L, 
    1301L, 1304L, 1305L, 1304L, 1303L, 1303L, 1303L, 1304L, 1305L, 
    1310L, 1306L, 1302L, 1299L, 1293L, 1289L, 1286L, 1287L, 1287L, 
    1290L, 1291L, 1292L, 1297L, 1301L, 1303L, 1302L, 1302L, 1301L, 
    1305L, 1305L, 1304L, 1304L, 1302L, 1302L, 1307L, 1310L, 1309L, 
    1305L, 1301L, 1299L, 1294L, 1289L, 1286L, 1287L, 1286L, 1289L, 
    1291L, 1295L, 1299L, 1305L, 1303L, 1303L, 1304L, 1301L, 1305L, 
    1305L, 1304L, 1304L, 1305L, 1304L, 1308L, 1312L, 1312L, 1306L, 
    1301L, 1298L, 1295L, 1289L, 1286L, 1286L, 1287L, 1293L, 1293L, 
    1292L, 1296L, 1301L, 1303L, 1303L, 1306L, 1304L, 1306L, 1303L, 
    1304L, 1305L, 1306L, 1309L, 1313L, 1315L, 1311L, 1304L, 1301L, 
    1298L, 1293L, 1288L, 1286L, 1286L, 1288L, 1292L, 1294L, 1294L, 
    1295L, 1305L, 1305L, 1305L, 1305L, 1306L, 1305L, 1304L, 1305L, 
    1305L, 1306L, 1314L, 1316L, 1314L, 1310L, 1305L, 1302L, 1299L, 
    1291L, 1287L, 1288L, 1287L, 1289L, 1290L, 1292L, 1292L, 1295L, 
    1306L, 1307L, 1307L, 1306L, 1306L, 1306L, 1305L, 1305L, 1308L, 
    1312L, 1314L, 1316L, 1315L, 1309L, 1307L, 1303L, 1298L, 1292L, 
    1288L, 1288L, 1288L, 1287L, 1289L, 1290L, 1295L, 1298L, 1305L, 
    1307L, 1306L, 1306L, 1306L, 1310L, 1309L, 1311L, 1314L, 1317L, 
    1319L, 1320L, 1315L, 1313L, 1311L, 1305L, 1301L, 1293L, 1289L, 
    1288L, 1288L, 1287L, 1288L, 1288L, 1289L, 1293L, 1306L, 1306L, 
    1306L, 1310L, 1309L, 1307L, 1310L, 1318L, 1319L, 1320L, 1324L, 
    1322L, 1318L, 1313L, 1311L, 1309L, 1304L, 1294L, 1291L, 1289L, 
    1288L, 1289L, 1289L, 1288L, 1288L, 1289L), offset = 0, gain = 1, 
        inmemory = TRUE, fromdisk = FALSE, isfactor = FALSE, 
        attributes = list(), haveminmax = TRUE, min = 1281L, 
        max = 1324L, band = 1L, unit = "", names = "nakhu_hosp_bend30"), 
    legend = new(".RasterLegend", type = character(0), values = logical(0), 
        color = logical(0), names = logical(0), colortable = logical(0)), 
    title = character(0), extent = new("Extent", xmin = 924530.12, 
        xmax = 925299.38, ymin = 3066804.08, ymax = 3067322.88), 
    rotated = FALSE, rotation = new(".Rotation", geotrans = numeric(0), 
        transfun = function () 
        NULL), ncols = 26L, nrows = 18L, crs = new("CRS", projargs = "+proj=utm +zone=44 +datum=WGS84 +units=m +no_defs"), 
    history = list(), z = list())

summary(nakhu_hosp)
        nakhu_hosp_bend30
Min.                 1281
1st Qu.              1288
Median               1293
3rd Qu.              1303
Max.                 1324
NA's                    0

我们选择一个合理的排水值,比如 1288m,然后创建一个栅格图层,将高于该值的值设置为 NA:

nakhu_drainage <- clamp(nakhu_hosp, upper=1288, useValues = FALSE)
# make matrix for raster::reclassify
m <- c(1281, 1288, 1288)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
drain_1288 <- reclassify(nakhu_drain, rclmat)
drain_1288_poly <- rasterToPolygons(drain_1288, n=8, dissolve=TRUE)
plot(drain_1288_poly)

好吧,关闭,但想去掉那个悬空的两个像素。所以,进入这个对象并获取我想要的东西(很高兴只是一个简短的多边形列表):

drain_poly_I_want <- drainage_poly@polygons[[2]]@Polygons[[1]]
# to SpatialPolygon
drain_SP <- SpatialPolygons(Srl=list(Polygons(list(srl = drain_poly_I_want), ID='1')), proj4string=CRS(as.character(nakhu_hosp@crs)))

library(deldir)
tessel1 <- deldir(x= drain_1288_poly@polygons[[2]]@Polygons[[1]]@coords[,1],
y= drain_1288_poly@polygons[[2]]@Polygons[[1]]@coords[,2]) 

# using clipp in `tile.list`, note doesn't like valid polygons
# i.e. drop either first or last @coord

drain_tiles <- tile.list(tessel1, clipp = list(x = drain_1288_poly@polygons[[2]]@Polygons[[1]]@coords[1:220,1], y = drain_1288_poly@polygons[[2]]@Polygons[[1]]@coords[1:220,2]))
plot(drain_tiles[1:131])

现在沿着外边缘找到一些点,然后画一条线,往返于上述端点。

【讨论】:

  • mamu_flow 只会给出一个带有编码流向的栅格。虽然这是朝着正确方向迈出的第一步,但它不是我正在寻找的最终结果(即使用最陡峭的上升连接起点和终点的折线)。不过谢谢。
  • 感谢您的详细解答。我还没有采用第二种方法。我正在尝试开发一个通用程序,但不确定是否可以调整此程序,因为某些步骤可能需要手动输入。顺便说一句, nakhu_hosp2 栅格在这里代表什么?是不是根据引流的某种编码细胞?
  • 我很确定这只是为了以防nakhu_hosp &lt;- raster::readALL(nakhu 以某种方式感到惊讶,但应该进行编辑。我没有像你那样将 30m 上采样到 10m。对于一般的解决方案,我肯定会建议一条坡度更大的河流。另外,topmodel 有一个sinkfill 函数,对于预处理来说非常方便。
猜你喜欢
  • 2021-09-19
  • 2018-03-20
  • 1970-01-01
  • 1970-01-01
  • 2014-11-28
  • 2016-01-01
  • 2019-05-27
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多