Как отфильтровать полигоны sfc по значениям long / lat?

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)

(Очевидно, это работает, если оба значения polys и «степень фильтрации» имеют значения lat / long, в противном случае Вы должны позаботиться о перепроектировании, и так далее.) ^{Создано в 2018-11-15 пакетом Представ (v0.2.1)}