wake-up-neo.net

Zuschneiden für SpatialPolygonsDataFrame

Ich habe zwei SpatialPolygonsDataFrame-Dateien: 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

Ich möchte die Datei dat1 mit dem Umfang von dat2 beschneiden. Ich weiß nicht, wie es geht. Ich verarbeite nur Rasterdateien mit der Funktion "Zuschneiden".

Wenn ich diese Funktion für meine aktuellen Daten verwende, tritt der folgende Fehler auf:

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

 unable to find an inherited method for function ‘crop’ for signature"SpatialPolygonsDataFrame"’
rcropspatial
25
Jian Zhang

Optimierte Methode hinzugefügt 2014-10-9 :

raster::crop() kann verwendet werden, um Spatial* (sowie Raster*) Objekte zu beschneiden.

So können Sie beispielsweise ein SpatialPolygons*-Objekt beschneiden:

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

Ursprüngliche (und immer noch funktionale) Antwort:

Die rgeos function gIntersection() macht dies ziemlich unkompliziert.

Verwenden Sie das clevere Beispiel von mnel als Absprungpunkt:

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

53
Josh O'Brien

Hier ein Beispiel, wie Sie dies mit rgeos am Beispiel der Weltkarte tun

Dies kommt von Roger Bivand in der R-sig-Geo-Mailingliste . Roger ist einer der Autoren des Pakets sp. 

Am Beispiel der Weltkarte

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

5
mnel

Sie können keinen Zuschnitt für sp Polygonobjekte verwenden. Sie müssen ein Polygon erstellen, das die Bbox-Koordinaten von dat2 darstellt. Anschließend können Sie gIntersects in der Rgeos-Bibliothek verwenden.

Bearbeiten: Dieser Kommentar bezieht sich auf die 2012 verfügbare Version und ist nicht mehr der Fall. 

1
Jeffrey Evans

sehen?

corp (x, y, dateiname = "", snap = 'near', Datentyp = NULL, ...)

x Raster * Objekt

y Ausdehnungsobjekt oder ein beliebiges Objekt, aus dem ein Ausdehnungsobjekt bestehen kann extrahiert (siehe Details

Sie müssen das erste SpatialPolygon mit der Funktion rasterize aus dem Rasterpaket rastern

Ich erstelle Daten, um zu zeigen, wie man rasterize verwendet:

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

wenn ich es versuche:

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

Jetzt rastern

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)
0
agstudy