Retain overlapping ellipses if the center is not overlapped

Asked 12/9, 2024 at 15:45 Answered 13/9, 2024 at 11:16

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.

Archaize answered 12/9, 2024 at 15:45 Comment(2)

What if in the figure with the green check, ID6 were rotated so that it overlapped the center of ID4, but ID4 did not overlap the center of ID6 (or vice versa)? – Calia 12/9, 2024 at 16:1

Rotation is not allowed in this case, all ellipses must have the same orientation. Generally, if one center is overlapped, the ellipse with the higher ID needs to be removed. – Archaize 12/9, 2024 at 16:7

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.

Whomp answered 12/9, 2024 at 17:20 Comment(1)

Thanks a lot sonshine, it appears that it does exactly what was desired. – Archaize 13/9, 2024 at 8:7

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.

Bowler answered 13/9, 2024 at 11:16 Comment(3)

Thanks again Friede for the answer. That does go into the right direction, but in direct comparsion to the answer of sonshine above, a few too many ellipses are deleted here. 53 ellipses remain here while with the solution above there are 59 left. – Archaize 13/9, 2024 at 11:31

Great, now it delivers the exact same result as above. – Archaize 13/9, 2024 at 13:18

Not too satisfying. I do not find the vectorised alternative. Seems like I am overlooking something. – Bowler 13/9, 2024 at 17:32

Hot tags

Godot Unity Godot Help Programming Godot 4.X GUI GDScript 3D 2D Physics CSharp Godot 3.X VR XR Projects C++

Creating ellipses using {sf} and sfdep::st_ellipse()