Przytnij obiekt prostych funkcji w R.

20

Czy istnieje funkcja przycinania obiektu mapy sf, podobna do tej maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))używanej w przypadku SpatialPolygon lub SpatialLine?

Rozważam, st_intersection()ale może być właściwy sposób.

Kazuhito
źródło

Odpowiedzi:

17

st_intersectionjest prawdopodobnie najlepszym sposobem. Znajdź sposób, który działa najlepiej, aby sfobiekt przecinał się z wprowadzonymi danymi. Oto sposób korzystania z wygody raster::extenti połączenia starego i nowego. ncjest tworzony przez example(st_read):

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

Nie sądzę, abyś mógł się przekonać, st_intersectionaby nie potrzebować dokładnie pasującego CRS, więc albo ustaw je na NA, albo upewnij się, że są takie same. Nie ma łatwych narzędzi do bbox / scope afaik, więc użycie rastra jest dobrym sposobem na uproszczenie.

mdsumner
źródło
Dzięki bardzo @mdsumner, działało to jak urok. Spędziłem godziny, st_intersectionale sam nie mogłem tego rozwiązać.
Kazuhito,
Możesz teraz użyć, spex::spexaby zastąpić st_as_sf(as(...))połączenie. Również tmaptools::crop_shape()mogę to zrobić.
AF7
1
sfzawiera teraz st_crop, zobacz moją odpowiedź, aby uzyskać szczegółowe informacje.
AF7,
23

Od dzisiaj istnieje st_cropfunkcja w wersji github sf( devtools::install_github("r-spatial/sf")prawdopodobnie na CRAN również w najbliższej przyszłości).

Wystarczy wydać:

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

Wektor musi być nazwany za pomocą xmin xmax ymin ymax(w dowolnej kolejności).

Możesz także użyć dowolnego obiektu, który można odczytać, st_bboxjako ograniczenia kadrowania, co jest bardzo przydatne.

AF7
źródło
5

Kolejne obejście, dla mnie było szybsze dla większych plików kształtów:

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
źródło
Dzięki. To ciekawy przepływ pracy, połączenie raster :: crop () i st_as_sf () ... + 1 ode mnie. Chciałbym, abyśmy mieli taką łatwo dostępną funkcję, jak crop () w przyszłych wersjach sf . Jeśli chodzi o szybkość, szybki bieg system.time z twoją funkcją zgłosił użytkownika: 5,42, system: 0,09, minęło 5,52 , podczas gdy st_intersection()podejście było użytkownikiem: 1,18, system: 0,05, minęło 1,23 w zbiorze danych. (Prawdopodobnie moje środowisko różni się od twojego ... nie jestem pewien.)
Kazuhito,
To ciekawe - podejście st_intersection zajmuje mi około 80 lat.
pbaylis,
Należy pamiętać, że funkcja raster :: crop, zastosowana do obiektów geometrii sp, pełni funkcję otoki dla funkcji rgeos. Chociaż bardzo wygodne opakowanie. Interfejs API GEOS działa na obiektach WKT, więc niezmiennie będzie standardem dla operacji nakładania sf.
Jeffrey Evans
1
I zmienia się z czasem, sf ma teraz wbudowane „indeksowanie przestrzenne” w wersji 0.5-1 cran.r-project.org/web/packages/sf/news.html, więc prawdopodobnie jest teraz szybszy niż sp / rgeos.
mdsumner
1
sfzawiera teraz st_crop, zobacz moją odpowiedź, aby uzyskać szczegółowe informacje.
AF7,
1

@ Rozwiązanie mdsumner jako funkcja. Działa, jeśli rastajest to RasterBrick, zakres, bbox itp.

# 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)))
}

Wyrzuca informacje CRS o rastrze, ponieważ nie wiem jak przekonwertować CRS raster () na st_crs ()

Na moim komputerze i dla mojej próbki danych ma to wydajność równoważną raster::cropz wersją danych SpatialLinesDataFrame.

Rozwiązanie @ pbaylis jest około 2,5 razy wolniejsze:

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

Edycja: Komentarz Somebodies sugeruje spex , który tworzy SpatialPolygons z crs z rasta, jeśli ma crs.

Ten kod używa tej samej metody co 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))
}
cmc
źródło
sf ma teraz st_cropfunkcję, którą prawdopodobnie warto sprawdzić.
cmc