你需要: 1)简化多边形,他们可能太复杂,特别是如果你想按地区聚合他们。使用rgeos软件包中的gSimplify
。没有数据,很难帮助你。 2)去除孔,而且占用了大量的空间,并导致烦恼,当你简化 3)做好工会和简化与允许数据
以下代码的严重简化合并了这一切,按国家这样做的国家也允许看到的事情的进展:
library(rgdal)
library(maptools)
library(rgeos)
# Remove all holes that are bigger than "limitholes", by default all of them
RemoveHoles <- function(SPol,limitHoles=+Inf){
fn <- function(subPol){
if([email protected] && [email protected] < limitHoles)
keep <- FALSE
else
keep <- TRUE
return(keep)
}
nPol <- length(SPol)
newPols <- list()
for(iPol in 1:nPol){
subPolygons <- list()
pol <- [email protected][[iPol]]
for(iSubPol in 1:length([email protected])){
subPol <- [email protected][[iSubPol]]
if(fn(subPol))
subPolygons[[length(subPolygons)+1]] <- subPol
}
newPols[[length(newPols)+1]] <- Polygons(subPolygons,[email protected])
}
newSPol <- SpatialPolygons(newPols,proj4string=CRS(proj4string(SPol)))
# SPolSimple <- gSimplify(newSPol,tol=0.01)
newSPol <- createSPComment(newSPol)
return(newSPol)
}
# Union Polygon (country) by polygon for a given region
UnionSimplifyPolByPol <- function(subReg,precision=0.2){
# if(length(subReg)>1){
# out <- unionSpatialPolygons(subReg[1:2,],rep(1,2),threshold=0.1)
# out <- RemoveHoles(out)
# }
out <- c()
for(iCountry in 1:length(subReg)){
cat("Adding:",[email protected][iCountry,"COUNTRY"],"\n")
toAdd <- gSimplify(as(subReg[iCountry,],"SpatialPolygons"),
tol=precision,topologyPreserve=TRUE)
toAdd <- RemoveHoles(toAdd)
if(is.null(out)){
out <- toAdd
}else{
toUnite <- rbind(out,toAdd)
out <- unionSpatialPolygons(toUnite,
IDs=rep(1,2),
threshold=precision)
}
out <- RemoveHoles(out)
# plot(out)
}
return(out)
}
# import the data
countries <- readOGR("regionscountriesSTACK.shp")
# you don't want to be bothered by factors
[email protected]$COUNTRY <- as.character([email protected]$COUNTRY)
[email protected]$REGION <- as.character([email protected]$REGION)
# aggregate region by region
vectRegions <- unique([email protected]$REGION)
world <- c()
for(iRegion in 1:length(vectRegions)){
regionName <- vectRegions[iRegion]
cat("Region:",regionName)
region <- UnionSimplifyPolByPol(countries[which(countries$REGION==regionName),])
region <- spChFIDs(region,regionName)
if(is.null(world)){
world <- region
}else{
world <- rbind(world,region)
}
plot(world)
}
该解决方案在包spatDataManagement实现。如果您只对区域感兴趣,也可以使用rworldmap
获取较轻的世界地图。然后,它变成了:
library("spatDataManagement")
library("rworldmap")
levels([email protected]$REGION)<-c(levels([email protected]$REGION),"Other")
[email protected]$REGION[which(is.na([email protected]$REGION))] <- "Other"
vectRegions <- unique([email protected]$REGION)
world <- c()
for(iRegion in 1:length(vectRegions)){
regionName <- vectRegions[iRegion]
cat("Region:",regionName)
region <- UnionSimplifyPolByPol(countriesLow[which(countriesLow$REGION==regionName),])
region <- spChFIDs(region,gsub(" ","",regionName)) # IDs can't have spaces
if(is.null(world)){
world <- region
}else{
world <- rbind(world,region)
}
}
world <- SpatialPolygonsDataFrame(world,data.frame(id=1:length(world),name=vectRegions),match.ID=FALSE)
plot(world,col=world$id)
而这个新地图是非常非常小(2.8兆字节)。
感谢您的建议。我刚刚添加链接到shapefile。 – Cecile