Cómo calcular los centroides poligonales en R (para formas no contiguas)

37

He pasado un poco de tiempo resolviendo la respuesta a esta pregunta. No es inmediatamente obvio de un Búsqueda en Google , por lo que pensé que puede ser útil publicar la respuesta aquí. También hay una pregunta adicional sobre polígonos no contiguos .

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

centroids <- getSpPPolygonsLabptSlots(polys)

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

Esto parece hacer exactamente lo mismo que

cents <- SpatialPointsDataFrame(coords=cents, [email protected], proj4string=CRS("+proj=longlat +ellps=clrk66"))

en el siguiente código, que debe ser replicable en cualquier instalación R (¡inténtalo!)

#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, [email protected], 
                  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 haber ejecutado el código):

Hastaahoratodobien.PerocuandosecalculanloscentroidespoligonalesenQGIS(menú:Vector|Geometry|PolygonCentroids),hayresultadosligeramentediferentesparapolígonosnocontiguos:

Entonces esta pregunta es 3 cosas:

  1. Una respuesta rápida y fácil
  2. Una advertencia para las personas que usan R para calcular los centroides para polígonos no contiguos
  3. Una pregunta sobre cómo se debe hacer en R para tener en cuenta adecuadamente los polígonos de varias partes (no contiguos)
pregunta RobinLovelace 08.12.2012 - 20:59

3 respuestas

50

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

Lo que desea para calcular el centroide como el centro de masa de una entidad es la función gCentroid del paquete rgeos . Al hacer help.search("centroid") habremos 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 el mismo que los centroides Qgis.

    
respondido por el Spacedman 09.12.2012 - 12:03
3

Lo que hice para superar este problema es generar una función que amortigua 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:
      [email protected][[i]] <- remove.holes([email protected][[i]])
      [email protected][[i]]@Polygons[[1]]@coords <- [email protected][[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([email protected],function(P)[email protected])
  is.hole <- unlist(is.hole)
  polys <- [email protected][!is.hole]
  Poly <- Polygons(polys,[email protected])
  # remove 'islands'
  max_area <- largest_area(Poly)
  is.sub <- lapply([email protected],function(P)[email protected]<max_area)  
  is.sub <- unlist(is.sub)
  polys <- [email protected][!is.sub]
  Poly <- Polygons(polys,[email protected])
  Poly
}
largest_area <- function(Poly){
  total_polygons <- length([email protected])
  max_area <- 0
  for (i in 1:total_polygons){
    max_area <- max(max_area,[email protected][[i]]@area)
  }
  max_area
}
    
respondido por el Gert 15.12.2017 - 10:22
3

aquí hay un enfoque usando 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') 

    
respondido por el sebdalgarno 02.03.2018 - 21:07

Lea otras preguntas en las etiquetas