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?

Well, I wouldn't intuit the placement of centroid as ArcMap does; that said, it appears you're between here [link](https://stackoverflow.com/questions/23613655/calculating-weighted-polygon-centroids-in-r?rq=1) and here [link](https://stackoverflow.com/questions/35720614/gcentroid-rgeos-r-vs-actual-centroid-in-python), which is somewhat surprising as R placement does look like geometric mean ( as though it had been passed points), so something like mapping the center of a C or an atoll. Perhaps ArcMap subsets convex polys to subareas prior to calculating centroid-ish. — Chris, Jun 03 '17 at 14:26
Thanks for the links so far, @Chris. I think they will add up to a sufficient solution. However, I think this is a basic but important GIS operation and I have not found a solution so far. I trust this community to rack their brains over a solution. — loki, Jun 07 '17 at 07:37
You'd think that computationally it would be point in polygon bounded by gWithin() as possible solution set, which sounds like pick a point, any point, as long as it is within; though it looks like ArcMap does something of the sort. Perhaps pop it over to r-sig-geo — Chris, Jun 08 '17 at 14:03
After a deeper look into the ArcMap manual it seems that you are right: *"For an input polygon: the output point will be inside the polygon."*. Therefore, is there a function in R that generates a point at a random point within polygons? — loki, Jun 08 '17 at 14:25
st_sample(pkg sf) spsample(sp) [link](https://www.rdocumentation.org/packages/sf/versions/0.4-2/topics/st_sample) but still check results with rgeos gWithin() to drop boundary samples. Perhaps noded vectors: you have a restricted solution space (the polygon, which remains infinite as to possible points) and calculate the "center" as an equilibrium of vectors' thetas from the points setting up the polygon to their intersection within. Sounds good anyway. — Chris, Jun 08 '17 at 14:52

score 7 · Accepted Answer · edited Aug 21 '19 at 21:10

`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 \n"))

    return(centsDF)
  }

Super useful! Thanks for chasing and sharing - after so many years of trusting 'inside' option of Arc I was surprised to see that it is simply a random point! — radek, Jan 16 '18 at 07:29
@loki huge thank for your answer. However, there is a little problem for me when I use the `gCentroidWithin` function. It works well for the majority of polygons but there are those with no point inside and others with more than one point. It seems that the coordinates of centroids are not corrected for all polygons. I cannot find the problems. Could I send you the shapefile of my polygons ? — Basilique, May 17 '20 at 20:56
@Basilique have you considered using the sf version? I think you should give it a try. The sf package in general is a better option for working with geodata in R anyways. — loki, May 19 '20 at 19:44
Unfortunately I could not install this package since it required the package units and I did not manage to install it either despite all my effort. — Basilique, May 19 '20 at 23:41

Calculate Centroid WITHIN / INSIDE a SpatialPolygon

1 Answers1

`sf` solution

`sp` solution

Linked

Calculate Centroid WITHIN / INSIDE a SpatialPolygon

1 Answers1

sf solution

sp solution

Linked

`sf` solution

`sp` solution