2

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)

enter image description here

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)

enter image description here

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?

3

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
...
...
4
  • 1
    Thanks - this is what I expected, however your solution st_distance(df_point, polys_single will return the distances to every single pixel for each pixel value. I'd like the distance from the point to every patch of a given pixel value. I think what I'll do is break this up into a loop to evaluate each landscape separately as I'm currently doing this across all of my landscapes/site points simultaneously. – Jeremy Hemberger Jan 5 at 18:19
  • oh I understand, i'll update the code – Elio Diaz Jan 5 at 18:30
  • Thank you! I'm struggling to find any other way to turn individual patches into polygons. Even by looping through each landscape (as I suggest above) I run into the same issue where all patches of pixel value i are unified into one large polygon (even if they aren't connected). Hopefully this further clarifies what I'm after... – Jeremy Hemberger Jan 5 at 18:39
  • 1
    Your update worked @Elio - thank you! – Jeremy Hemberger Jan 5 at 18:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.