CODEWITHSUNDEEP
rgispolygon
luciano
luciano

Reputation: 13850

Find direction that closest point of a polygon is from another polygon

How do I find the direction that the closest point of polygon_sf2 is from polygon_sf1 divided into one of eight directions (north, north-east, east, south-east, south, south-west, west, north-west)? In this case it is either north-east or east.

library(sf)

# make polygons
coords <- matrix(c(
  -1.6522, 55.571,
  -1.6487, 55.568,
  -1.6522, 55.565,
  -1.6550, 55.568,
  -1.6522, 55.571
), ncol = 2, byrow = TRUE)

polygon <- st_sfc(st_polygon(list(coords)), crs = 4326)

polygon_sf1 <- st_sf(geometry = polygon)

coords <- matrix(c(
  -1.645, 55.575,
  -1.6449, 55.571,
  -1.6472, 55.569,
  -1.649, 55.571,
  -1.645, 55.575
), ncol = 2, byrow = TRUE)

polygon2 <- st_sfc(st_polygon(list(coords)), crs = 4326)

polygon_sf2 <- st_sf(geometry = polygon2)

# closest points between two polygons
st_distance(polygon_sf1, polygon_sf2)

# I need the direction that the closest point of polygon_sf2 is from polygon_sf1 divided into one of eight directions (north, north-east, east, south-east, south, south-west, west, north-west). In this case it is either north-east or east

Upvotes: 3

Views: 114

Answers (1)

Friede
Friede

Reputation: 7978

Initialise

sf-objects as

 library(sf)
 p1 = st_as_sf(as.data.frame(coords1), coords = c('V1', 'V2'), crs = 4326) 
 polygon_sf1 = p1 |> st_combine() |> st_cast('POLYGON')
 p2 = st_as_sf(as.data.frame(coords2), coords = c('V1', 'V2'), crs = 4326) 
 polygon_sf2 = p2 |> st_combine() |> st_cast('POLYGON')

We might need p1, p2 later on.

Functions

Lets borrow orient() from this answer.

 orient = \(lines) {
  lines |>
    sf::st_coordinates() |>
    data.frame() |>
    split(~L1) |>
    sapply(\(p) geosphere::bearing(p[1, 1:2], p[nrow(p), 1:2]))
  }

Notice that geosphere::bearing() returns azimuths (angles) expressed as degrees between -180, ..., 180. That is the relative direction you would turn when facing "[N]orth".

To return cardinal direction we can create a wrapper ang2cd().

ang2cd = \(a) {
  a = (a + 360) %% 360
  cd = c('north', 'north-east', 'east',' south-east', 
         'south','south-west', 'west', 'north-west', 'north')
  d = seq.int(0, 360, 22.5)[c(FALSE, TRUE)]
  cd[findInterval(a , d, rightmost.closed = TRUE) + 1]
}

Apply

> o = st_nearest_points(polygon_sf1, polygon_sf2)
> ang2cd(orient(o))
[1] "south-west"
> os = st_nearest_points(polygon_sf2, polygon_sf1)
> ang2cd(orient(os))
[1] "north-east"

(Fine-tuning. The set-up of cd and d in ang2cd() is left to you.)

Plot

library(ggplot2)
ggplot() +
  geom_sf(data = polygon_sf2) +
  geom_sf(data = polygon_sf1) +
  geom_sf(data = o) +
  ggspatial::annotation_north_arrow(
    location = 'tl', which_north = 'true',
    style = ggspatial::north_arrow_nautical(
      fill = c('grey40', 'white'),
      line_col = 'grey20')) +
  theme_bw()

Given typically north-oriented plots, it seems more logical to use 'south-west' instead of the 'north-east' or 'east' you requested. However, which answer is more appropriate depends on a person's perspective represented by the order of arguments to st_nearest_points(). We can read those as from point1 to point2, see demonstration above. A standardised approach would be to use directions expressed in degrees from 0, ..., 360.

Data

 coords1 = matrix(c(
   -1.6522, 55.571,
   -1.6487, 55.568,
   -1.6522, 55.565,
   -1.6550, 55.568,
   -1.6522, 55.571
 ), ncol = 2, byrow = TRUE)

 coords2 = matrix(c(
   -1.645, 55.575,
   -1.6449, 55.571,
   -1.6472, 55.569,
   -1.649, 55.571,
   -1.645, 55.575
 ), ncol = 2, byrow = TRUE)

Upvotes: 4

Related Questions