r - Find direction that closest point of a polygon is from another polygon - Stack Overflow

admin管理员组

文章数量:1287595

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

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

• r
• gis
• polygon

Share Improve this question Follow asked Feb 25 at 7:26 lucianoluciano 13.9k3737 gold badges9494 silver badges134134 bronze badges 2

• 2 What are the points of your polygons? The five points from which the polygons are created or all points on the circumference of the newly created polygons? And from which perspective? – Friede Commented Feb 25 at 7:43
• All points on the circumference of the newly created polygons. I need the direction that the closest point of polygon_sf2 is from the closest point of polygon_sf1. – luciano Commented Mar 1 at 6:24

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 4

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)

本文标签： rFind direction that closest point of a polygon is from another polygonStack Overflow