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

41

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

Zentroide berechnet nach R

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:

QGIS generierte 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
RobinLovelace
quelle
Ich muss wissen, wie ich den oben erläuterten Funktionsschwerpunkt zitieren kann. Vielen Dank
Santiago Fernandez
Willkommen bei GIS StackExchange! Als neuer Benutzer nehmen Sie bitte an der Tour teil . Dies scheint eher eine neue Frage als eine Antwort auf diese Frage zu sein. Bitte poste als neue Frage.
Smiller

Antworten:

56

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.

Raumfahrer
quelle
3
Laut Roger Bivand, Entwickler einer Reihe von räumlichen Paketen von R, gilt Folgendes: "Ja. Die Klassendokumentation unter?" Polygons-class "besagt nicht, dass dies der Fall ist, da andere Punkte möglicherweise als Bezeichnungspunkte gültig eingefügt werden. Der Standardkonstruktor verwendet den Schwerpunkt des größten Rings ohne Loch im Polygons-Objekt. " - Erklärt die Nichtzusammenhängigkeit. stat.ethz.ch/pipermail/r-help/2009-February/187436.html . Bestätigt: gCentroid (sids, byid = TRUE) löst das Problem tatsächlich.
Robin Lovelace
funktioniert bei mir nicht ... selbst wenn gCentroid (polygon, byid = TRUE) angewendet wird, befindet sich mein Schwerpunkt zwischen zwei Polygonen. Ich gehe also davon aus, dass diese als mehrteilige Polygone betrachtet werden. Wie kann ich sie aufteilen? Die Punkte (Koordinaten (SC.tracks), pch = 16, col = "blue", cex = 0.4) erzeugen jedoch keinen Schwerpunkt aus Polygon ... Danke!
Maycca
Der Link zu stat.ethz.ch funktioniert nicht mehr. Der Vollständigkeit halber bin ich mir fast sicher, dass die Antwort jetzt hier zu finden ist: r.789695.n4.nabble.com/…
Exocom
8

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

Bildbeschreibung hier eingeben

sebdalgarno
quelle
3

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
}
Gert
quelle
Langsam ergibt sich aber ein sehr schönes Ergebnis. Es ist gut zentriert und liefert ein gutes Ergebnis für die Platzierung der Etiketten
Bastien