Fusion de polygones d'un objet SpatialPolygonDataframe

Question

J'ai créé un SpatialPolygonsDataFrame dans lequel des pays sont associés à des régions. Voici à quoi il ressemble https://cloudstor.aarnet.edu.au/plus/index.php/s/RpYr3xyMrmhaGKA (lien EDU avec des données):

Comme l'objet est trop grand pour travailler avec et comme je ne suis intéressé par les régions d'information et non les pays, je tiens à fusionner les pays de chaque région. En d'autres termes, je souhaite créer un seul polygone pour chaque région. J'ai essayé d'exécuter les commandes suivantes, mais les exécutions prennent des heures et je n'ai jamais pu voir si elles fonctionnaient correctement.

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) 

Un conseil sur un moyen efficace de traiter? Merci beaucoup!

Answer 1

Vous devez: 1) simplifier les polygones, ils sont probablement trop complexes pour commencer, surtout si vous voulez les agréger par région. Utilisez gSimplify à partir du paquetage rgeos. Difficile de vous aider plus loin sans les données. 2) enlever les trous, ils prennent beaucoup d'espace et conduisent à des problèmes lorsque vous simplifiez 3) font l'union et la simplification avec permettant une simplification sérieuse des données

Les codes suivants combine tout cela, faire pays par pays permet également de voir les progrès de la chose:

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

Cette solution est implémentée dans le package spatDataManagement. Aussi, vous pouvez utiliser rworldmap pour une carte du monde plus légère si vous êtes seulement intéressé par les régions. Il devient alors:

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) 

Et cette nouvelle carte est beaucoup, beaucoup plus petites (2,8 Mégaoctets).