Recortar objeto de características simples en R

17

¿Existe una función para recortar el objeto del mapa sf, similar al maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15)) utilizado para SpatialPolygon o SpatialLine?

Estoy considerando st_intersection() pero puede haber una forma adecuada.

    
pregunta Kazuhito 06.03.2017 - 09:01

4 respuestas

14

st_intersection es probablemente la mejor manera. Encuentre la forma que mejor funcione para obtener un objeto sf que se interseca con su entrada. Aquí hay una manera de usar la conveniencia de raster::extent y una mezcla de lo antiguo y lo nuevo. nc es creado por example(st_read) :

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

No creo que puedas convencer a st_intersection para que no necesites un CRS exacto, así que configura ambos en NA o asegúrate de que sean iguales. No hay herramientas fáciles para bbox / extension afaik, por lo que usar raster es una buena manera de hacer las cosas más fáciles.

    
respondido por el mdsumner 06.03.2017 - 12:22
17

Desde hoy , hay una función st_crop en la versión de github de sf ( devtools::install_github("r-spatial/sf") , probablemente también en CRAN en un futuro cercano).

Solo problema:

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

El vector debe tener un nombre con xmin xmax ymin ymax (en el orden que sea).

También puedes usar cualquier objeto que se pueda leer con st_bbox como límites de recorte, lo que es muy útil.

    
respondido por el AF7 23.04.2018 - 09:33
4

Otra solución alternativa, para mí fue más rápido para archivos de formas más grandes:

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)
    
respondido por el pbaylis 01.04.2017 - 08:40
0

La solución de @ mdsumner como una función. Funciona si rasta es un RasterBrick, extension, bbox, etc.

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

Tira la información crs de la trama porque no sé cómo convertir una raster crs () en st_crs ()

En mi máquina y para mi muestra de datos, esto tiene un rendimiento equivalente a raster::crop con una versión de SpatialLinesDataFrame de los datos.

La solución de

@ pbaylis es aproximadamente 2.5 veces más lenta:

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

Edición: el comentario de Somebodies sugiere spex , que produce SpatialPolygons con los crs desde el rasta, si tiene un crs.

Este código utiliza el mismo método que 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))
}
    
respondido por el cmc 11.04.2018 - 17:18

Lea otras preguntas en las etiquetas