cortar - Recorte para SpatialPolygonsDataFrame

cortar raster en r

Question

Tengo dos archivos SpatialPolygonsDataFrame : dat1, dat2

extent(dat1)
class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 


extent(dat2)
class       : Extent 
xmin        : -120.0014 
xmax        : -109.9997 
ymin        : 48.99944 
ymax        : 60 

Quiero recortar el archivo dat1 utilizando la extensión de dat2. No se como hacerlo Yo solo manejo archivos raster usando la función "recortar" antes.

Cuando uso esta función para mis datos actuales, ocurre el siguiente error:

> r1 <- crop(BiomassCarbon.shp,alberta.shp)
Error in function (classes, fdef, mtable)  : 

 unable to find an inherited method for function ‘crop’ for signature"SpatialPolygonsDataFrame"’

Answer 1

Aquí hay un ejemplo de cómo hacer esto con rgeos usando el mapa del mundo como ejemplo

Esto viene de Roger Bivand en la lista de correo R-sig-Geo . Roger es uno de los autores del paquete sp .

Usando el mapa del mundo como ejemplo

library(maptools)

data(wrld_simpl)

# interested in the arealy bounded by the following rectangle
# rect(130, 40, 180, 70)

library(rgeos)
# create  a polygon that defines the boundary
bnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))
# convert to a spatial polygons object with the same CRS
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),
proj4string=CRS(proj4string(wrld_simpl)))
# find the intersection with the original SPDF
gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
# create the new spatial polygons object.
out <- vector(mode="list", length=length(which(gI)))
ii <- 1
for (i in seq(along=gI)) if (gI[i]) {
  out[[ii]] <- gIntersection(wrld_simpl[i,], SP)
  row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1
}
# use rbind.SpatialPolygons method to combine into a new object.
out1 <- do.call("rbind", out)
# look here is Eastern Russia and a bit of Japan and China.
plot(out1, col = "khaki", bg = "azure2")

Answer 2

No se puede usar el recorte en objetos de polígono sp. Necesitará crear un polígono que represente las coordenadas bbox de dat2 y luego puede usar gIntersects en la biblioteca rgeos.

Edición: este comentario fue en relación con la versión disponible en 2012 y ya no es así.

Answer 3

ver? cosechar

corp (x, y, filename = "", snap = ''near'', datatype = NULL, ...)

x Raster * objeto

y Objeto de extensión, o cualquier objeto del que se pueda extraer un objeto de Extensión (ver Detalles

Debe rasterizar el primer SpatialPolygon, usando la función rasterize del paquete raster

Creo algunos datos para mostrar cómo usar rasterizar:

n <- 1000
x <- runif(n) * 360 - 180
y <- runif(n) * 180 - 90
xy <- cbind(x, y)
vals <- 1:n
p1 <- data.frame(xy, name=vals)
p2 <- data.frame(xy, name=vals)
coordinates(p1) <- ~x+y
coordinates(p2) <- ~x+y

si lo intento

 crop(p1,p2)
 unable to find an inherited method for function ‘crop’ for signature ‘"SpatialPointsDataFrame"’

Ahora usando rasterizar

r <- rasterize(p1, r, ''name'', fun=min)
crop(r,p2)

class       : RasterLayer 
dimensions  : 18, 36, 648  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : 1, 997  (min, max)

Answer 4

Método simplificado añadido 2014-10-9 :

raster::crop() se puede usar para recortar objetos Spatial* (así como Raster* ).

Por ejemplo, aquí es cómo podría usarlo para recortar un objeto SpatialPolygons* :

## Load raster package and an example SpatialPolygonsDataFrame
library(raster) 
data("wrld_simpl", package="maptools")

## Crop to the desired extent, then plot
out <- crop(wrld_simpl, extent(130, 180, 40, 70))
plot(out, col="khaki", bg="azure2")

Respuesta original (y aún funcional):

La función gIntersection() hace que esto sea bastante sencillo.

Usando el ingenioso ejemplo de mnel como punto de partida:

library(maptools)
library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
library(rgeos)
data(wrld_simpl)

## Create the clipping polygon
CP <- as(extent(130, 180, 40, 70), "SpatialPolygons")
proj4string(CP) <- CRS(proj4string(wrld_simpl))

## Clip the map
out <- gIntersection(wrld_simpl, CP, byid=TRUE)

## Plot the output
plot(out, col="khaki", bg="azure2")