Pivot sf geometry column such that each row contains a single polygon

Question

I am trying to calculate the average distance to all habitat patches within a landscape from a site point (sf object). To do this, I have polygonized my habitat raster (dissolving to get patches) and am hoping to use st_distance() to calculate the distance from my points to each habitat patch (with 4 possible habitat types).

When I polygonize my raster and convert to an sf object, the resulting data frame looks like this:

habitats.sf <- rasterToPolygons(habitat.raster, dissolve = TRUE) %>%
  st_as_sf()
view(habitats.sf)

Each layer has a number of geometries corresponding to each patch within the layer.

When I calculate the distance between my points and these habitat polygons, the result looks like this:

st_distance(sites.sf, habitats.sf, by_element = FALSE)

with each matrix row corresponding to a site, and each column corresponding to a habitat type.

Is the st_distance() calculation calculating the distance to the closest polygon of each class, or the mean distance of all of the polygons of a given class? If the former, how can I pivot my habitats.sf data frame such that each row only contains a single polygon geometry rather than a list so that I can calculate the distance to ALL patch geometries within each habitat type?

Answer 1

You are passing dissolve = TRUE to the rasterToPolygons function, which means all pixels of a same value will be merged into a multipolygon; as you asked, st_distance returns the shortest distance.

In the following example a projected CRS is used, although it is noteworthy (an really useful) that sf will return an ellipsoidal distance in meters if CRS is lat long

library(raster)
library(sf)
library(dplyr)

r = raster(xmn = 500000, xmx = 505000, ymn = 2000000, ymx = 2005000, ncol = 20, nrow = 20)
values(r) = sample(1:5, 10*10, replace = T)

polys_diss = raster::rasterToPolygons(r, dissolve = T) %>% st_as_sf() %>%
  st_set_crs(32615)
polys_single = raster::rasterToPolygons(r, dissolve = F) %>% st_as_sf() %>%
  st_set_crs(32615)

df_point = data.frame(wkt = c("POINT (503250 2001250)","POINT (501250 2001200)" ) ) %>% st_as_sf(wkt = "wkt") %>%
  st_set_crs(32615)

st_distance(df_point, polys_single)
Units: [m]
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]
[1,] 4257.347 3482.097 3250.000 2850.439 3181.981 2371.708 2371.708 2573.908 1903.943
[2,] 3384.154 3309.456 3735.305 3924.602 2313.547 2419.194 3585.038 3981.520 1950.000
        [,10]    [,11]    [,12]    [,13]     [,14]     [,15] [,16] [,17] [,18] [,19]
[1,] 3020.761 2150.581 1274.755 1767.767 1274.7549  353.5534  2750  2250   250     0
[2,] 1500.833 1300.000 2598.557 3500.357  390.5125 1285.4960   750   250  1250  1750
         [,20]     [,21]
[1,]  353.5534 1457.7380
[2,] 2258.8714  743.3034

And the dissolved layer:

st_distance(df_point, polys_diss[polys_diss$layer == 3,])
Units: [m]
      [,1]
[1,]    0
[2,]  250

EDIT

As you mentioned, st_distance will give you the closest part of the 5 multipolygons and you are looking for the distance to each patch, this can be done with a call to st_cast:

# the outcome from `rasterToPolygons` is a MULTIPOLYGON object consisting of
# all the same class polygons
polys_simple = polys_diss %>% st_cast("POLYGON") 
st_distance(df_point, polys_simple[polys_simple$layer == 3,])

Units: [m]
         [,1]     [,2]     [,3]     [,4] [,5]     [,6]     [,7]      [,8]     [,9]
[1,] 1677.051 1250.000 1520.691 2136.001 1500  500.000 3010.399  707.1068  500.000
[2,] 3569.314 2909.467 3512.834  800.000   50 1001.249 1044.031 2559.7851 2074.247
         [,10] [,11]     [,12]     [,13]    [,14]    [,15]     [,16] [,17]    [,18]
[1,] 2549.5098     0  353.5534 2371.7082 3041.381  559.017 1581.1388   500 1767.767
[2,]  743.3034  1750 1328.5330  538.5165 1096.586 2294.559  514.7815  2500  200.000
...
...