r - 从数据框中同一行中的两个点创建线串

Posted 2023-02-16

【中文标题】r - 从数据框中同一行中的两个点创建线串

【英文标题】：r - Create linestring from two points in same row in dataframe

【发布时间】：2019-01-25 20:42:21

【问题描述】：

我想知道是否有一种方法可以从新几何列中数据框中同一行中给定的两个点创建linestring。换句话说，两个点的经度和纬度在数据框中给出，如下所示：

df <- data.frame(id = c("a", "b"), lon1 = c(1,2), lat1 = c(3,4), lon2 = c(5,6), lat2 = c(7,8))  

其中lon1 和lat1 代表第一个点的坐标，lon2 和lat2 是第二个点的坐标。所需的数据框将有两行和两列 - id 列和 geometry 列。

我试过sf::st_linestring，但似乎这个函数只适用于矩阵。

所需的数据框：

desired_df <- data.frame(id = c("a", "a", "b", "b"), lon = c(1,2,5,6), lat = c(3,4,7,8)) %>% st_as_sf(coords = c("lon", "lat"), dim = "XY") %>% st_set_crs(4236) %>% group_by(id) %>% summarise(geometry = st_union(geometry), do_union = FALSE) %>% st_cast("LINESTRING")

【参考方案1】：

更新 - 2021 年 1 月 30 日

我原来的答案的问题是它没有正确设置边界框。

今天我将使用 sfheaders 和 data.table 的这种方法

library(data.table)
library(sfheaders)

dt <- as.data.table(df)

## To use `sfheaders` the data needs to be in long form

dt1 <- dt[, .(id, lon = lon1, lat = lat1)]
dt2 <- dt[, .(id, lon = lon2, lat = lat2)]

## Add on a 'sequence' variable so we know which one comes first
dt1[, seq := 1L ]
dt2[, seq := 2L ]

## put back together
dt <- rbindlist(list(dt1, dt2), use.names = TRUE)
setorder(dt, id, seq)

sf <- sfheaders::sf_linestring(
  obj = dt
  , x = "lon"
  , y = "lat"
  , linestring_id = "id"
)

sf

# Simple feature collection with 2 features and 1 field
# geometry type:  LINESTRING
# dimension:      XY
# bbox:           xmin: 1 ymin: 3 xmax: 6 ymax: 8
# CRS:            NA
#   id              geometry
# 1  a LINESTRING (1 3, 5 7)
# 2  b LINESTRING (2 4, 6 8)

---

原答案

另一种使用data.table的方法

需要（data.table）

dt <- as.data.table(df)

sf <- dt[
    , 
        geometry <- sf::st_linestring(x = matrix(c(lon1, lon2, lat1, lat2), nrow = 2, ncol = 2))
        geometry <- sf::st_sfc(geometry)
        geometry <- sf::st_sf(geometry = geometry)
    
    , by = id
]

sf::st_as_sf(sf)
# Simple feature collection with 2 features and 1 field
# geometry type:  LINESTRING
# dimension:      XY
# bbox:           xmin: 1 ymin: 3 xmax: 5 ymax: 7
# epsg (SRID):    NA
# proj4string:    NA
# id              geometry
# 1  a LINESTRING (1 3, 5 7)
# 2  b LINESTRING (2 4, 6 8)

【参考方案2】：

我们可以循环使用pmap 并在创建的matrix 上应用st_linestring

library(tidyverse)
library(sf)
out <- pmap(df[-1], ~
               c(...) %>%
                matrix(., , ncol=2, byrow = TRUE) %>% 
                st_linestring) %>%
          reduce(st_sfc) %>%
          mutate(df, geometry = .)

out$geometry
#Geometry set for 2 features 
#geometry type:  LINESTRING
#dimension:      XY
#bbox:           xmin: 1 ymin: 3 xmax: 6 ymax: 8
#epsg (SRID):    NA
#proj4string:    NA
#LINESTRING (1 3, 5 7)
#LINESTRING (2 4, 6 8)

【讨论】：

线串列是否有可能是几何列而不是列表？

@adl 你能检查一下更改是否正常

结果仍然是一个列表。每行在 lstring 列中都有一个列表，而不仅仅是 linestring

@adl 现在怎么样

@adl 我想我之前错过了st_sfc。

【参考方案3】：

df = data.frame(id = c("a", "b"), lon1 = c(1,2), lat1 = c(3,4), lon2 = c(5,6), lat2 = c(7,8))  
df

##   id lon1 lat1 lon2 lat2
## 1  a    1    3    5    7
## 2  b    2    4    6    8

这是另一种方式，通过 WKT：

library(sf)
df$geom = sprintf("LINESTRING(%s %s, %s %s)", df$lon1, df$lat1, df$lon2, df$lat2)
df = st_as_sf(df, wkt = "geom")
df

## Simple feature collection with 2 features and 5 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 1 ymin: 3 xmax: 6 ymax: 8
## CRS:            NA
##   id lon1 lat1 lon2 lat2                  geom
## 1  a    1    3    5    7 LINESTRING (1 3, 5 7)
## 2  b    2    4    6    8 LINESTRING (2 4, 6 8)

【讨论】：

感谢 WKT 的实施。它有点不直观，但非常容易阅读，并且如果您从带有开始/结束列的 df 开始，则需要 FAR FAR FAR 更少的代码。

【参考方案4】：

这个方案也是使用purrr的pmap，得到想要的格式的结果

library(tidyverse)
library(sf) 

df <- data.frame(id = c("a", "b"), lon1 = c(1,2), lat1 = c(3,4), lon2 = c(5,6), lat2 = c(7,8))  

make_line <- function(lon1, lat1, lon2, lat2) 
    st_linestring(matrix(c(lon1, lon2, lat1, lat2), 2, 2))


df %>%
    select(-id) %>% 
    pmap(make_line) %>% 
    st_as_sfc(crs = 4326) %>% 
    tibble(id = df$id, geometry = .) %>% 
    st_sf() 

结果：

Simple feature collection with 2 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 1 ymin: 3 xmax: 6 ymax: 8
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
# A tibble: 2 x 2
  id            geometry
  <fct> <LINESTRING [°]>
1 a           (1 3, 5 7)
2 b           (2 4, 6 8)