时间戳

首页 > 解决方案 > 增加折线长度 r sf

r - 增加折线长度 r sf

问题描述

我想找到一种方法将折线从起点增加 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")
)

标签： rsfrgdal

解决方案

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

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

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)
}

因此，在这种情况下，将 sf 问题中的 sf 对象用作line_feature：

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

推荐阅读