2012-12-21 13 views
24

Ho due SpatialPolygonsDataFrame file: DAT1, DAT2Crop per SpatialPolygonsDataFrame

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 

voglio ritagliare il DAT1 file utilizzando l'estensione di DAT2. Non so come farlo Ho appena gestito i file raster utilizzando la funzione "ritaglio" prima.

Quando uso questa funzione per i miei dati attuali, verifica il seguente errore:

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

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

questo è senza speranza, il lavoro sui dettagli per una domanda che coinvolge DAT1 e DAT2, o di quelle altre cose – mdsumner

risposta

48

metodo semplificato aggiunto 2014/10/09:

raster::crop() può essere utilizzato per ritagliare Spatial* (così come Raster*) oggetti.

Ad esempio, ecco come si potrebbe usarlo per ritagliare un oggetto 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") 

originale (ed ancora funzionale) risposta:

I rgeos funzione gIntersection() marche questo piuttosto semplice.

Utilizzando esempio elegante di mnel come un punto di partenza:

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

enter image description here

+0

Molto più semplice. – mnel

+1

Grazie per aver fornito un esempio di codice. Non ho avuto tempo quando ho postato la mia risposta. +1 per la tua soluzione intelligente nell'utilizzo del pacchetto raster per creare il poligono di delimitazione. –

+1

@JeffreyEvans - Sì. Trovo che ** raster ** svolga un * molto * migliore lavoro nel fornire agli utenti metodi e metodi di conversione facili da usare rispetto a ** sp **. Apprezzo che gli autori ** sp ** vogliano mantenere il loro pacco di ricambio e leggero (dal momento che è inteso a sopportare così tanti altri pacchetti) ma spero ancora che alla fine aggiungano qualche altra funzione di utilità. –

0

vedono raccolto

corp (x, y, filename = "", scatto = 'vicino? ', datatype = NULL, ...)

x raster * oggetto

oggetto y Estensione, o qualsiasi oggetto da cui un oggetto punto può essere estratto (vedi dettagli

È necessario rasterizzare il primo SpatialPolygon, utilizzando rasterize funzione dal pacchetto raster

creo alcuni dati per mostrare come utilizzare 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 provo:

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

Ora utilizzando 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) 
1

Non è possibile utilizzare delle colture su oggetti poligonali sp. Dovrai creare un poligono che rappresenti le coordinate bbox di dat2 e quindi utilizzare gIntersects nella libreria rgeos.

Modifica: questo commento era in relazione alla versione disponibile nel 2012 e non è più il caso.

5

Ecco un esempio di come fare questo con rgeos utilizzando la mappa del mondo come un esempio

Questo viene da Roger Bivand su R-sig-Geo mailing list. Roger è uno degli autori del pacchetto sp.

Utilizzando la mappa del mondo come un esempio

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

enter image description here

+0

risposta Cool. 'gIntersection()' lo rende ancora più semplice di 'gIntersects()'. –

Problemi correlati