Comment calculer les centroïdes de polygones en R (pour les formes non contiguës)

41

J'ai passé un peu de temps à trouver la réponse à cette question. Ce n'est pas immédiatement évident à partir d'une recherche Google. Nous avons donc pensé qu'il serait utile de poster la réponse ici. Il y a aussi une question supplémentaire sur les polygones non contigus .

Réponse facile instantanée: utilisez la commande:

centroids <- getSpPPolygonsLabptSlots(polys)

(Cela a été trouvé dans la description de classe de la classe de données SpatialPolygonsDataFrame R pour le package spatial global dans R, sp )

Cela semble faire exactement la même chose que

cents <- SpatialPointsDataFrame(coords=cents, data=sids@data, proj4string=CRS("+proj=longlat +ellps=clrk66"))

dans le code suivant, qui devrait être réplicable sur toute installation de R (essayez-le!)

#Rcentroids
install.packages("GISTools")
library(GISTools)
sids <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], 
                      proj4string=CRS("+proj=longlat +ellps=clrk66"))
class(sids)
plot(sids)
writeSpatialShape(sids, "sids")
cents <- coordinates(sids)
cents <- SpatialPointsDataFrame(coords=cents, data=sids@data, 
                  proj4string=CRS("+proj=longlat +ellps=clrk66"))
points(cents, col = "Blue")
writeSpatialShape(cents, "cents")

centroids <- getSpPPolygonsLabptSlots(sids)
points(centroids, pch = 3, col = "Red")

Où les cent (bleu) et les centroïdes (rouge) sont des centroïdes identiques (le tracé devrait apparaître après l'exécution du code):

centroïdes calculés par R

Jusqu'ici tout va bien. Mais lorsque vous calculez des centroïdes de polygones dans QGIS (menu: Vecteur | Géométrie | Centroïdes polygonaux), les résultats obtenus pour les polygones non contigus sont légèrement différents:

Polygones générés par QGIS

Donc, cette question est 3-choses:

  1. Une réponse rapide et facile
  2. Avertissement pour les personnes utilisant R pour calculer les centroïdes de polygones non contigus
  3. Une question sur la manière de procéder en R pour prendre en compte correctement les polygones à plusieurs parties (non contigus)
RobinLovelace
la source
J'ai besoin de savoir comment citer la fonction centroïde expliquée ci-dessus. Merci
Santiago Fernandez le
Bienvenue sur GIS StackExchange! En tant que nouvel utilisateur, veuillez faire le tour . Cela semble être une nouvelle question plutôt qu'une réponse à cette question. S'il vous plaît poster comme une nouvelle question.
smiller

Réponses:

56

Premièrement, je ne trouve pas de documentation indiquant cela coordinatesou getSpPPolygonsLabptSlotsrenvoyant le centre de masse du centre de masse. En fait, cette dernière fonction apparaît désormais comme «Obsolète» et devrait émettre un avertissement.

Ce que vous voulez pour calculer le centroïde en tant que centre de gravité d'une fonction est la gCentroidfonction du rgeospackage. Faire help.search("centroid")aura trouvé cela.

trueCentroids = gCentroid(sids,byid=TRUE)
plot(sids)
points(coordinates(sids),pch=1)
points(trueCentroids,pch=2)

devrait montrer la différence, et être le même que les centroïdes Qgis.

Spacedman
la source
3
Selon Roger Bivand, développeur de plusieurs paquets spatiaux de R, il déclare: "Oui. La documentation de la classe sur?" Polygons-class "n'indique pas que c'est le cas, car d'autres points pourraient être valablement insérés en tant que points d'étiquette. Le constructeur par défaut utilise le centroïde du plus grand anneau non-trou de l’objet Polygons. " - Explique la non-contiguïté. stat.ethz.ch/pipermail/r-help/2009-Février/187436.html . Confirmé: gCentroid (Sid, byid = TRUE) résout effectivement le problème.
RobinLovelace
cela ne fonctionne pas pour moi ... même si j'applique gCentroid (polygone, byid = TRUE), mon centroïde est un espace situé entre deux polygones. Je suppose donc que ceux-ci sont considérés comme des polygones en plusieurs parties? Comment puis-je les séparer? les points (coordonnées (SC.tracks), pch = 16, col = "blue", cex = 0.4), cependant, produire ne produit pas de centre de gravité en dehors du polygone ... merci!
Maycca
Le lien vers stat.ethz.ch ne fonctionne plus. Juste par souci d'exhaustivité, je suis presque sûr que la réponse peut maintenant être trouvée ici: r.789695.n4.nabble.com/…
Exocom
8

Voici une approche utilisant SF. Comme je le démontre, les résultats de sf :: st_centroid et de rgeos :: gCentroid sont identiques.

library(sf)
library(ggplot2)

# I transform to utm because st_centroid is not recommended for use on long/lat 
nc <- st_read(system.file('shape/nc.shp', package = "sf")) %>% 
  st_transform(32617)

# using rgeos
sp_cent <- gCentroid(as(nc, "Spatial"), byid = TRUE)

# using sf
sf_cent <- st_centroid(nc)

# plot both together to confirm that they are equivalent
ggplot() + 
  geom_sf(data = nc, fill = 'white') +
  geom_sf(data = sp_cent %>% st_as_sf, color = 'blue') + 
  geom_sf(data = sf_cent, color = 'red') 

entrez la description de l'image ici

sebdalgarno
la source
3

Ce que j’ai fait pour résoudre ce problème est de générer une fonction qui tamponne négativement le polygone jusqu’à ce qu’il soit suffisamment petit pour s’attendre à un polygone convexe. La fonction à utiliser estcentroid(polygon)

#' find the center of mass / furthest away from any boundary
#' 
#' Takes as input a spatial polygon
#' @param pol One or more polygons as input
#' @param ultimate optional Boolean, TRUE = find polygon furthest away from centroid. False = ordinary centroid

require(rgeos)
require(sp)

centroid <- function(pol,ultimate=TRUE,iterations=5,initial_width_step=10){
  if (ultimate){
    new_pol <- pol
    # For every polygon do this:
    for (i in 1:length(pol)){
      width <- -initial_width_step
      area <- gArea(pol[i,])
      centr <- pol[i,]
      wasNull <- FALSE
      for (j in 1:iterations){
        if (!wasNull){ # stop when buffer polygon was alread too small
          centr_new <- gBuffer(centr,width=width)
          # if the buffer has a negative size:
          substract_width <- width/20
          while (is.null(centr_new)){ #gradually decrease the buffer size until it has positive area
            width <- width-substract_width
            centr_new <- gBuffer(centr,width=width)
            wasNull <- TRUE
          }
          # if (!(is.null(centr_new))){
          #   plot(centr_new,add=T)
          # }
          new_area <- gArea(centr_new)
          #linear regression:
          slope <- (new_area-area)/width
          #aiming at quarter of the area for the new polygon
          width <- (area/4-area)/slope
          #preparing for next step:
          area <- new_area
          centr<- centr_new
        }
      }
      #take the biggest polygon in case of multiple polygons:
      d <- disaggregate(centr)
      if (length(d)>1){
        biggest_area <- gArea(d[1,])
        which_pol <- 1                             
        for (k in 2:length(d)){
          if (gArea(d[k,]) > biggest_area){
            biggest_area <- gArea(d[k,])
            which_pol <- k
          }
        }
        centr <- d[which_pol,]
      }
      #add to class polygons:
      new_pol@polygons[[i]] <- remove.holes(new_pol@polygons[[i]])
      new_pol@polygons[[i]]@Polygons[[1]]@coords <- centr@polygons[[1]]@Polygons[[1]]@coords
    }
    centroids <- gCentroid(new_pol,byid=TRUE)
  }else{
    centroids <- gCentroid(pol,byid=TRUE)  
  }  
  return(centroids)
}

#Given an object of class Polygons, returns
#a similar object with no holes


remove.holes <- function(Poly){
  # remove holes
  is.hole <- lapply(Poly@Polygons,function(P)P@hole)
  is.hole <- unlist(is.hole)
  polys <- Poly@Polygons[!is.hole]
  Poly <- Polygons(polys,ID=Poly@ID)
  # remove 'islands'
  max_area <- largest_area(Poly)
  is.sub <- lapply(Poly@Polygons,function(P)P@area<max_area)  
  is.sub <- unlist(is.sub)
  polys <- Poly@Polygons[!is.sub]
  Poly <- Polygons(polys,ID=Poly@ID)
  Poly
}
largest_area <- function(Poly){
  total_polygons <- length(Poly@Polygons)
  max_area <- 0
  for (i in 1:total_polygons){
    max_area <- max(max_area,Poly@Polygons[[i]]@area)
  }
  max_area
}
Gert
la source
Lent mais donne de très bons résultats. C'est bien centré et donne un bon résultat pour le placement de l'étiquette
Bastien