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.
2
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.
1
Version 2
of an approach mainly written in base R and {sf}, which iterates over the created ellipses and their centres (xy-coords), from lowest to hightest row number (id:= row number), to verify intersections.
Creating ellipses using {sf} and sfdep::st_ellipse()
coordinates.sf =
cbind(coordinates.sf, st_geometry(coordinates.sf) |>
lapply((i) sfdep::st_ellipse(geometry=i, sx=1400L, sy=1000L, rotation=225L, n=100L) |>
st_sf() |> st_cast("POLYGON")) |>
do.call(what="rbind") |>
`st_geometry<-`("geometryE") |>
`st_crs<-`(st_crs(coordinates.sf)))
Since my educational background is not programming, I am regularly lacking
an adequate choice of words for the approach I suggest (recommendations are more than welcome):
Iteratively (from lowest to highest id) identifying indices of ellipses which do not overlap any ellipses’s centres with lower ids:
Identification
idx = 1L
for(i in seq(nrow(coordinates.sf))[-1L]) {
q = length(unlist(with(coordinates.sf, st_intersects(geometry[idx], geometryE[i]))))
if(!q) idx = c(idx, i)
}
Plotting all together with {ggplot2}:
library(ggplot2)
ggplot(data=coordinates.sf[idx, ]) +
# geom_sf(aes(geometry=geometry), size=1L, shape=4L) +
geom_sf(aes(geometry=geometryE), fill="blue", alpha=.2) +
geom_sf_text(aes(label=id)) +
theme_minimal() +
labs(x="Easting", y="Northing")
giving
Note
(1) Input
Data given, re-shaped creation code
coordinates.df = data.frame(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))
coordinates.df$id = seq(nrow(coordinates.df)) # id
coordinates.df = coordinates.df[c("id", "x", "y")] # re-order
coordinates.sf = sf::st_as_sf(coordinates.df, coords = c("x", "y"), crs=25832) # create sf-object with crs
(2) References
This answer is based on a previous of mine.
3