Cómo calcular los centros de los polígonos en R (para formas no contiguas)

Question

He pasado un poco de tiempo averiguando la respuesta a esta pregunta. No es inmediatamente obvia de un Búsqueda en Google así que pensé que podría ser útil poner la respuesta aquí. También hay una pregunta adicional sobre polígonos no contiguos .

Respuesta fácil e instantánea: usa el comando:

centroids <- getSpPPolygonsLabptSlots(polys)

(Esto se encontró en el descripción de la clase de la clase de datos SpatialPolygonsDataFrame R para el paquete espacial global en R, sp )

Esto parece hacer exactamente lo mismo que

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

en el siguiente código, que debería ser replicable en cualquier instalación R (¡pruébalo!)

#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")

Donde los centavos (azul) y los centroides (rojo) son centroides idénticos (esto debería aparecer después de ejecutar el código):

Hasta ahora todo bien. Pero cuando se calculan los centros de los polígonos en QGIS (menú: Vector | Geometría | Centros de los polígonos), hay resultados ligeramente diferentes para los polígonos no contiguos:

Así que esta pregunta es tres cosas:

1. Una respuesta rápida y fácil
2. Una advertencia para las personas que usan R para calcular los centros de los polígonos no contiguos
3. Una pregunta sobre cómo debe hacerse en R para contabilizar adecuadamente los polígonos de varias partes (no contiguos)

Answer 1

En primer lugar, no puedo encontrar ninguna documentación que diga que coordinates o getSpPPolygonsLabptSlots devuelve el centro de la masa. De hecho, esta última función aparece ahora como 'Deprecated' y debería emitir una advertencia.

Lo que se quiere para calcular el centroide como centro de masa de una característica es el gCentroid de la función rgeos paquete. Haciendo help.search("centroid") habrá encontrado esto.

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

debería mostrar la diferencia, y ser la misma que los centroides de Qgis.

Answer 2

aquí hay un enfoque que utiliza sf. Como demuestro, los resultados de sf::st_centroid y rgeos::gCentroid son los mismos.

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

Answer 3

Lo que hice para superar este problema es generar una función que amortigüe negativamente el polígono hasta que sea lo suficientemente pequeño como para esperar un polígono convexo. La función a utilizar es centroid(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
}