R中数据框的空间子集

Posted 2023-02-14

【中文标题】R中数据框的空间子集

【英文标题】：Spatial subset of data frame in R

【发布时间】：2022-01-09 04:06:42

【问题描述】：

我有一个使用plot() 绘制的大数据框。然后我用了：

library(splancs) 

polygon_xy = getpoly(quiet=FALSE) 

并在绘图上绘制点以选择我感兴趣的区域。这生成了我绘制的多边形的 x,y 坐标。

我想提取位于多边形内的数据，或将我的 df 子集化以仅包含位于多边形内的点。有什么建议吗？

【问题讨论】：

你从polygon_xy得到什么样的输出？它与点在同一坐标系中吗？还是相对于情节（即0.3205943、0.9671482）？

这就是它的样子： > polygon_xy [,1] [,2] [1,] 2.483359 0.7856676 [2,] 1.773868 1.1975744 [3,] 1.773868 2.5936972 [4,] 2.714676 3.914031 5,] 3.713715 2.2312045 [6,] 3.652968 1.2096540

我觉得和剧情有关，但不知道如何检验这个假设

我提供了几个提取选项。似乎主要问题是点坐标的格式是否不正确。这些点需要位于具有 x 和 y 列以及坐标的矩阵中。

【参考方案1】：

当您使用plot 和getpoly 识别坐标时，点数据需要采用特定格式（即带有x 和y 的矩阵）。

library(splancs)
library(tidyverse)
library(sf)

set.seed(543)
xy <-
  cbind(x = runif(n = 25, min = -118, max = -117),
        y = runif(n = 25, min = 40, max = 42))

plot(xy)

# Draw a polygon for study area.
poly <- getpoly()

# Convert to sf objects.
polysf <- st_as_sf(as.data.frame(poly), coords = c("V1", "V2"), crs = 4326) %>% 
  dplyr::summarise() %>%
  st_cast("POLYGON") %>% 
  st_convex_hull()

xysf <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = 4326)

# Do an intersection to keep only points inside the drawn polygon.
xy_intersect <- st_intersection(polysf, xysf)

输出

Simple feature collection with 9 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -117.7913 ymin: 40.82405 xmax: -117.4264 ymax: 41.7448
Geodetic CRS:  WGS 84
                    geometry
1 POINT (-117.4264 41.18712)
2  POINT (-117.5756 41.7448)
3 POINT (-117.7913 40.82405)
4 POINT (-117.7032 41.15077)
5 POINT (-117.5634 41.23936)
6 POINT (-117.7441 40.84163)
7  POINT (-117.692 41.27514)
8 POINT (-117.6864 40.98462)
9 POINT (-117.5759 40.88477)

使用来自library(mapview) 的mapview::mapview(xy_intersect) 绘制

但是，如果您想从原始数据框中提取行，那么这里有另一个技巧，用于提取落在绘制多边形内的点（例如，当多边形坐标看起来像 0.003456 时）。

library(splancs)
library(tidyverse)

set.seed(543)
xy <-
  cbind(x = runif(n = 25, min = -118, max = -117),
        y = runif(n = 25, min = 40, max = 42))

plot(xy)

# Draw a polygon for study area.
poly <- getpoly()

# Plot the results.
plot(xy)
polygon(poly)

# This will return a logical vector for points in the polygon
io <- inout(xy, poly)
points(xy[io,], pch = 16, col = "blue")

# Then, can use the index from io to extract the points that 
# are inside the polygon from the original set of points.
extract_points <- as.data.frame(xy)[which(io == TRUE),]

extract_points

输出

           x        y
2  -117.4506 41.17794
3  -117.4829 40.71030
8  -117.4679 40.71702
19 -117.3354 40.53687
21 -117.5219 40.47077
22 -117.4876 40.18188
25 -117.2015 40.86243

【讨论】：

非常感谢！我使用了第二个选项，它奏效了。 inout() 正是我需要的魔法！