Here is a code that creates ellipses of 1400 metres length in the main axis and 1000 metres length in the side axis. Each ellipse is assigned with an ID
.
library(sf)
library(ggplot2)
library(dplyr)
x <- c(611547.6411, 589547.6411, 611447.6411, 609847.6411, 606347.6411, 611447.6411, 613547.6411,642747.6411, 589647.6411, 606447.6411, 613547.6411, 640347.6411, 642847.6411, 612147.6411, 613847.6411, 640247.6411, 642947.6411, 584347.6411, 587747.6411, 606447.6411, 614247.6411, 640447.6411, 642747.6411, 584447.6411, 608647.6411, 612047.6411, 612747.6411,
613847.6411, 643147.6411, 583147.6411, 608747.6411, 611847.6411, 609647.6411, 610047.6411, 613747.6411, 586247.6411, 588647.6411, 643147.6411, 584347.6411, 606447.6411, 610147.6411, 613347.6411, 614647.6411, 586047.6411, 587247.6411, 611547.6411, 640347.6411, 643147.6411, 587147.6411, 583047.6411, 608747.6411, 612047.6411, 613947.6411, 587647.6411, 588547.6411, 586847.6411, 611247.6411, 643247.6411, 587247.6411, 590347.6411, 582947.6411, 608947.6411, 611847.6411, 613447.6411, 614647.6411, 585147.6411, 587647.6411, 588547.6411, 586947.6411, 611247.6411, 643047.6411, 587147.6411, 583947.6411, 587747.6411, 608547.6411, 611747.6411, 614047.6411, 585247.6411, 586247.6411, 588447.6411, 589147.6411, 611347.6411, 642447.6411, 586947.6411, 585847.6411, 587747.6411, 581447.6411, 612447.6411, 611947.6411, 600547.6411,
612047.6411, 610347.6411, 614147.6411, 582847.6411, 588547.6411, 589247.6411, 611247.6411, 638147.6411, 640547.6411, 642947.6411, 587047.6411, 585947.6411, 587647.6411, 600447.6411, 611347.6411, 612347.6411, 610347.6411, 587747.6411, 579747.6411, 583847.6411, 586847.6411, 588447.6411, 589347.6411, 643347.6411, 589347.6411, 586947.6411, 588247.6411, 588847.6411, 585847.6411, 590847.6411, 589447.6411, 590947.6411, 581347.6411, 611847.6411, 600647.6411, 610347.6411, 615947.6411, 613947.6411, 586347.6411, 579647.6411, 584047.6411, 586347.6411, 587747.6411, 587947.6411, 586547.6411, 587647.6411, 614047.6411, 643047.6411, 587947.6411, 585747.6411, 584947.6411, 600547.6411, 611947.6411, 606847.6411, 600847.6411, 612847.6411, 615747.6411, 620747.6411, 614047.6411, 632947.6411, 588147.6411, 579747.6411, 582747.6411)
y <- c(5272140.5728, 5271740.5728, 5271640.5728, 5267440.5728, 5271540.5728, 5272040.5728, 5272340.5728, 5268540.5728, 5271240.5728, 5271640.5728, 5272140.5728, 5272240.5728, 5272240.5728, 5277940.5728, 5278040.5728, 5278040.5728, 5266940.5728, 5267040.5728, 5267440.5728, 5268140.5728, 5268640.5728, 5271140.5728, 5271740.5728, 5271740.5728, 5271940.5728, 5272140.5728, 5272240.5728, 5272040.5728, 5272140.5728, 5272140.5728, 5272140.5728, 5272240.5728, 5272340.5728, 5277240.5728, 5278040.5728, 5268540.5728, 5271240.5728, 5271340.5728, 5272240.5728, 5271940.5728, 5271940.5728, 5272040.5728, 5272040.5728, 5272040.5728, 5272040.5728, 5272040.5728, 5272140.5728, 5272140.5728,
5272140.5728, 5272240.5728, 5272240.5728, 5277240.5728, 5278040.5728, 5278140.5728, 5278140.5728, 5265540.5728, 5266840.5728, 5266940.5728, 5267040.5728, 5268540.5728, 5272240.5728, 5272340.5728, 5272040.5728, 5272040.5728, 5277340.5728, 5278140.5728, 5278140.5728, 5265640.5728, 5266840.5728, 5267240.5728, 5268440.5728, 5271540.5728, 5272140.5728, 5271840.5728, 5271940.5728, 5271940.5728, 5271940.5728, 5272040.5728, 5272040.5728, 5272140.5728, 5272140.5728, 5272140.5728, 5272340.5728, 5277140.5728, 5277240.5728, 5277340.5728, 5277740.5728, 5277740.5728, 5278040.5728, 5278140.5728, 5278140.5728, 5278240.5728, 5278240.5728, 5264940.5728, 5265040.5728, 5265140.5728, 5266740.5728, 5266840.5728, 5266940.5728, 5267040.5728, 5267140.5728, 5267340.5728, 5267440.5728, 5268340.5728,
5271240.5728, 5271840.5728, 5271940.5728, 5272040.5728, 5272040.5728, 5272340.5728, 5271840.5728, 5271840.5728, 5272140.5728, 5272140.5728, 5272240.5728, 5272240.5728, 5272340.5728, 5274340.5728, 5274440.5728, 5274640.5728, 5285140.5728, 5285240.5728, 5277340.5728, 5277540.5728, 5277840.5728, 5278040.5728, 5278040.5728, 5278140.5728, 5278140.5728, 5265540.5728, 5265640.5728, 5266740.5728, 5266740.5728, 5266940.5728, 5268340.5728, 5268440.5728, 5271440.5728, 5271540.5728, 5271540.5728, 5271740.5728, 5272040.5728, 5272340.5728, 5271740.5728, 5272240.5728, 5272240.5728, 5274540.5728, 5275040.5728, 5275340.5728, 5284840.5728, 5284940.5728, 5284940.5728, 5285040.5728, 5285040.5728)
# create data frame
coordinates.df <- data.frame(x = x, y = y)
# add ID column
coordinates.df$ID <- 1:nrow(coordinates.df)
coordinates.df <- coordinates.df[c(3,1:2)]
# convert data frame to sf object
coordinates.sf <- st_as_sf(coordinates.df, coords = c("x", "y"), crs = 25832)
# function for creating ellipses
create_ellipse <- function(center, a = 1400, b = 1000, angle = 225, n = 100) {
angle_rad <- angle * pi / 180
angles <- seq(0, 2 * pi, length.out = n)
ellipse_coords <- cbind(a * cos(angles), b * sin(angles))
rotation_matrix <- matrix(c(cos(angle_rad), -sin(angle_rad),
sin(angle_rad), cos(angle_rad)),
nrow = 2)
rotated_ellipse <- as.matrix(ellipse_coords) %*% rotation_matrix
x_center <- center[1]
y_center <- center[2]
rotated_ellipse <- rotated_ellipse + matrix(c(x_center, y_center), nrow = n, ncol = 2, byrow = TRUE)
rotated_ellipse <- rbind(rotated_ellipse, rotated_ellipse[1, ])
st_polygon(list(rotated_ellipse))
}
# aplly on coordinates
ellipses <- st_coordinates(coordinates.sf) %>%
apply(1, function(p) create_ellipse(p)) %>%
st_sfc(crs = st_crs(coordinates.sf))
# convert ellipses into sf objects
ellipses.sf <- st_sf(geometry = ellipses)
# plot
ggplot() +
geom_sf(data = coordinates.sf, color = "black", size = 2) +
geom_sf(data = ellipses.sf, fill = "blue", alpha = 0.2) +
theme_minimal() +
labs(x = "Easting", y = "Northing")
resulting in
There are overlapping ellipses. Overlaps are ok as long as the center of an ellipse is not overlapped:
In the first case the upper ellipse is not overlapping the center of the lower one, which is fine. In this case both ellipses should be kept. In the second case, the upper ellipse is overlapping the lower one including its center, which is not ok. In this case, the ellipse with the higher ID
should be removed, i.e. the ellipse with the lower ID
remains. Ellipses that do not touch another ellipses should also be retained.
I think this is what you're after. I'm using data.table
out of habit but could be adapted to dplyr
or base. In essence, this just uses sf::st_intersects()
but turns it into something more usable IMO.
library(data.table)
# intersects
i <- sf::st_intersects(coordinates.sf, ellipses.sf)
i <- as.matrix(i)
rownames(i) <- seq_len(nrow(i))
colnames(i) <- seq_len(ncol(i))
i <- i |>
as.data.frame.table(stringsAsFactors = FALSE) |>
data.table::as.data.table() |>
data.table::setnames(c("I1", "I2", "Intersect")) |>
_[Intersect == TRUE] |>
_[, c("I1", "I2") := lapply(.SD, as.integer), .SDcols = c("I1", "I2")] |>
_[order(I1, I2)]
# remove identity
i <- i[I1 != I2]
# ensure that I2 > I1
# if 1 intersects 3 then 3 intersects 1
# this makes sure that only polygon 3 is removed, keeping polygon 1
i <- i[I2 > I1]
# now subset ellipses
ellipses.sf.subset <- ellipses.sf[!seq_len(nrow(ellipses.sf)) %in% i$I2, ]
# plot
ggplot() +
geom_sf(data = coordinates.sf, color = "black", size = 2) +
geom_sf(data = ellipses.sf.subset, fill = "blue", alpha = 0.2) +
theme_minimal() +
labs(x = "Easting", y = "Northing")
Using mapview::mapview
to verify the suspicious looking ones circled in red.
Edit
Upon further investigation, that doesn't quite work out. i[I2 > I1]
removes a few ellipses that it shouldn't. I knew that was too easy... It fails in the case where an ellipse is removed by an ellipse that would've previously been removed. The code is now much nastier (and could probably be cleaned up substantially) but I believe this produces the actual desired result.
# ensure that I2 > I1
# if 1 intersects 3 then 3 intersects 1
# this makes sure that only polygon 3 is removed, keeping 1
rm <- c()
for (ix in unique(i$I1)) {
x <- i[I1 == ix & I2 > I1 & !I1 %in% rm, I2]
x <- x[!x %in% rm]
if (length(x) > 0) {
rm <- c(rm, x)
}
}
# now subset ellipses
ellipses.sf.subset <- ellipses.sf[!seq_len(nrow(ellipses.sf)) %in% rm, ]
The red circles here indicate the ellipses that should not have been removed by the first approach, but were.