如何修复使用Riverdist时“修剪器中的错误 - 没有剩余的线段”?

问题描述 投票:1回答:2

我一直在寻找解决方案来计算采样点和河口/河口之间的距离,跟随河流线。看起来“Riverdist”可以解决我所有的问题,但我从一开始就遇到了麻烦。我在国家层面(南非)工作,但由于河流系统非常庞大和复杂,所以我们一直在尝试单一集水区。

我做了什么 :

  • 我使河流shp文件尽可能简单,我删除了所有不必要的支流
  • 使用预计的CRS
  • 试图用“line2network”加载它
> library(riverdist)

> limpopo <- line2network(path = "/Volumes/Shadowfax/Distribution/simple rivers/Limpopo.shp", layer="Limpopo", tolerance = 100, reproject = NULL,
+              supplyprojection = NULL)

这是我得到的错误消息:trimriver错误(trim =问题,河流=河流):错误 - 导致河网没有剩余的线段

我的shp文件看起来对我很好。

我想在编辑我的shp或保存时我做错了但我不知道是什么。这里是shp的一瞥:Limpopo.shp我必须补充一点,我对R或任何编码都缺乏经验。这是我第一次去“Riverdist”。

希望我能在这里找到帮助!谢谢 !

[编辑]使用MM回答。但我还有一个问题。河段不被认为是连续的(如果有意义的话),当计算点与河口之间的距离时,它会让我遇到麻烦。我还试图让QPI上的shp变得更加简单。

library(riverdist)

# Create your custom function to read the file. 
# Examine the line2network function and modify the lines that cause an issue with your specific shape file

my.custom.line2network =  function (path = ".", layer, tolerance = 100, reproject = NULL, supplyprojection = NULL) {
  sp <- suppressWarnings(rgdal::readOGR(dsn = path, layer = layer, verbose = F))
  if (class(sp) != "SpatialLinesDataFrame") 
    stop("Specified shapefile is not a linear feature.")
  if (is.na(sp@proj4string@projargs) & !is.null(supplyprojection)) 
    sp@proj4string@projargs <- supplyprojection
  if (is.na(sp@proj4string@projargs)) 
    stop("Shapefile projection information is missing.  Use supplyprojection= to specify a Proj.4 projection to use.  If the input shapefile is in WGS84 geographic (long-lat) coordinates, this will be +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 (in double-quotes).  If so, it must also be reprojected using reproject=.")
  proj4 <- strsplit(sp@proj4string@projargs, split = " ")
  projected <- sp::is.projected(sp)
  if (is.null(reproject) & !projected) 
    stop("Distances can only be computed from a projected coordinate system. Use reproject= to specify a Proj.4 projection to use.")
  if (!is.null(reproject)) {
    sp <- sp::spTransform(sp, sp::CRS(reproject))
    proj4 <- strsplit(sp@proj4string@projargs, split = " ")
  }
  units <- "unknown"
  for (i in 1:length(proj4[[1]])) {
    if (proj4[[1]][i] != "") {
      proj4arg <- strsplit(proj4[[1]][i], split = "=")
      if (proj4arg[[1]][1] == "+units") {
        units <- proj4arg[[1]][2]
        cat("\n", "Units:", proj4arg[[1]][2], "\n")
      }
    }
  }
  if (length(sp@lines) > 1) {
    sp_line <- NA
    sp_seg <- NA
    lines <- list()
    j <- 1
    for (i in 1:length(sp@lines)) {
      for (k in 1:length(sp@lines[i][[1]]@Lines)) {
        lines[[j]] <- sp@lines[i][[1]]@Lines[[k]]@coords
        sp_line[j] <- i
        sp_seg[j] <- k
        j <- j + 1
      }
    }
  }
  if (length(sp@lines) == 1) {
    lines <- sp@lines[1][[1]]@Lines
    length <- length(lines)
    lines.new <- list()
    for (i in 1:length) {
      lines.new[[i]] <- lines[[i]]@coords
    }
    lines <- lines.new
    sp_line <- rep(1, length)
    sp_seg <- 1:length
  }
  length <- length(lines)
  rivID <- 1:length
  lineID <- data.frame(rivID, sp_line, sp_seg)
  connections <- calculateconnections(lines = lines, tolerance = tolerance)
  if (any(connections %in% 5:6)) 
    braided <- TRUE
  lengths <- rep(NA, length)
  for (i in 1:length) {
    lengths[i] <- pdisttot(lines[[i]])
  }
  names <- rep(NA, length)
  mouth.seg <- NA
  mouth.vert <- NA
  mouth <- list(mouth.seg, mouth.vert)
  names(mouth) <- c("mouth.seg", "mouth.vert")
  sequenced <- FALSE
  braided <- NA
  cumuldist <- list()
  for (i in 1:length) {
    xy <- lines[[i]]
    n <- dim(xy)[1]
    cumuldist[[i]] <- c(0, cumsum(sqrt(((xy[1:(n - 1), 1] - 
                                           xy[2:n, 1])^2) + ((xy[1:(n - 1),    2] - xy[2:n, 2])^2))))
  }
  out.names <- c("sp", "lineID", "lines", "connections", "lengths", 
                 "names", "mouth", "sequenced", "tolerance", "units", 
                 "braided", "cumuldist")
  out <- list(sp, lineID, lines, connections, lengths, names, 
              mouth, sequenced, tolerance, units, braided, cumuldist)
  names(out) <- out.names
  class(out) <- "rivernetwork"
  length1 <- length(out$lengths)
  suppressMessages(out <- removeduplicates(out))
  length2 <- length(out$lengths)
  if (length2 < length1) 
    cat("\n", "Removed", length1 - length2, "duplicate segments.", "\n")

  # THIS LINE CAUSES ISSUES COMENT IT OUT
  # THIS LINE CAUSES ISSUES COMENT IT OUT
  # suppressMessages(out <- removemicrosegs(out))

  length3 <- length(out$lengths)
  if (length3 < length2) 
    cat("\n", "Removed", length2 - length3, "segments with lengths shorter than the connectivity tolerance.", "\n")
  return(out)
}

limpopo <- my.custom.line2network(path = "/Volumes/Shadowfax/Distribution/simple rivers/Limpopo.shp", 
                                  layer = "Limpopo", 
                                  tolerance = 100, 
                                  reproject = NULL,
                                  supplyprojection = NULL)
#displaying
plot(limpopo)

#Convert XY to river location
library(readxl)
limpopo_eel <- read_excel("/Volumes/Shadowfax/Distribution/Eel_data/limpopo_eel.xlsx")
View(limpopo_eel)

Limpopo_eel_riv <- xy2segvert(x=limpopo_eel$X, y=limpopo_eel$Y, rivers=limpopo)
hist(Limpopo_eel_riv$snapdist, main="snapping distance (m)")


#displaying point data in river location 
zoomtoseg(seg=c(4,3), rivers=limpopo)
points(limpopo_eel$X, limpopo_eel$Y, pch=16, col="red")
riverpoints(seg=Limpopo_eel_riv$seg, vert=Limpopo_eel_riv$vert, rivers=limpopo, pch=15, 
            col="blue")

#Computing a a matrix between all observations
dmat <- riverdistancemat(Limpopo_eel_riv$seg,Limpopo_eel_riv$vert,limpopo)

# starting location: segment 5, vertex 93
# ending location: segment 4, vertex 2420
detectroute(start=4, end=4, rivers=limpopo)
riverdistance(startseg=4, startvert=93, endseg=4, endvert=2420, rivers=limpopo, map=TRUE)
nter code here

它让我这个Erro riverdist when calculating distance

r qgis
2个回答
0
投票

我不熟悉riverdist包或line2network()函数,但是,有可能稍微修改line2network()以强制它在这种情况下工作。函数中有一行导致问题。

# Create your custom function to read the file. 
# Examine the line2network function and modify the lines that cause an issue with your specific shape file

my.custom.line2network =  function (path = ".", layer, tolerance = 100, reproject = NULL, supplyprojection = NULL) {
sp <- suppressWarnings(rgdal::readOGR(dsn = path, layer = layer, verbose = F))
  if (class(sp) != "SpatialLinesDataFrame") 
    stop("Specified shapefile is not a linear feature.")
  if (is.na(sp@proj4string@projargs) & !is.null(supplyprojection)) 
    sp@proj4string@projargs <- supplyprojection
  if (is.na(sp@proj4string@projargs)) 
    stop("Shapefile projection information is missing.  Use supplyprojection= to specify a Proj.4 projection to use.  If the input shapefile is in WGS84 geographic (long-lat) coordinates, this will be +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 (in double-quotes).  If so, it must also be reprojected using reproject=.")
  proj4 <- strsplit(sp@proj4string@projargs, split = " ")
  projected <- sp::is.projected(sp)
  if (is.null(reproject) & !projected) 
    stop("Distances can only be computed from a projected coordinate system. Use reproject= to specify a Proj.4 projection to use.")
  if (!is.null(reproject)) {
    sp <- sp::spTransform(sp, sp::CRS(reproject))
    proj4 <- strsplit(sp@proj4string@projargs, split = " ")
  }
  units <- "unknown"
  for (i in 1:length(proj4[[1]])) {
    if (proj4[[1]][i] != "") {
      proj4arg <- strsplit(proj4[[1]][i], split = "=")
      if (proj4arg[[1]][1] == "+units") {
        units <- proj4arg[[1]][2]
        cat("\n", "Units:", proj4arg[[1]][2], "\n")
      }
    }
  }
  if (length(sp@lines) > 1) {
    sp_line <- NA
    sp_seg <- NA
    lines <- list()
    j <- 1
    for (i in 1:length(sp@lines)) {
      for (k in 1:length(sp@lines[i][[1]]@Lines)) {
        lines[[j]] <- sp@lines[i][[1]]@Lines[[k]]@coords
        sp_line[j] <- i
        sp_seg[j] <- k
        j <- j + 1
      }
    }
  }
  if (length(sp@lines) == 1) {
    lines <- sp@lines[1][[1]]@Lines
    length <- length(lines)
    lines.new <- list()
    for (i in 1:length) {
      lines.new[[i]] <- lines[[i]]@coords
    }
    lines <- lines.new
    sp_line <- rep(1, length)
    sp_seg <- 1:length
  }
  length <- length(lines)
  rivID <- 1:length
  lineID <- data.frame(rivID, sp_line, sp_seg)
  connections <- calculateconnections(lines = lines, tolerance = tolerance)
  if (any(connections %in% 5:6)) 
    braided <- TRUE
  lengths <- rep(NA, length)
  for (i in 1:length) {
    lengths[i] <- pdisttot(lines[[i]])
  }
  names <- rep(NA, length)
  mouth.seg <- NA
  mouth.vert <- NA
  mouth <- list(mouth.seg, mouth.vert)
  names(mouth) <- c("mouth.seg", "mouth.vert")
  sequenced <- FALSE
  braided <- NA
  cumuldist <- list()
  for (i in 1:length) {
    xy <- lines[[i]]
    n <- dim(xy)[1]
    cumuldist[[i]] <- c(0, cumsum(sqrt(((xy[1:(n - 1), 1] - 
                                           xy[2:n, 1])^2) + ((xy[1:(n - 1),    2] - xy[2:n, 2])^2))))
  }
  out.names <- c("sp", "lineID", "lines", "connections", "lengths", 
             "names", "mouth", "sequenced", "tolerance", "units", 
             "braided", "cumuldist")
  out <- list(sp, lineID, lines, connections, lengths, names, 
          mouth, sequenced, tolerance, units, braided, cumuldist)
  names(out) <- out.names
  class(out) <- "rivernetwork"
  length1 <- length(out$lengths)
  suppressMessages(out <- removeduplicates(out))
  length2 <- length(out$lengths)
  if (length2 < length1) 
    cat("\n", "Removed", length1 - length2, "duplicate segments.", "\n")

  # THIS LINE CAUSES ISSUES COMENT IT OUT
  # THIS LINE CAUSES ISSUES COMENT IT OUT
  # suppressMessages(out <- removemicrosegs(out))

  length3 <- length(out$lengths)
  if (length3 < length2) 
    cat("\n", "Removed", length2 - length3, "segments with lengths shorter than the connectivity tolerance.", "\n")
  return(out)
}

现在,您可以使用刚刚创建的新函数读取形状文件

limpopo <- my.custom.line2network(path = "/Volumes/Shadowfax/Distribution/simple rivers/Limpopo.shp", 
                               layer = "Limpopo", 
                               tolerance = 100, 
                               reproject = NULL,
                               supplyprojection = NULL)

您应该获得以下图表:

Limpopo river plot


0
投票

问题解决了。使用预计的CRS保存我的河流时发生了一些错误。投影正确,问题不再出现。

© www.soinside.com 2019 - 2024. All rights reserved.