Search code examples
r

Retain overlapping ellipses if the center is not overlapped


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

enter image description here

There are overlapping ellipses. Overlaps are ok as long as the center of an ellipse is not overlapped: enter image description here

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.


Solution

  • 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")
    

    ggplot2

    Using mapview::mapview to verify the suspicious looking ones circled in red.

    verification

    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.

    approach 2