Làm thế nào để clip WorldMap với đa giác trong R?

Question

Tôi đã nhập một tập dữ liệu bản đồ thế giới từ www.GADM.org bằng cách sử dụng trình bao gói R. Tôi muốn cắt nó thành đa giác mà tôi tạo để giảm kích thước của bản đồ. Tôi có thể lấy dữ liệu và tôi có thể tạo ra đa giác không có vấn đề, nhưng khi tôi sử dụng lệnh 'gIntersection' tôi nhận được một thông báo lỗi tối nghĩa.

Bất kỳ đề xuất nào về cách nén tập dữ liệu Bản đồ thế giới của tôi?

library(raster) 
library(rgeos) 

## Download Map of the World ## 
WorldMap <- getData('countries') 

## Create the clipping polygon 
clip.extent <- as(extent(-20, 40, 30, 72), "SpatialPolygons") 
proj4string(clip.extent) <- CRS(proj4string(WorldMap)) 

## Clip the map 
EuropeMap <- gIntersection(WorldMap, clip.extent, byid = TRUE) 

Thông báo lỗi:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : 
    Geometry collections may not contain other geometry collections 
In addition: Warning message: 
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : 
    spgeom1 and spgeom2 have different proj4 strings 

Answer 1

Bạn không cần phải sử dụng PBS (Tôi đã học được điều này một cách khó khăn, như liên kết r-sig-geo đăng bởi @flowla là một câu hỏi ban đầu được đăng bởi tôi!). Mã này cho thấy cách thực hiện tất cả trong rgeos, nhờ vào số various khác nhau postings từ Roger Bivand. Đây sẽ là cách rút kinh điển hơn, mà không cần phải ép buộc đối tượng PolySet.

Lý do cho thông báo lỗi là bạn không thể gIntersection trên một bộ sưu tập của SpatialPolygons, bạn cần phải làm chúng riêng lẻ. Tìm hiểu xem bạn muốn sử dụng cái nào bằng cách sử dụng gIntersects. Sau đó, tôi đặt từng đa giác theo quốc gia bằng cách sử dụng gIntersection. Tôi gặp vấn đề khi chuyển một danh sách các đối tượng SpatialPolygons quay lại SpatialPolygons, biến các shapefile bị cắt xén thành SpatialPolygons, đó là vì không phải tất cả các đối tượng đã cắt là classSpatialPolygons. Khi chúng tôi loại trừ tất cả mọi thứ hoạt động tốt.

# This very quick method keeps whole countries 
gI <- gIntersects(WorldMap , clip.extent , byid = TRUE) 
Europe <- WorldMap[ which(gI) , ] 
plot(Europe) 


#If you want to crop the country boundaries, it's slightly more involved: 
# This crops countries to your bounding box 
gI <- gIntersects(WorldMap , clip.extent , byid = TRUE) 
out <- lapply(which(gI) , function(x){ gIntersection(WorldMap[x,] , clip.extent) }) 

# But let's look at what is returned 
table(sapply(out , class)) 
SpatialCollections SpatialPolygons 
       2     63 

# We want to keep only objects of class SpatialPolygons     
keep <- sapply(out, class) 
out <- out[keep == "SpatialPolygons"] 


# Coerce list back to SpatialPolygons object 
Europe <- SpatialPolygons(lapply(1:length(out) , function(i) { Pol <- slot(out[[i]], "polygons")[[1]]; slot(Pol, "ID") <- as.character(i) 
    Pol 
})) 

plot(Europe) 

Nếu có thể, tôi khuyên bạn nên nhìn vào naturalearthdata. Họ có shapefiles chất lượng cao được giữ up-to-date và liên tục kiểm tra lỗi (vì chúng là mã nguồn mở nếu bạn tìm thấy một lỗi làm email cho họ). Ranh giới quốc gia nằm dưới các nút Văn hóa. Bạn sẽ thấy chúng cũng nhẹ hơn một chút và bạn có thể chọn độ phân giải phù hợp với nhu cầu của mình.

Answer 2

Làm thế nào về một bước trung gian ít? Tôi đã sử dụng mã sau chủ yếu từ R-sig-Geo và tôi nghĩ rằng nó nên làm các trick. Bạn sẽ cần cả hai gói 'maptools' và 'PBSmapping', vì vậy hãy đảm bảo bạn đã cài đặt chúng. Đây là mã của tôi:

# Required packages 
library(raster) 
library(maptools) 
library(PBSmapping) 

# Download world map 
WorldMap <- getData('countries') 
# Convert SpatialPolygons to class 'PolySet' 
WorldMap.ps <- SpatialPolygons2PolySet(WorldMap) 
# Clip 'PolySet' by given extent 
WorldMap.ps.clipped <- clipPolys(WorldMap.ps, xlim = c(-20, 40), ylim = c(30, 72)) 
# Convert clipped 'PolySet' back to SpatialPolygons 
EuropeMap <- PolySet2SpatialPolygons(WorldMap.ps.clipped, close_polys=TRUE) 

Tôi vừa thử nghiệm và nó hoạt động mà không có bất kỳ vấn đề gì. Phải mất một số thời gian tính toán để biến SpatialPolygons thành PolySet.

Chúc mừng, Florian