Cara menghitung poligon centroid di R (untuk bentuk yang tidak bersebelahan)

Saya telah menghabiskan sedikit waktu mencari tahu jawaban untuk pertanyaan ini. Ini tidak segera terlihat dari pencarian Google , jadi mungkin bermanfaat untuk mengirim jawabannya di sini. Ada juga pertanyaan tambahan tentang poligon yang tidak berdekatan .

Jawaban mudah instan: gunakan perintah:

centroids <- getSpPPolygonsLabptSlots(polys)

(Ini ditemukan dalam deskripsi kelas dari kelas data SpatialPolygonsDataFrame R untuk paket spasial menyeluruh di R, sp )

Ini tampaknya melakukan hal yang persis sama dengan

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

dalam kode berikut, yang seharusnya dapat direplikasi pada instalasi R apa saja (coba!)

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

Di mana sen (biru) dan centroid (merah) adalah centroid identik (plot ini seharusnya muncul setelah Anda menjalankan kode):

Sejauh ini baik. Tetapi ketika Anda menghitung poligon centroid di QGIS (menu: Vektor | Geometri | Poligon Centroid), ada hasil yang sedikit berbeda untuk poligon yang tidak berdekatan:

Jadi pertanyaan ini adalah 3 hal:

1. Jawaban yang cepat dan mudah
2. Peringatan untuk orang yang menggunakan R untuk menghitung centroid untuk poligon yang tidak berdekatan
3. Sebuah pertanyaan tentang bagaimana hal itu harus dilakukan dalam R untuk memperhitungkan poligon multi-bagian (tidak berdampingan) dengan benar

Pertama, saya tidak dapat menemukan dokumentasi yang mengatakan itu coordinatesatau getSpPPolygonsLabptSlotsmengembalikan centroid pusat massa. Bahkan fungsi yang terakhir sekarang muncul sebagai 'Usang' dan harus mengeluarkan peringatan.

Apa yang Anda inginkan untuk menghitung centroid sebagai pusat fitur adalah gCentroidfungsi dari rgeospaket. Melakukan help.search("centroid")akan menemukan ini.

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

harus menunjukkan perbedaan, dan harus sama dengan centroid Qgis.

di sini adalah pendekatan menggunakan sf. Seperti yang saya tunjukkan, hasil dari sf :: st_centroid dan rgeos :: gCentroid adalah sama.

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

Apa yang saya lakukan untuk mengatasi masalah ini adalah menghasilkan fungsi yang secara negatif menyangga poligon hingga cukup kecil untuk mengharapkan poligon cembung. Fungsi yang digunakan adalahcentroid(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
}