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.

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.


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.


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)

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

sf ma teraz st_cropfunkcję, którą prawdopodobnie warto sprawdzić.
— cmc
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.