Colheita para SpatialPolygonsDataFrame

Eu tenho dois arquivos 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 

Eu quero cortar o arquivo dat1 usando a extensão de dat2. Eu não sei como fazer isso. Eu apenas manipulo arquivos raster usando a function “crop” antes.

Quando eu uso essa function para meus dados atuais, ocorre o seguinte erro:

 > r1 <- crop(BiomassCarbon.shp,alberta.shp) Error in function (classes, fdef, mtable) : unable to find an inherited method for function 'crop' for signature"SpatialPolygonsDataFrame"' 

Método simplificado adicionado em 2014-10-9 :

raster::crop() pode ser usado para cortar objects Spatial* (bem como Raster* ).

Por exemplo, aqui está como você pode usá-lo para cortar um object 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") 

Resposta original (e ainda funcional):

A function gIntersection() torna isso bastante simples.

Usando o exemplo bacana do mnel como um ponto 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") 

insira a descrição da imagem aqui

Aqui está um exemplo de como fazer isso com rgeos usando o mapa do mundo como um exemplo

Isto vem de Roger Bivand na lista de discussão R-sig-Geo . Roger é um dos autores do pacote sp .

Usando o mapa do mundo como um exemplo

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

insira a descrição da imagem aqui

Você não pode usar o crop em objects de polígonos sp. Você precisará criar um polígono que represente as coordenadas bbox de dat2 e, em seguida, usar gIntersects na biblioteca rgeos.

Edit: Este comentário foi em relação à versão disponível em 2012 e isso não é mais o caso.

ver?

corp (x, y, nome do arquivo = “”, snap = ‘near’, datatype = NULL, …)

Objeto raster * x

y Objeto de extensão ou qualquer object a partir do qual um object de extensão pode ser extraído (consulte Detalhes

Você precisa rasterizar o primeiro SpatialPolygon, usando a function rasterize do pacote de rasterização

Eu crio alguns dados para mostrar como usar o rasterize:

 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 

se eu tentar

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

Agora usando rasterize

 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) 
    Intereting Posts