Berechnung von Polygonschwerpunkten in R (für nicht zusammenhängende Formen)

Ich habe ein wenig Zeit damit verbracht, die Antwort auf diese Frage zu finden. Aus einer Google-Suche ist dies nicht sofort ersichtlich. Daher ist es möglicherweise hilfreich, die Antwort hier zu veröffentlichen. Es gibt auch eine zusätzliche Frage zu nicht zusammenhängenden Polygonen .

Sofortige einfache Antwort: benutze den Befehl:

centroids <- getSpPPolygonsLabptSlots(polys)

(Dies wurde in der Klassenbeschreibung der R-Datenklasse SpatialPolygonsDataFrame für das übergeordnete räumliche Paket in R, sp gefunden. )

Dies scheint genau dasselbe zu tun wie

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

im folgenden Code, der bei jeder R-Installation replizierbar sein sollte (probieren Sie es aus!)

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

Wo Cent (blau) und Centroide (rot) identische Centroide sind (dieser Plot sollte erscheinen, nachdem Sie den Code ausgeführt haben):

So weit, ist es gut. Wenn Sie jedoch in QGIS Polygonschwerpunkte berechnen (Menü: Vektor | Geometrie | Polygonschwerpunkte), ergeben sich geringfügig andere Ergebnisse für nicht zusammenhängende Polygone:

Diese Frage besteht also aus drei Dingen:

1. Eine schnelle und einfache Antwort
2. Eine Warnung für Benutzer, die R verwenden, um Schwerpunkte für nicht zusammenhängende Polygone zu berechnen
3. Eine Frage, wie in R mehrteilige (nicht zusammenhängende) Polygone richtig berücksichtigt werden sollen

Erstens kann ich keine Dokumentation finden, die dies sagt coordinatesoder getSpPPolygonsLabptSlotsden Schwerpunkt zurückgibt. Tatsächlich wird die letztere Funktion jetzt als "Veraltet" angezeigt und sollte eine Warnung ausgeben.

Was Sie für die Berechnung des Schwerpunkts als Massenmittelpunkt eines Features benötigen, ist die gCentroidFunktion aus dem rgeosPaket. Tun help.search("centroid")wird dies gefunden haben.

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

sollte den Unterschied zeigen und mit den Qgis-Zentroiden identisch sein.

Hier ist ein Ansatz mit sf. Wie ich demonstriere, sind die Ergebnisse von sf :: st_centroid und rgeos :: gCentroid identisch.

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

Um dieses Problem zu lösen, habe ich eine Funktion generiert, die das Polygon negativ puffert, bis es klein genug ist, um ein konvexes Polygon zu erwarten. Die zu verwendende Funktion istcentroid(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
}