r - 提取简单特征数据框中点的坐标
问题描述
我有 2 条线,我需要提取这 2 条线相交的坐标。我正在尝试使用 R 中的 sf 包来做到这一点。我可以做交叉点,我可以绘制它们,但我无法提取坐标。当每条线只有 1 个相交点时,它确实有效,但是当有多个(并且在实际数据中有很多)时,相交点的坐标存储在每条线的单个字段中。
inters_pt 看起来像:类几何 A c(3, 3, 3.2, 4.2) B c(1.5, 1.5) A c(2, 4)
我收到消息:
st_coordinates.sfc(st_geometry(x)) 中的错误:未针对 sfc_GEOMETRY 类的对象实现
我想提取每个交叉点的坐标并将属性保留在“类”列中。
输出应如下所示: 几何类 A (3, 3.2) A (3, 4.2) B (1.5, 1.5) A (2, 4)
虽然最初我认为它看起来像这样:几何类 A (3, 3) A (3.2, 4.2)
我错了。
搜索让我来到这里: https ://github.com/r-spatial/sf/issues/114
但是那里的解决方案要么掉点(st_cast(mp,“POINT”),失去属性'class'或抛出以下错误:
vapply(lst, class, rep(NA_character_, 3)) 中的错误:值的长度必须为 3,但 FUN(X[ 1 ]) 结果的长度为 2
#First create my data, I create points, and convert them into SpatialLines
lines1.df <- data.frame(
x = c(1,2,3,3,3,3,1,2,3),
y = c(1,2,2,3,4,5,4,4,4),
id = c(rep("A",6), rep("B",3))
)
#with Kyle Walker's functions I convert the points to lines
#https://rpubs.com/walkerke/points_to_line
library(sp)
library(maptools)
points_to_line <- function(data, long, lat, id_field = NULL, sort_field = NULL) {
# Convert to SpatialPointsDataFrame
coordinates(data) <- c(long, lat)
# If there is a sort field...
if (!is.null(sort_field)) {
if (!is.null(id_field)) {
data <- data[order(data[[id_field]], data[[sort_field]]), ]
} else {
data <- data[order(data[[sort_field]]), ]
}
}
# If there is only one path...
if (is.null(id_field)) {
lines <- SpatialLines(list(Lines(list(Line(data)), "id")))
return(lines)
# Now, if we have multiple lines...
} else if (!is.null(id_field)) {
# Split into a list by ID field
paths <- sp::split(data, data[[id_field]])
sp_lines <- SpatialLines(list(Lines(list(Line(paths[[1]])), "line1")))
# I like for loops, what can I say...
for (p in 2:length(paths)) {
id <- paste0("line", as.character(p))
l <- SpatialLines(list(Lines(list(Line(paths[[p]])), id)))
sp_lines <- spRbind(sp_lines, l)
}
return(sp_lines)
}
}
lines1 <- points_to_line(data = lines1.df,
long = "x",
lat = "y",
id_field = "id")
proj4string(lines1) <- CRS("+proj=utm +zone=32 +datum=WGS84")
library(sf)
lines1_sf <- st_as_sf(lines1,
crs = "+proj=utm +zone=32 +datum=WGS84")
lines2.df <- data.frame(
x = c(2,2,2.5,3.5,3.5,2.5,1,2),
y = c(3,4.2,4.2,4.2,3.2,3.2,2,1),
id = c(rep("A",6), rep("B",2))
)
lines2 <- points_to_line(data = lines2.df,
long = "x",
lat = "y",
id_field = "id")
proj4string(lines2) <- CRS("+proj=utm +zone=32 +datum=WGS84")
lines2_sf <- st_as_sf(lines2,
crs = "+proj=utm +zone=32 +datum=WGS84")
#add attributes
lines2_sf$class <- c("A", "B")
#plot both lines
library(ggplot2)
ggplot(lines1_sf) +
geom_sf(aes(color="red")) +
geom_sf(data=lines2_sf) +
coord_sf(crs = "+proj=utm +zone=32 +datum=WGS84") +
theme_bw()
#find intersecting points
#intersect
inters_pt <- st_intersection(lines2_sf, lines1_sf)
#plot shows the intersecting points
ggplot(lines1_sf) +
geom_sf(aes(color="red")) +
geom_sf(data=lines2_sf) +
geom_sf(data=inters_pt) +
coord_sf(crs = "+proj=utm +zone=32 +datum=WGS84") +
theme_bw()
#extract the coordinates from the geometry column
st_coordinates(inters_pt) #<-- gives error
st_cast(inters_pt, "POINT") #removes points
y = as(inters_pt, "Spatial") #loses the 'class' attribute
SpatialPoints(y)
st_cast(st_sfc(inters_pt), "POINT", group_or_split = FALSE) #error
解决方案
首先,值得使用 sf 从点创建线,这减少了所需的代码量:
library(sf)
library(dplyr)
#First create my data, I create points, and convert them into SpatialLines
lines1.df <- data.frame(
x = c(1,2,3,3,3,3,1,2,3),
y = c(1,2,2,3,4,5,4,4,4),
id = c(rep("A",6), rep("B",3))
)
lines2.df <- data.frame(
x = c(2,2,2.5,3.5,3.5,2.5,1,2),
y = c(3,4.2,4.2,4.2,3.2,3.2,2,1),
id = c(rep("A",6), rep("B",2))
)
points_to_line <- function(data, group = "id"){
data <- data %>%
group_by_at(group) %>%
summarise(do_union = FALSE) %>%
st_cast("LINESTRING") %>%
ungroup %>%
select(-do_union)
}
crs <- "+proj=utm +zone=32 +datum=WGS84"
coords <- c("x", "y")
lines1_sf <- st_as_sf(lines1.df, coords = coords, crs = crs) %>%
points_to_line()
lines2_sf <- st_as_sf(lines2.df, coords = coords, crs = crs) %>%
points_to_line()
#add attributes
lines2_sf$class <- c("A", "B")
#find intersecting points
#intersect
inters_pt <- st_intersection(lines1_sf, lines2_sf)
在这里,您需要先转换为 MULTIPOINT 然后再转换为 POINT
#extract the coordinates from the geometry column
inters_pt <- inters_pt %>%
st_cast("MULTIPOINT") %>%
st_cast("POINT")
> st_coordinates(inters_pt)
# X Y
#1 3.0 3.2
#2 3.0 4.2
#3 2.0 4.0
#4 1.5 1.5
推荐阅读
- networking - 无法从单独的 LAN 机器访问 Kafka Broker
- apache - 如何使用 Apache/httpd.conf 为整个域重定向 410(在一行中)?
- python - 在Python中,如果该元素存在于文件中,我如何打印列表的元素
- javascript - 如何在 axios-http 中发布对象数组
- json - 在 github 操作中创建 JSON 字符串时出错
- json - Spotify API 如何匹配非英语数据(例如艺术家姓名)?
- scala - Apache Spark - foreachRDD 中的异常处理
- javascript - 如何在变量打字稿中动态设置当前时间
- json - 使用 json.load 时出现 json.decoder.JSONDecodeError 错误
- javascript - 函数未在 Jest 测试中调用模拟构造函数