I have a data objects that contains a sample of measurements obtained for two air pollutants at a series of monitoring stations. I also have the spatial coordinates of these stations. For example:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
library(ggplot2)

monitoring_stations <- st_sf(
  data.frame(ID = c("A", "B", "C")),
  geometry = st_sfc(
    st_point(c(1, 1)),
    st_point(c(2, 2)),
    st_point(c(3, 3))
  )
)
values <- data.frame(
  ID = sample(c("A", "B", "C"), size = 100, replace = TRUE), 
  value1 = runif(100), 
  value2 = runif(100)
)
(full_data <- merge(monitoring_stations, values, by = "ID"))
#> Simple feature collection with 100 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 3 ymax: 3
#> CRS:           NA
#> First 10 features:
#>    ID     value1     value2    geometry
#> 1   A 0.13221654 0.76275479 POINT (1 1)
#> 2   A 0.67336289 0.29259564 POINT (1 1)
#> 3   A 0.25384031 0.35145494 POINT (1 1)
#> 4   A 0.99337433 0.30237615 POINT (1 1)
#> 5   A 0.68685394 0.42133826 POINT (1 1)

I can create a plot displaying the data ellipse for each monitoring station,
for example

ggplot(full_data, aes(value1, value2)) + 
  geom_point() + 
  stat_ellipse(aes(group = ID, col = ID))

and I can also display the spatial locations of the monitoring stations

ggplot(full_data) + 
  geom_sf()

but I would like to merge the two previous plots creating a map that
represent 1) the monitoring stations (as dots); 2) the data ellipse at each
monitoring station overlapped on that map. How can I accomplish this task using R and ggplot2 or any other package?

^{Created on 2024-05-14 with reprex v2.0.2}