如何通过长/纬度值过滤 sfc 多边形？

Question

假设t是：

t <- structure(list(structure(list(structure(c(-89.990791, -89.990772, 
-89.990901, -89.99092, -89.990791, 30.727025, 30.727083, 30.727114, 
30.727057, 30.727025), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(-89.991691, -89.991755, -89.991755, 
-89.991691, -89.991691, 30.716004, 30.716004, 30.715916, 30.715916, 
30.716004), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"
))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = -89.991755, 
ymin = 30.715916, xmax = -89.990772, ymax = 30.727114), class = "bbox"), crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)
> t
Geometry set for 2 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -89.99176 ymin: 30.71592 xmax: -89.99077 ymax: 30.72711
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((-89.99079 30.72703, -89.99077 30.7270...
POLYGON ((-89.99169 30.716, -89.99175 30.716, -...

如何按长/纬度边界过滤多边形？假设我想过滤掉任何 hax lat>30.72 的多边形（因此只保留第二个多边形）。有什么特定的函数可以用来过滤多边形吗？

Answer 1

不知道有没有现成的函数，但是一个长方形的“空间过滤器”很容易搭建。只需定义“角”，从它们创建一个 bbox，转换为多边形并找到包含/重叠“过滤器区域”的原始多边形。

这是一个快速而肮脏的例子：

library(sf)
polys_sf<-st_read(system.file("shape/nc.shp", package="sf") ) %>% 
  st_transform(crs="+init=epsg:4326")
plot(st_geometry(polys_sf))

定义一个“空间过滤器”

xmin <- -80
xmax <- -76
ymin <- 34
ymax <- 36

根据过滤器创建一个多边形。（您可以对某些值使用“NA”，例如，如果您只想过滤“左侧”，则可以将 xmax 设置为 NA）

filt_bbox <- sf::st_bbox(c(xmin = ifelse(is.na(xmin), -180, xmin), 
                           ymin = ifelse(is.na(ymin),  -90,  ymin), 
                           xmax = ifelse(is.na(xmax), +180, xmax), 
                           ymax = ifelse(is.na(ymax),  +90, ymax)), 
                         crs = st_crs(4326)) %>% 
  sf::st_as_sfc(.)

现在根据 bbox 多边形“过滤”原始数据集：st_within如果您只想保留完全包含在定义区域中的多边形，请使用

find_data <- sf::st_within(polys_sf, filt_bbox)
#> although coordinates are longitude/latitude, st_within assumes that they are planar
filt_data <- polys_sf[which(lengths(find_data) != 0), ]

plot(filt_bbox)
plot(st_geometry(filt_data), add = TRUE, reset = FALSE)

st_intersects如果要保留与定义区域相交的所有多边形，请使用

find_data <- sf::st_intersects(polys_sf, filt_bbox)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
filt_data <- polys_sf[which(lengths(find_data) != 0), ]

plot(st_geometry(filt_data))
plot(filt_bbox, add = TRUE)

（显然，如果您的多边形和“过滤范围”都是纬度/经度，则此方法有效，否则您必须注意重新投影等。） ^{由reprex 包(v0.2.1)于 2018 年 11 月 15 日创建}