RのポリゴンでWorldMapをクリップする方法は？

Question

私は、Rパッケージラスタを使用してwww.GADM.orgから世界地図データセットをインポートしました。マップのサイズを小さくするために作成したポリゴンにクリップしたいと思います。私はデータを取得することができ、問題なくポリゴンを作成できますが、 'gIntersection'コマンドを使用するとわかりにくいエラーメッセージが表示されます。

私のワールドマップデータセットをクリップする方法に関する提案はありますか？

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) 

エラーメッセージ：少しの中間ステップに関する

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

あなたはPBSを使用する必要はありませんが（@flowlaによって投稿r-sig-geoリンクは、もともと私が掲示する質問だったので、私は、このハードな方法を学びました！）。このコードは、Roger Bivandのvariousと異なるpostingsのおかげで、rgeosですべてを行う方法を示しています。これはPolySetオブジェクトへの強制に頼ることなく、より標準的なサブセット化の方法になります。

エラーメッセージが表示されるのは、SpatialPolygonsのコレクション上でgIntersectionを実行できないため、個別に行う必要があるということです。 gIntersectsを使用して、必要なものを見つけます。次に、gIntersectionを使用して各国のポリゴンをサブセット化します。 SpatialPolygonsオブジェクトのリストをSpatialPolygonsに戻して問題を起こしました。これは、クロップされたシェープファイルをSpatialPolygonsに変換します。これは、切り取られたオブジェクトのすべてがclassSpatialPolygonsではないためです。これらを除外すると、すべて正常に動作します。

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

あなたができるなら、私はあなたがnaturalearthdataを見てお勧めします。彼らは最新の状態に保たれた高品質のシェイプファイルを持っています。エラーが見つかった場合はオープンソースなので、エラーがないか常にチェックされます。国境はCulturalボタンの下にあります。それらはさらに軽量であることがわかり、必要に応じて適切な解像度を選択することができます。

Answer 2

どのように？私はR-sig-Geoから主に以下のコードを採用しました。 'maptools'と 'PBSmapping'の両方のパッケージが必要なので、それらがインストールされていることを確認してください。ここに私のコード：

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

私はそれをテストしただけで問題なく動作しました。しかし、SpatialPolygonsをPolySetに変換するには、計算時間がかかりました。

乾杯、 フロリアン