r - Calculate Centroid WITHIN / INSIDE a SpatialPolygon

Question

In Software like ArcMap one can create centroids for polygons within a polygon. In cases like the one shown below this is necessary.

In R it is possible to calculate centroids of spatial polygons with rgeos::gCentroid(). However there is no way to force the calculation of centroids within the polygon.

library(rgdal)
library(rgeos)

x <- readWKT("POLYGON ((1441727.5096940901130438 6550163.0046194596216083, 
             1150685.2609429201111197 6669225.7427449300885201, 
             975398.4520359700545669 6603079.7771196700632572, 
             866257.6087542800232768 6401334.5819626096636057, 
             836491.9242229099618271 6106985.0349301798269153, 
             972091.1537546999752522 5835786.5758665995672345, 
             1547561.0546945100650191 5782869.8033663900569081, 
             1408654.5268814601004124 5600968.3978968998417258, 
             720736.4843787000281736 5663807.0652409195899963, 
             598366.4479719599476084 6001151.4899297598749399, 
             654590.5187534400029108 6341803.2128998702391982, 
             869564.9070355399744585 6784981.1825891500338912, 
             1451649.4045378800947219 6788288.4808704098686576, 
             1441727.5096940901130438 6550163.0046194596216083))")
plot(x)

This is the polygon x

gCentroid() creates a centroid which in this specific case is located outside of the polygon. Despite being geometrically correct, some applications require centroids within the polygon, as they can be calculated by ArcMap.

xCent <- gCentroid(x, byid = TRUE)
points(xCent, col = "red", pch = 16)

A desired output (from ArcMap) looks like this:

Is there any possibility to generate centroids like this in R?

EDIT:

After some digging, it turns out that ArcMap picks a random point within the Polygon:

"For an input polygon: the output point will be inside the polygon."

Thus the question has to be: is there a function that creates a point at any random position WITHIN the polygons?

Answer 1

sf solution

With the advent of the sf package, things got a bit easier. Just use:

library(sf)
y <- st_as_sf(x) # only necessary when you don't already have an sf object
st_point_on_surface(y)

It "returns a point guaranteed to be on the (multi)surface."

sp solution

As pointed out in the updates of the Question, it seems that ArcMap is just putting a point at a random location within the polygon. This can be achieved by gPointsOnSurface(..., n = 1, type = 'random') as well.

xCent2 <- gPointOnSurface(x, byid = T)
points(xCent2, col = "blue", pch = 16)

I wrote this function which first finds the centroid and, if it is not on within (i.e. it does not overlap / intersect the polygon), it is substituted by a point on the surface. Furhtermore, it returns a new column which indicates if a point is the real centroid or not.

gCentroidWithin <- function(pol) {
  require(rgeos)

  pol$.tmpID <- 1:length(pol)
  # initially create centroid points with gCentroid
  initialCents <- gCentroid(pol, byid = T)

  # add data of the polygons to the centroids
  centsDF <- SpatialPointsDataFrame(initialCents, pol@data)
  centsDF$isCentroid <- TRUE

  # check whether the centroids are actually INSIDE their polygon
  centsInOwnPoly <- sapply(1:length(pol), function(x) {
    gIntersects(pol[x,], centsDF[x, ])
  })

  if(all(centsInOwnPoly) == TRUE){
        return(centsDF)
    }

  else {
    # substitue outside centroids with points INSIDE the polygon
    newPoints <- SpatialPointsDataFrame(gPointOnSurface(pol[!centsInOwnPoly, ], 
                                                        byid = T), 
                                        pol@data[!centsInOwnPoly,])
    newPoints$isCentroid <- FALSE
    centsDF <- rbind(centsDF[centsInOwnPoly,], newPoints)

    # order the points like their polygon counterpart based on `.tmpID`
    centsDF <- centsDF[order(centsDF$.tmpID),]

    # remove `.tmpID` column
    centsDF@data <- centsDF@data[, - which(names(centsDF@data) == ".tmpID")]

    cat(paste(length(pol), "polygons;", sum(centsInOwnPoly), "actual centroids;", 
              sum(!centsInOwnPoly), "Points corrected 
"))

    return(centsDF)
  }