4

I have a bunch of variously shaped and sized polygons with centroids. I want to calculate distance from each centroid to the furthest point of its respective polygon.

This question has been resolved here using the package::sp and package::rgeos.

According to its vignette, the sf package "aims at succeeding sp in the long term." Looking through the documentation I was unable to find a solution, but I'm no expert on simple features. Is there a good way to complete this operation using the sf package, or should I stick with sf and rgeos for now?

John J.
  • 1,450
  • 1
  • 13
  • 28

2 Answers2

6

Casting the polygon to POINT (thus getting the vertexes) and then computing distances wrt the centroid should work. Something like: 

library(sf)

# build a test poly
geometry <- st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,3),  c(0,0))))) 
pol <- st_sf(r = 5, geometry)

# compute distances 
distances <- pol %>% 
  st_cast("POINT") %>% 
  st_distance(st_centroid(pol))

distances
#>          [,1]
#> [1,] 1.201850
#> [2,] 1.054093
#> [3,] 2.027588
#> [4,] 1.201850

# maximum dist:
max_dist <- max(distances)
max_dist
#> [1] 2.027588

# plot to see if is this correct: seems so.
plot(st_geometry(pol))
plot(st_centroid(pol), add = T)
plot(st_cast(pol, "POINT")[which.max(distances),],
     cex =3, add = T, col = "red")

You get the same distance twice, since the first and last vertex are the same, but since you are interested in the maximum it should not matter.

HTH

lbusett
  • 5,801
  • 2
  • 24
  • 47
  • I didn't knew `st_cast`. Interesting ! However when I try to use it in my answer, it works correctly for one polygon but not in the function `furthest` for several polygons. It seems that this is because the `lapply` part at the end does : `st_cast(st_geometry(poly)[1], "POINT")` which gives the coordinates of the first point only instead of `st_cast(st_geometry(poly)[[1]], "POINT")` which gives the expected coordinates of all vertices. – Gilles San Martin Jan 25 '18 at 23:35
1

I'm not certain to understand what you want exactly : the distance to the furthest point (this what you ask) or the coordinates of the furthest point (this is what is provided by the answer you point to).

Here is a solution to calculate the distance (and that can easily be changed to extract the coordinates)

# This is an example for one polygon.
# NB the polygon is the same as in the answer pointed to in the question

library(sf)
sfPol <- st_sf(st_sfc(st_polygon(list(cbind(c(5,4,2,5),c(2,3,2,2))))))

center <- st_centroid(sfPol)
vertices <-  st_coordinates(sfPol)[,1:2]
vertices <-  st_as_sf(as.data.frame(vertices), coords = c("X", "Y"))
furthest <- max(st_distance(center, vertices))
furthest

## [1] 1.699673



# You can adapt this code into a function that will work 
# for multiple polygons

furthest <- function(poly) {
    # tmpfun find the furthest point from the centroid of one unique polygon
    tmpfun <- function(x) {
        center <- st_centroid(x)
        vertices <-  st_coordinates(x)[,1:2]
        vertices <-  st_as_sf(as.data.frame(vertices), coords = c("X", "Y"))
        furthest <- max(st_distance(center, vertices))
        return(furthest)
    }

    # We apply tmpfun on each polygon
    return(lapply(st_geometry(poly), tmpfun))
}


poly1 <- cbind(c(5,4,2,5),c(2,3,2,2))
poly2 <- cbind(c(15,10,8,15),c(5,4,12,5))

sfPol <- st_sf(st_sfc(list(st_polygon(list(poly1)), 
                           st_polygon(list(poly2)))))

furthest(sfPol)

## [[1]]
## [1] 1.699673
## 
## [[2]]
## [1] 5.830952
Gilles San Martin
  • 4,224
  • 1
  • 18
  • 31