sf解决方案

随着sf软件包的出现，事情变得简单了一些.只需使用:

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)

返回保证在(多)曲面上的点."

正如Question更新中所指出的那样，ArcMap似乎只是将点放置在多边形内的任意位置.这也可以通过gPointsOnSurface(..., n = 1, type = 'random')来实现.

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