CODEWITHSUNDEEP
rgisr-sf
Philippe Massicotte
Philippe Massicotte

Reputation: 1529

How to calculate the distance to the closest polygon in R with the sf package

I have a point for which I would like to calculate the distance to the closest polygon.

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(rnaturalearth)
library(ggplot2)

pts <- st_as_sf(
    data.frame(longitude = 3.7, latitude = 52.3667),
    coords = c("longitude", "latitude"),
    crs = 4326
  )

ne_land <- ne_download(
  category = "physical",
  type = "land",
  returnclass = "sf",
  scale = "small"
)
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmplq0wPd", layer: "ne_110m_land"
#> with 127 features
#> It has 3 fields

This graph shows where the point is located in relation to the land polygon.

ggplot() +
  geom_sf(data = ne_land) +
  geom_sf(data = pts) +
  coord_sf(xlim = c(2, 5), ylim = c(50, 54))

Then I am trying to get the distance to the closest polygon with:

nearest <- st_nearest_feature(pts, ne_land)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
nearest
#> [1] 113

dist <-  st_distance(pts, ne_land[nearest, ], by_element=TRUE)
dist
#> 0 [m]

As seen, the dist is equal to 0 meter. Based on the previous graph, this is not the result I was expecting. Any help is appreciated. EDITS: Based on the last answer, I have tried with some other locations and higher land resolution and also get 0 meters distance.

ne_land <-
  rnaturalearth::ne_download(
    category = "physical",
    type = "land",
    returnclass = "sf",
    scale = "large"
  )
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmplq0wPd", layer: "ne_10m_land"
#> with 10 features
#> It has 3 fields

pts <- st_as_sf(
  data.frame(longitude = 4.1833, latitude = 52.1833),
  coords = c("longitude", "latitude"),
  crs = 4326
)

st_distance(pts, st_union(ne_land))
#> although coordinates are longitude/latitude, st_union assumes that they are planar
#> Units: [m]
#>      [,1]
#> [1,]    0

ggplot() + 
  geom_sf(data = st_nearest_points(pts, st_union(ne_land))) + 
  geom_sf(data = ne_land) +
  geom_sf(data = pts) +
  coord_sf(xlim = c(2, 5), ylim = c(50, 54))
#> although coordinates are longitude/latitude, st_union assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar

Created on 2021-04-06 by the reprex package (v2.0.0)

Upvotes: 1

Views: 2357

Answers (1)

mrhellmann
mrhellmann

Reputation: 5489

You need to use st_union() on the ne_land data to find the single closest point.

st_distance(pts, st_union(ne_land))

although coordinates are longitude/latitude, st_union assumes that they are planar
Units: [m]
         [,1]
[1,] 36319.29

ggplot() + 
 geom_sf(data = st_nearest_points(pts, st_union(ne_land))) + 
 geom_sf(data = ne_land) +
 geom_sf(data = pts) +
 coord_sf(xlim = c(2, 5), ylim = c(50, 54))

enter image description here

Edit: The following distances are much different from the above distance since pts has changed as in the edited question.

I don't know why it isn't working with the large naturalearth data, but here is a workaround:

st_nearest_points(pts, ne_land) %>% 
  st_length()

#although coordinates are longitude/latitude, st_nearest_points assumes #that they are planar
#11481.03 [m]

edit2: This also works, but returns a different distance. This probably has to do with the projection (crs) used or the simplification of polygons to lines.

st_combine(ne_land) %>% st_cast('MULTILINESTRING') %>% st_distance(pts)
Units: [m]
         [,1]
[1,] 10136.87

Upvotes: 2

Related Questions