Search code examples
rr-sfcentroid

r sf package centroid within polygon


I need to add labels to polygons and I normally use the centroid, however the centroid does not fall inside the polygon. I found this question Calculate Centroid WITHIN / INSIDE a SpatialPolygon but I'm using the sf package.

Below is a toy data

rm(list = ls(all = TRUE)) #start with empty workspace

library(sf)
library(tidyverse)
library(ggrepel)

pol <- st_polygon(list(rbind(c(144, 655),c(115, 666)
                         ,c(97, 660),c(86, 640)
                         ,c(83, 610),c(97, 583)
                         ,c(154, 578),c(140, 560)
                         ,c(72, 566),c(59, 600)
                         ,c(65, 634),c(86, 678)
                         ,c(145, 678),c(144, 655)))) %>%
  st_sfc()

a = data.frame(NAME = "A")
st_geometry(a) = pol

a <- a  %>% 
  mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
     lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

ggplot() +
  geom_sf(data = a, fill = "orange") +
  geom_label_repel(data = a, aes(x = lon, y = lat, label = NAME)) 

which results in the following

enter image description here


Solution

  • The simple answer is to replace st_centroid with st_point_on_surface. This won't return the true centroid in cases where the centroid is inside the polygon.

    a2 <- a  %>% 
      mutate(lon = map_dbl(geometry, ~st_point_on_surface(.x)[[1]]),
             lat = map_dbl(geometry, ~st_point_on_surface(.x)[[2]]))
    
    ggplot() +
      ggplot2::geom_sf(data = a2, fill = "orange") +
      geom_label_repel(data = a2, aes(x = lon, y = lat, label = NAME))
    

    Alternatively

    If the polygon has a centroid that is inside the polygon, use that, otherwise, find a point within the polygon.

    st_centroid_within_poly <- function (poly) {
    
      # check if centroid is in polygon
      centroid <- poly %>% st_centroid() 
      in_poly <- st_within(centroid, poly, sparse = F)[[1]] 
    
      # if it is, return that centroid
      if (in_poly) return(centroid) 
    
      # if not, calculate a point on the surface and return that
      centroid_in_poly <- st_point_on_surface(poly) 
      return(centroid_in_poly)
    }
    
    a3 <- a  %>% 
      mutate(lon = map_dbl(geometry, ~st_centroid_within_poly(.x)[[1]]),
             lat = map_dbl(geometry, ~st_centroid_within_poly(.x)[[2]]))
    
    ggplot() +
      ggplot2::geom_sf(data = a3, fill = "orange") +
      geom_label_repel(data = a3, aes(x = lon, y = lat, label = NAME)) 
    

    The function above st_centroid_within_polygon is adapted from the question you reference for the sf package. A more thorough review of how st_point_on_surface works can be found here.