cắt xén để SpatialPolygonsDataFrame

Question

Tôi có hai SpatialPolygonsDataFrame file: 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 

Tôi muốn cắt dat1 tập tin bằng cách sử dụng mức độ dat2. Tôi không biết phải làm thế nào. Tôi chỉ xử lý các tập tin raster bằng cách sử dụng chức năng "cắt" trước đây.

Khi tôi sử dụng chức năng này cho dữ liệu hiện tại của tôi, các lỗi sau xảy ra:

> 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

phương pháp Streamlined thêm 2014/10/09:

raster::crop() thể được sử dụng để cắt Spatial* (cũng như Raster*) các đối tượng.

Ví dụ, dưới đây là cách bạn có thể sử dụng nó để cắt một đối tượng 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") 

---

gốc (và vẫn hoạt động) trả lời:

Các rgeos chức năng gIntersection() làm cho điều này khá đơn giản.

Sử dụng ví dụ tiện lợi mnel như một nhảy tắt điểm:

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

Answer 2

thấy cây trồng

corp (x, y, filename = "", chụp = 'gần? ', datatype = NULL, ...)

x Raster * đối tượng

y Mức đối tượng, hoặc bất kỳ đối tượng mà từ đó một đối tượng Mức độ có thể chiết xuất (xem chi tiết

Bạn cần phải rasterize SpatialPolygon đầu tiên, sử dụng rasterize chức năng từ gói raster

tôi tạo ra một số dữ liệu để hiển thị cách sử dụng 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 

nếu tôi cố gắng:

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

Bây giờ sử dụng 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) 

Answer 3

Bạn không thể sử dụng cây trồng trên các đối tượng sp đa giác. Bạn sẽ cần phải tạo một đa giác đại diện cho các tọa độ bbox của dat2 và sau đó có thể sử dụng gIntersects trong thư viện rgeos.

Chỉnh sửa: Nhận xét này liên quan đến phiên bản có sẵn trong năm 2012 và điều này không còn là trường hợp.

Answer 4

Dưới đây là một ví dụ về cách để làm điều này với rgeos sử dụng bản đồ thế giới là một ví dụ

này xuất phát từ Roger Bivand trên R-sig-Geo mailing list. Roger là một trong những tác giả của gói sp.

Sử dụng bản đồ thế giới là một ví dụ

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