合并SpatialPolygonDataframe的多边形

Question

我已经创建了一个SpatialPolygonsDataFrame，其中的国家与区域相关联。这里是什么样子https://cloudstor.aarnet.edu.au/plus/index.php/s/RpYr3xyMrmhaGKA（有数据EDU链接）：

由于对象太大用，并作为我只是在信息的区域，而不是国家的相关工作，我要合并各地区内的国家。换句话说，我希望为每个区域创建一个多边形。 我试图运行下面的命令，但运行需要几个小时，我从来没有能够看到他们是否正常工作。

library(rgdal) 
regions = readOGR("./regionscountriesSTACK.shp") 
library(maptools) 
regions <- unionSpatialPolygons(regions, IDs=regions@data$REGION) 
library(rgeos) 
gUnionCascaded(regions, id = regions@data$REGION) 
gUnaryUnion(regions, id = regions@data$REGION) 

任何建议有效的方式来处理？非常感谢！

Answer 1

你需要： 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(subPol@hole && subPol@area < limitHoles) 
      keep <- FALSE 
     else 
      keep <- TRUE 
     return(keep) 
    } 
    nPol <- length(SPol) 
    newPols <- list() 
    for(iPol in 1:nPol){ 
     subPolygons <- list() 
     pol <- SPol@polygons[[iPol]] 
     for(iSubPol in 1:length(pol@Polygons)){ 
      subPol <- pol@Polygons[[iSubPol]] 

      if(fn(subPol)) 
       subPolygons[[length(subPolygons)+1]] <- subPol 
     } 
     newPols[[length(newPols)+1]] <- Polygons(subPolygons,pol@ID) 
    } 
    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:",subReg@data[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 
countries@data$COUNTRY <- as.character(countries@data$COUNTRY) 
countries@data$REGION <- as.character(countries@data$REGION) 

# aggregate region by region 
vectRegions <- unique(countries@data$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(countriesLow@data$REGION)<-c(levels(countriesLow@data$REGION),"Other") 
countriesLow@data$REGION[which(is.na(countriesLow@data$REGION))] <- "Other" 

vectRegions <- unique(countriesLow@data$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兆字节）。