Objekt für einfache Features in R zuschneiden

20

Gibt es eine Funktion zum Zuschneiden von SF-Kartenobjekten, ähnlich wie maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))sie für SpatialPolygon oder SpatialLine verwendet werden?

Ich überlege, st_intersection()aber es kann einen richtigen Weg geben.

Kazuhito
quelle

Antworten:

17

st_intersectionist wohl der beste weg. Finden Sie heraus, wie es am besten funktioniert, um ein sfObjekt mit Ihrer Eingabe zu kreuzen. Hier finden Sie eine Möglichkeit, die Bequemlichkeit von raster::extentund eine Mischung aus Altem und Neuem zu nutzen. ncwird erstellt von example(st_read):

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

Ich glaube nicht, dass Sie dazu überreden können st_intersection, ein genau passendes CRS nicht zu benötigen. Setzen Sie daher entweder beide auf NA oder stellen Sie sicher, dass sie gleich sind. Es gibt keine einfachen Tools für Bbox / Extent Afaik, daher ist die Verwendung von Raster eine gute Möglichkeit, die Dinge einfacher zu gestalten.

mdsumner
quelle
Vielen Dank @mdsumner, es hat wie ein Zauber funktioniert. Ich habe Stunden damit verbracht st_intersection, konnte es aber nicht selbst lösen.
Kazuhito
Sie können jetzt verwenden spex::spex, um den st_as_sf(as(...))Anruf zu ersetzen . Auch tmaptools::crop_shape()kann dies tun.
AF7
1
sfjetzt enthält st_crop, siehe meine Antwort für Details.
AF7
22

Seit heute gibt es eine st_cropFunktion in der Github-Version von sf( devtools::install_github("r-spatial/sf")wahrscheinlich auch auf CRAN in naher Zukunft).

Einfach ausgeben:

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

Der Vektor muss mit xmin xmax ymin ymax(in welcher Reihenfolge auch immer) benannt werden.

Sie können auch jedes Objekt verwenden, das st_bboxals Beschneidungsgrenze gelesen werden kann , was sehr praktisch ist.

AF7
quelle
5

Eine andere Problemumgehung, für mich war es bei größeren Shapefiles schneller:

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
quelle
Vielen Dank. Es ist ein interessanter Workflow, eine Kombination aus raster :: crop () und st_as_sf () ... + 1 von mir. Ich wünschte, wir könnten eine solche leicht zugängliche Funktion wie crop () in zukünftigen Versionen von sf haben . Was die Geschwindigkeit anbelangt, so meldete ein schneller Durchlauf von system.time mit Ihrer Funktion Benutzer: 5.42, System: 0.09, verstrichene 5.52 , während st_intersection()Annäherung war Benutzer: 1.18, System: 0.05, verstrichene 1.23 in Ihrem Datensatz war. (Wahrscheinlich ist meine Umgebung bei Ihnen anders ... nicht sicher.)
Kazuhito
Das ist interessant - für mich dauert der Ansatz von st_intersection ungefähr 80s.
Pbaylis
Beachten Sie, dass die Funktion raster :: crop bei der Anwendung auf sp-Geometrieobjekte einen Wrapper für Rgeos-Funktionen darstellt. Ein sehr praktischer Wrapper. Die GEOS-API verarbeitet WKT-Objekte und wird daher immer ein Standard für SF-Overlay-Operationen sein.
Jeffrey Evans
1
Und es ändert sich mit der Zeit, sf hat jetzt eine "räumliche Indizierung" in 0.5-1 eingebaut. Cran.r-project.org/web/packages/sf/news.html ist also wahrscheinlich schneller als sp / rgeos.
mdsumner
1
sfjetzt enthält st_crop, siehe meine Antwort für Details.
AF7
1

@ mdsumner's Lösung als Funktion. Funktioniert, wenn rastaes sich um einen RasterBrick, einen Extent, eine Bbox usw. handelt.

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

Es wirft die crs-Informationen des Rasters weg, weil ich nicht weiß, wie man ein Raster crs () in ein st_crs () konvertiert

Auf meinem Computer und für mein raster::cropDatenbeispiel entspricht dies der Leistung einer SpatialLinesDataFrame-Version der Daten.

Die Lösung von @pbaylis ist etwa 2,5-mal langsamer:

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

Edit: Ein Kommentar von jemandem schlägt Spex vor , das SpatialPolygons mit dem Crs aus dem Rasta erzeugt, wenn es ein Crs hat.

Dieser Code verwendet dieselbe Methode wie 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
quelle
sf hat jetzt eine st_cropfunktion, die es wahrscheinlich wert ist, ausgecheckt zu werden.
23.