CODEWITHSUNDEEP
rgispolygonr-sf
KNieder
KNieder

Reputation: 25

Select a random point along the longest edges of a polygon

I have a shapefile that consists of many polygons representing fields on a farm. For field sampling in real life, I need to walk a 50m transect along one of the 2 longest sides of each field (polygon). My original goal was to create a script that would place a random 50m line somewhere along the longest side of the polygon (like st_sample, but with a line instead of a point), but that seems to be much more difficult than just finding a point on the edge of a polygon.

What I have attempted: I haven't thought of a way to select a 50m line segment along the edge of the polygon, so my next thought was to use st_segmentize from the sf package, and then maybe try to filter for points that are 50m apart using st_buffer and use the two points as starting/ending points for the transect line. However, there is no guarantee that the points are on the longest sides of the polygon, nor that there will be any that are exactly 50m apart (given the nature of st_buffer, they could be more than 50m apart). This is what I've gotten so far:

library(sf)
library(tidyverse)

field <- st_polygon(x=list(cbind(x=c(8.069730,  8.069209, 8.069521, 8.070026, 8.069730),
                                 y=c(46.94751, 46.94925, 46.94922, 46.94754, 46.94751)))) %>% 
  st_sfc(crs = 4326)

edge_pts <- st_segmentize(field, dfMaxLength = 50) %>% 
    st_coordinates() %>% 
    as.data.frame() %>% 
    dplyr::select(X, Y) %>% 
    st_as_sf(coords = c("X", "Y")) %>% 
    st_set_crs(4326)

plot(field, reset = FALSE)
plot(edge_pts, reset = FALSE, add = TRUE)

enter image description here

Honestly, at this point, if I could even find a way to just get a single point near the middle of the two longest sides of the polygon (and not on the vertices of the field), I would be fine with that.

Upvotes: 0

Views: 135

Answers (2)

mrhellmann
mrhellmann

Reputation: 5529

The commented example below should split a polygon into segments, select the two longest segments, and return a point along each of the two segments.

library(sf)
library(tidyverse)
library(nngeo)

set.seed(6)                      # for reproducibility due to st_sample

field <- st_polygon(x=list(cbind(x=c(8.069730,  8.069209, 8.069521, 8.070026, 8.069730),
                                 y=c(46.94751, 46.94925, 46.94922, 46.94754, 46.94751)))) %>% 
  st_sfc(crs = 4326)

points <- 
field %>%
  st_as_sf() %>%
  st_segments() %>%                 # from nngeo package, splits polygon into lines
  mutate(len = st_length(.)) %>%    # length of each edge
  top_n(n = 2, wt = len) %>%        # select two longest sides from polygon
  st_sample(size = c(1,1)) 

#plotting:
ggplot() + 
  geom_sf(data = field) +
  geom_sf(data = points, color = 'red')

Created on 2023-05-25 by the reprex package (v2.0.1)

Upvotes: 2

user21508463
user21508463

Reputation:

As I am not R-speaking, I will just give you the maths.

Find the longest edge by computing Li:= (Xi+1 - Xi)² + (Yi+1 - Yi)² for all edges (the index wraps around for the last edge), and keeping the i that gives the largest value.

To get the length of this side, take the square root (so replace Li -> √Li).

If Li < 50, there is no solution. Otherwise, the coordinates are

X = Xi + 50/Li (Xi+1 - Xi)
Y = Yi + 50/Li (Yi+1 - Yi)

Addendum:

If your polygons are 3D, you simply add the interpolation on Z

Z = Zi + 50/Li (Zi+1 - Zi)

This assumes that the 50 m is an horizontal length. If you need the oblique length, use

Li:= (Xi+1 - Xi)² + (Yi+1 - Yi)² + (Zi+1 - Zi)²

instead.

Upvotes: 0

Related Questions