Wie WorldMap mit Polygon in R schneiden?

Question

Ich habe einen Weltkartendatensatz von www.GADM.org mit dem R-Paket-Raster importiert. Ich möchte es auf ein Polygon abschneiden, das ich erstellt habe, um die Größe der Karte zu reduzieren. Ich kann die Daten abrufen und ich kann das Polygon kein Problem erstellen, aber wenn ich den 'gIntersection' Befehl verwende, erhalte ich eine obskure Fehlermeldung.

Haben Sie Vorschläge, wie Sie meinen World Map-Datensatz schneiden können?

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) 

Fehlermeldung:

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

Sie brauchen nicht PBS zu verwenden (Ich habe dies auf die harte Tour gelernt, wie die r-sig-geo Link von @flowla gepostet eine Frage ursprünglich von mir geschrieben war!). Dieser Code zeigt, wie man alles in RGEOS macht, dank various verschiedener postings von Roger Bivand. Dies wäre die kanonischere Art der Teilmengenbildung, ohne auf PolySet-Objekte umzuschalten.

Der Grund für die Fehlermeldung ist, dass Sie nicht auf eine Sammlung von SpatialPolygons GIntersection können, müssen Sie sie einzeln ausführen. Finden Sie heraus, welche Sie mit gIntersects verwenden möchten. Ich unterteile dann jedes Landpolygon unter Verwendung gIntersection. Ich hatte ein Problem damit, eine Liste von SpatialPolygons-Objekten zurück an SpatialPolygons zu übergeben, wodurch die gecroppten Shapefiles in SpatialPolygons umgewandelt werden, da nicht alle ausgeschnittenen Objekte von class waren. Sobald wir diese ausgeschlossen haben, funktioniert alles gut.

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

Wenn Sie können, empfehle ich Ihnen bei naturalearthdata suchen. Sie haben qualitativ hochwertige Shapefiles, die auf dem neuesten Stand gehalten werden und ständig auf Fehler überprüft werden (da sie Open Source sind, wenn Sie einen Fehler finden, mailen Sie sie). Ländergrenzen sind unter den kulturellen Schaltflächen. Sie werden sehen, dass sie auch ein wenig leichter sind und Sie können eine Auflösung wählen, die Ihren Bedürfnissen entspricht.

Answer 2

Wie wärs mit einem kleinen Zwischenschritt? Ich nahm den folgenden Code hauptsächlich von R-sig-Geo an und ich denke, dass es den Trick tun sollte. Sie benötigen sowohl 'maptools' als auch 'PBSmapping' Pakete, also stellen Sie sicher, dass Sie sie installiert haben. Hier ist mein Code:

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

Ich habe es gerade getestet und es funktionierte ohne Probleme. Es dauerte jedoch einige Zeit, um SpatialPolygons in PolySet zu konvertieren.

Cheers, Florian