增加折线长度 r sf

Question

我想找到一种方法将多段线从起点增加 1000 英尺，从终点增加 1000 英尺。任何帮助是极大的赞赏！下面是一个示例。

structure(
    list(
        TECHNICAL_ = NA_character_,
        geometry = structure(
            list(
                structure(
                    c(812697.360851467,
                      813792.18162311,
                      430678.939150205,
                      425750.102767913),
                    .Dim = c(2L, 2L),
                    class = c("XY", "LINESTRING", "sfg")
                )
            ),
            n_empty = 0L,
            crs = structure(
                list(
                    input = "+init=epsg:2257",
                    wkt = "PROJCRS[\"NAD83 / New Mexico East (ftUS)\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"SPCS83 New Mexico East zone (US Survey feet)\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",31,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-104.333333333333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.999909091,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",541337.5,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8807]],\n        ID[\"EPSG\",15339]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"USA - New Mexico - SPCS - E\"],\n        BBOX[32,-105.72,37,-102.99]]]"
                ),
                class = "crs"
            ),
            class = c("sfc_LINESTRING", "sfc"),
            precision = 0,
            bbox = structure(
                c(xmin = 812697.360851467,
                  ymin = 425750.102767913,
                  xmax = 813792.18162311,
                  ymax = 430678.939150205),
                class = "bbox"
            )
        )
    ),
    row.names = 1L,
    sf_column = "geometry",
    agr = structure(c(TECHNICAL_ = NA_integer_),
                    .Label = c("constant", "aggregate", "identity"),
                    class = "factor"),
    class = c("sf", "tbl_df", "tbl", "data.frame")
)

Answer 1

下面是一个更通用的方法：

首先，您需要一个函数来查找每个端点指向的方向：

st_ends_heading <- function(line)
{
  M <- sf::st_coordinates(line)
  i <- c(2, nrow(M) - 1)
  j <- c(1, -1)
  
  headings <- mapply(i, j, FUN = function(i, j) {
    Ax <- M[i-j,1]
    Ay <- M[i-j,2]
    Bx <- M[i,1]
    By <- M[i,2]
    unname(atan2(Ay-By, Ax-Bx))
  })
  
  return(headings)
}

然后是将线（一端或两端）延伸给定距离（单位距离）的主要功能：

st_extend_line <- function(line, distance, end = "BOTH")
{
  if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
  
  M <- sf::st_coordinates(line)[,1:2]
  keep <- !(end == c("TAIL", "HEAD"))
  
  ends <- c(1, nrow(M))[keep]
  headings <- st_ends_heading(line)[keep]
  distances <- if (length(distance) == 1) rep(distance, 2) else rev(distance[1:2])
  
  M[ends,] <- M[ends,] + distances[keep] * c(cos(headings), sin(headings))
  newline <- sf::st_linestring(M)
  
  # If input is sfc_LINESTRING and not sfg_LINESTRING
  if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
  
  return(newline)
}

所以在这种情况下，在 OP 的问题中使用 sf 对象作为

line_feature

：

sf::st_geometry(line_feature) <- st_extend_line(sf::st_geometry(line_feature), 1000)

Answer 2

我找到了解决方案。首先，使用 st_length() 计算原始线的长度，并附加类似于下面数据表的原始坐标。构建表格后，以下公式将起到作用。

df<-structure(list(x= "abcd", surf_long = 569270.67464245, 
    surf_lat = 645729.112942885, bh_long = 569667.630944828, 
    bh_lat = 645187.374426766, length = 671.606228265491), row.names = c(NA, 
-1L), groups = structure(list(x= "abcd", .rows = structure(list(
    1L), ptype = integer(0), class = c("vctrs_list_of", "vctrs_vctr", 
"list"))), row.names = 1L, class = c("tbl_df", "tbl", "data.frame"
), .drop = TRUE), class = c("grouped_df", "tbl_df", "tbl", "data.frame"
))

现在我们可以按比例增加线长X距离。

df<- df%>%mutate(surf_lat_ext=surf_lat-(bh_lat-surf_lat)*1250/length,
                             surf_long_ext=surf_long-(bh_long-surf_long)*1250/length,
                             bh_lat_ext=bh_lat+(bh_lat-surf_lat)*1250/length,
                             bh_long_ext=bh_long+(bh_long-surf_long)*1250/length)

Answer 3

对于努力让 R 代码工作的人来说，这是另一个选择：使用允许您运行 QGIS 命令的 r 包。

我能够在几行代码中扩展 10,000 多行：

  #Install Required Packages
    install.packages("remotes")
        remotes::install_github("JanCaha/r_package_qgis")
        
  #Import your shapefile with the lines you want to extend
        YOUR_SHAPEFILE_WITH_LINES<-sf::st_read("PATH_TO_SHAPEFILE.shp")
        
  #Extend the lines your desired start and end distance
        result<-qgis::qgis_extendlines(INPUT = YOUR_SHAPEFILE_WITH_LINES,
                                       START_DISTANCE = 100,
                                       END_DISTANCE = 100) %>% st_as_sf()
  #Plot the before and after
    bbox<-sf::st_bbox(result[1:2,])
    ggplot2::ggplot(result)+
      geom_sf(color="blue")+
        xlim(bbox[1], bbox[3])+
        ylim(bbox[2], bbox[4])+
      geom_sf(data=YOUR_SHAPEFILE_WITH_LINES, color="red")