Search code examples
r

Draw ellipses around a point and select an ellipse if they overlap


Here is some code to reproduce a data.frame containing projected coordinates and to add an ID column:

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 the data frame
coordinates.df <- as.data.frame(cbind(x,y))

# add an ID column and rearrange the columns
coordinates.df$ID <- 1:nrow(coordinates.df)
coordinates.df <- coordinates.df[c(3,1:2)]

Each pair of x and y is a geographic point. I would like to accomplish:

  1. Draw an ellipse around each point. The main axis is supposed to be 700 m long, the minor axis 500 m. The orientation of all ellipses should be 225 degrees.

  2. There will be overlaps of some ellipses. In case of overlaps, the ellipe with the lower or lowest ID value should be kept, while the others should be removed. So at the end there should be no overlapping ellipses.

EDIT: The Projected Coordinate System is ETRS 1989 UTM Zone 32N.


Solution

  • sfdep::st_ellipse() and a simple iteration should help here.

    Create ellipses

    coordinates.sf$E = sf::st_geometry(coordinates.sf) |>
      lapply(\(i) sfdep::st_ellipse(geometry=i, sx=1400L, sy=1000L, rotation=225L, n=100L) |> 
               sf::st_sf() |> sf::st_cast("POLYGON")) |> 
      do.call(what="rbind") |> 
      sf::`st_geometry<-`("geometryE") |>
      `st_crs<-`(st_crs(coordinates.sf))
    

    and loop over ellipses, starting with row number 1

    idx = 1L 
    for(i in seq(nrow(coordinates.sf))[-1L]) {
      q = length(unlist(sf::st_intersects(coordinates.sf$E$geometry[idx], coordinates.sf$E$geometry[i])))
      if(!q) idx = c(idx, i)
    }
    

    Plot

    library(ggplot2)
    ggplot(data=coordinates.sf[idx, ]) +
      geom_sf(aes(geometry=E$geometryE), fill="blue", alpha=.2) +
      geom_sf_text(aes(label=id)) + 
      theme_minimal() +
      labs(x="Easting", y="Northing")
    

    gives

    enter image description here


    Input Data

    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