code_01_ellipse_example_with_sample_data.R

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2

# Sample data
dat = data.frame(
  x = c(1, 1, -1, 3, 3),
  y = c(0, -3, 2, -2, 0),
  ex = c(0.5, 2, 2, 0.3, 0.6),
  ey = c(0.5, 0.2, 1, 1, 0.3),
  stringsAsFactors = FALSE
)
dat = st_as_sf(dat, coords = c("x", "y"))
dat

## Simple feature collection with 5 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -1 ymin: -3 xmax: 3 ymax: 2
## epsg (SRID):    NA
## proj4string:    NA
##    ex  ey     geometry
## 1 0.5 0.5  POINT (1 0)
## 2 2.0 0.2 POINT (1 -3)
## 3 2.0 1.0 POINT (-1 2)
## 4 0.3 1.0 POINT (3 -2)
## 5 0.6 0.3  POINT (3 0)

# Plot 1
plot(dat %>% st_geometry, graticule = TRUE, axes = TRUE, main = "Input")
text(dat %>% st_coordinates, as.character(1:nrow(dat)), pos = 2)

# Creating ellipse with x, y, ex, ey
st_ellipse = function(pnt, ex, ey, res = 30) {
  
  # Checks
  stopifnot(length(st_geometry(pnt)) == length(ex) & length(ex) == length(ey)) 
  stopifnot(class(ex) == "numeric" & class(ey) == "numeric") 
  stopifnot(class(st_geometry(dat))[1] == "sfc_POINT")
  
  # Point coordinates
  coords = st_coordinates(pnt)
  
  # Create ellipses
  result = list()
  for(i in 1:nrow(pnt)) {
    theta = seq(0, 2 * pi, length = res)
    x = coords[, 1][i] - ex[i] * cos(theta)
    y = coords[, 2][i] - ey[i] * sin(theta)
    elp = cbind(na.omit(x), na.omit(y))
    pnt = st_multipoint(elp)
    pol = st_cast(pnt, "POLYGON")
    result[[i]] = pol
  }

  # Combine  
  result = st_sfc(result)

  # Return result
  return(result)

}

# Calculate ellipses
el = st_ellipse(pnt = dat, ex = dat$ex, ey = dat$ey)

# Plot 2
plot(el, graticule = TRUE, axes = TRUE, main = "Output")
plot(dat %>% st_geometry, pch = 3, add = TRUE)
text(dat %>% st_coordinates, as.character(1:nrow(dat)), pos = 2)