r - R - SF(同心圆多边形)中两个多边形之间的过滤/子集数据

Question

不确定是否缺少某个功能，但我在过滤/检查点几何图形是否落在两个多边形(同心圆)之间时遇到问题。

您可以在两个同心多边形之间创建一个掩模，然后使用它来过滤掉包含感兴趣的要素数据的点几何图形吗？

我尝试使用 R 中的 sf_filter 包在两个多边形之间进行子集化。这不起作用。

可重现的代码如下:

library(sf)
library(tidyverse)
library(sp)

#Create fake data
my.data <- data.frame(replicate(2,sample(-88: -14,100,rep=TRUE))) # Point data
d <- cbind(seq(-180,180,length.out=360),rep(-88,360))
e <- cbind(seq(-180,180,length.out=360),rep(-30,360))

#Project fake data
d = SpatialPoints(cbind(d[,1], d[,2]), proj4string=CRS("+proj=longlat"))
d <- spTransform(d, CRS("+init=epsg:3976"))

e = SpatialPoints(cbind(e[,1], e[,2]), proj4string=CRS("+proj=longlat"))
e <- spTransform(e, CRS("+init=epsg:3976"))

my.data = SpatialPoints(cbind(my.data[,1], my.data[,2]), proj4string=CRS("+proj=longlat"))
my.data <- spTransform(my.data, CRS("+init=epsg:3976"))

d <- sf::st_as_sf(d, coords = c("X1", "X2"),
                  remove = FALSE,
                  crs = st_crs("epsg:3976"))

e <- sf::st_as_sf(e, coords = c("X1", "X2"),
                  remove = FALSE,
                  crs = st_crs("epsg:3976"))

my.data <- sf::st_as_sf(my.data, coords = c("X1", "X2"),
                        remove = FALSE,
                        crs = st_crs("epsg:3976"))
# Create linestrings from circle
d <- d %>% 
  summarise(do_union = FALSE) %>%
  st_cast("LINESTRING")

e <- e %>% 
  summarise(do_union = FALSE) %>%
  st_cast("LINESTRING")

#Join geometries
nst <- rbind(d,e) 

#Create polygon
nst <- nst %>% 
  st_cast("POLYGON")

#Filtering between polygons doesn't return anything
PFz <- st_filter(my.data,nst)

Answer 1

考虑这种方法；它建立在北卡罗来纳州的三个半随机城市上(因为我喜欢随 {sf} 一起提供的 nc.shp)

它的作用是构建两个缓冲区作为 sf 对象，然后构建两个逻辑向量 - 用于大圆和小圆的 sf::st_contains()。

然后这是一个检查点的简单逻辑操作:

1. 都包含在大圆圈内，并且同时
2. 不包含在小圆圈内

如果你想变得更奇特，你可以运行 sf::st_difference()在两个缓冲区对象上，直接获取掩码并仅对“rim”对象检查一次 sf::st_contains() 。

library(sf)
library(dplyr)

# 3 semi rancom cities in NC (because I *deeply love* the nc.shp file)
cities <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
                     x = c(-78.633333, -79.819444, -77.912222),
                     y = c(35.766667, 36.08, 34.223333)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326)

# small buffer - Greensboro will be in; Wilmington not
small_buffer <- cities %>% 
  filter(name == "Raleigh") %>% 
  st_geometry() %>% 
  st_buffer(units::as_units(100, "mile"))

# big buffer - both Greensboro & Wilmington are  in
big_buffer <- cities %>% 
  filter(name == "Raleigh") %>% 
  st_geometry() %>% 
  st_buffer(units::as_units(150, "mile"))

# a visual overview
mapview::mapview(list(big_buffer, small_buffer, cities))

# vector of cities in big buffer
in_big_buffer <- st_contains(big_buffer,
                             cities,
                             sparse = F) %>% 
  t()

# vector of cities in small buffer
in_small_buffer <- st_contains(small_buffer,
                               cities,
                               sparse = F) %>% 
  t()

# cities in concentric circle = in big, and not in small
in_concentric_circle <- in_big_buffer & !in_small_buffer

# check - subset of cities by logical vector
cities %>% 
  filter(in_concentric_circle)

# Simple feature collection with 1 feature and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.91222 ymin: 34.22333 xmax: -77.91222 ymax: 34.22333
# Geodetic CRS:  WGS 84
#         name                   geometry
# 1 Wilmington POINT (-77.91222 34.22333)