Обрезать простые объекты объекта в R

20

Есть ли функция обрезки sf карты, похожая на maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))используемую для SpatialPolygon или SpatialLine?

Я рассматриваю, st_intersection()но может быть правильный путь.

Казухито
источник

Ответы:

17

st_intersectionэто, наверное, лучший способ. Найдите любой способ, который лучше всего подходит sfдля пересечения объекта с вашим входом. Вот способ, использующий удобство raster::extentи сочетание старого и нового. ncсоздается example(st_read):

st_intersection(nc, st_set_crs(st_as_sf(as(raster::extent(-82, -80, 35, 36), "SpatialPolygons")), st_crs(nc)))

Я не думаю, что вы можете уговорить st_intersectionне требовать точного соответствия CRS, поэтому либо установите оба параметра на NA, либо убедитесь, что они одинаковые. Для bbox / экстента afaik нет простых инструментов, поэтому использование растра - хороший способ облегчить задачу.

mdsumner
источник
Большое спасибо @mdsumner, это работает как шарм. Я провел часы, st_intersectionно не мог решить это сам.
Казухито
Теперь вы можете использовать spex::spexдля замены st_as_sf(as(...))вызова. Кроме того, tmaptools::crop_shape()можете сделать это.
AF7
1
sfтеперь включает st_crop, см. мой ответ для деталей.
AF7
23

С сегодняшнего дняst_crop в github-версии есть функция sf( devtools::install_github("r-spatial/sf")возможно, на CRAN и в ближайшем будущем).

Просто выпустите:

st_crop(nc, c(xmin=-82, xmax=-80, ymin=35, ymax=36))

Вектор должен быть назван с xmin xmax ymin ymax(в любом порядке).

Вы также можете использовать любой объект, который может быть прочитан st_bboxкак предел обрезки, что очень удобно.

AF7
источник
5

Другой обходной путь, для меня это было быстрее для больших шейп-файлов:

library(sf)
library(raster)
library(rgeos)
library(ggplot2)

# Load National Forest shapefile
# https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.AdministrativeForest.zip
nf.poly <- st_read("S_USA.AdministrativeForest"), "S_USA.AdministrativeForest")

crop_custom <- function(poly.sf) {
  poly.sp <- as(poly.sf, "Spatial")
  poly.sp.crop <- crop(poly.sp, extent(c(-82, -80, 35, 36)))
  st_as_sf(poly.sp.crop)
}

cropped <- crop_custom(nf.poly)
pbaylis
источник
Благодарю. Это интересный рабочий процесс, комбинация raster :: crop () и st_as_sf () ... +1 от меня. Хотелось бы, чтобы в будущих версиях sf у нас была такая легкодоступная функция, как crop () . Что касается скорости, быстрый запуск system.time с вашей функцией сообщил о пользователе: 5.42, система: 0.09, прошло 5.52 , а st_intersection()подход был пользователь: 1.18, система: 0.05, прошло 1.23 в вашем наборе данных. (Возможно , мое окружение отличается с вашим ... не уверен.)
Казухито
Это интересно - подход st_intersection занимает у меня около 80-х годов.
pbaylis
Имейте в виду, что функция raster :: crop применительно к объектам sp geometry выполняет функцию оболочки для функций rgeos. Хотя и очень удобная обертка. GEOS API работает с объектами WKT, поэтому он всегда будет стандартом для операций наложения sf.
Джеффри Эванс
1
И со временем оно меняется, теперь в sf встроена «пространственная индексация» в 0.5-1 cran.r-project.org/web/packages/sf/news.html, поэтому, вероятно, теперь она быстрее, чем sp / rgeos.
mdsumner
1
sfтеперь включает st_crop, см. мой ответ для деталей.
AF7
1

Решение @ mdsumner как функция. Работает, если rastaесть RasterBrick, экстент, bbox и т. Д.

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf = function(sfdf, rasta) {
  st_intersection(sfdf, st_set_crs(st_as_sf(as(extent(rasta), "SpatialPolygons")), st_crs(sfdf)))
}

Он выбрасывает информацию о растровом crs, потому что я не знаю, как преобразовать растр crs () в st_crs ()

На моем компьютере и для моего образца данных это эквивалентно производительности raster::cropс версией данных SpatialLinesDataFrame.

Решение @ pbaylis примерно в 2,5 раза медленнее:

# Slower option.
crop.sf2 = function(sfdf, rasta) {
  st_as_sf(crop(as(sfdf, "Spatial"), rasta))
}

Редактировать: Комментарий Somebody предлагает spex , который производит SpatialPolygons с crs из раста, если у него есть crs.

Этот код использует тот же метод, что и spex:

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf3 <- function(sfdf, rasta) {
  # Get extent and crs
  ext.sp <- as(extent(rasta), "SpatialPolygons")
  crs(ext.sp) <- crs(rasta)

  # crop
  st_intersection(sfdf, st_as_sf(ext.sp))
}
КМЦ
источник
sf теперь имеет st_cropфункцию, которую, вероятно, стоит проверить.
CMC