合并 SpatialPolygonsDataframe 的多边形

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 兆字节）。