Rotate around the centroid in st_geotransform

[marine-ecologist]

I'm not sure whether this is an issue or just my misunderstanding of the st_geotransform() function. I'd like to use an affine transformation to translate and rotate a stars object around the centroid. As far as I understand, the general approach is to i) translate the object so that the centroid coincides with the origin, 2) perform the affine transformation, 3) translate it back.

Following these steps with an example stars object:

library(stars)
library(dplyr)
library(ggplot2)

# example with Landsat image from stars package
l <-  st_downsample(st_as_stars(L7_ETMs), 9) 

ggplot() + theme_bw() + 
  geom_stars(data=l, alpha=1) +
  scale_fill_viridis_c()

1. Translate the stars object so that the centroid is in the upper left corner:

# get the the step (pixel, cell) size for x and y dimensions
scale_x <- st_dimensions(l)$x$delta
scale_y <- st_dimensions(l)$y$delta

# get the translation from the lower left corner of the bbox
trans_x = st_bbox(l)[1] |> as.numeric()
trans_y = st_bbox(l)[4] |> as.numeric()

# set angle for rotation and convert to radians
theta_deg <- 0
theta_rad <- theta_deg * (pi / 180)

# define affine parameters for rotation around centroid
a = scale_x * cos(theta_rad)
b = scale_x * sin(theta_rad)     
c = -scale_y * sin(theta_rad)   
d = scale_y * cos(theta_rad)

# Calculate offset for the centroid given trans_x and trans_y
offset_x = (st_dimensions(l)$x$to * scale_x) / 2
offset_y = (st_dimensions(l)$y$to * scale_y) / 2

# Use the offsets in the affine transformation
e = trans_x - offset_x
f = trans_y - offset_y

# make new stars object and apply st_geotransform
l_affine <- l

st_geotransform(l_affine) = c(e, b, a, f, d, c)

l_affine_p1 <- l_affine |> st_transform(crs=st_crs(l)) |> 
  st_set_crs(st_crs(l)) 

#plot
ggplot() + theme_bw() + 
  geom_stars(data=l, alpha=0.2) +
  geom_stars(data=l_affine_p1, alpha=0.8) +
  scale_fill_viridis_c()

2. Rotate around the upper left corner:

theta_deg <- 45 # set angle for rotation
theta_rad <- theta_deg * (pi / 180) # convert to radians:

# define affine parameters for rotation around centroid
a = scale_x * cos(theta_rad)
b = scale_x * sin(theta_rad)     
c = -scale_y * sin(theta_rad)   
d = scale_y * cos(theta_rad)
e = trans_x - offset_x
f = trans_y - offset_y

# make new stars object and apply st_geotransform
l_affine2 <- l_affine1

st_geotransform(l_affine2) = c(e, b, a, f, d, c)

l_affine_p2 <- l_affine2 |> st_transform(crs=st_crs(l)) |> 
  st_set_crs(st_crs(l)) 

l_affine_w2 <- st_warp(l_affine2, crs=st_crs(l), cellsize=c(st_dimensions(l)$x$to,st_dimensions(l)$y$to))


#plot
ggplot() + theme_bw() + 
  geom_stars(data=l, alpha=0.2) +
  geom_stars(data=l_affine_p1, alpha=0.2) +
  geom_stars(data=l_affine_p2, alpha=0.8, na.rm=TRUE) +
  scale_fill_viridis_c(na.value="white")

3. translate the object back to remove the offset:

# set angle for rotation and convert to radians
theta_deg <- 45
theta_rad <- theta_deg * (pi / 180)

# define affine parameters for rotation around centroid
a = scale_x * cos(theta_rad)
b = scale_x * sin(theta_rad)     
c = -scale_y * sin(theta_rad)   
d = scale_y * cos(theta_rad)
e = trans_x + offset_x
f = trans_y + offset_y

# make new stars object and apply st_geotransform
l_affine_p3 <- l_affine_p2

st_geotransform(l_affine_p3) = c(e, b, a, f, d, c)

l_affine_p3 <- l_affine_p3 |> st_transform(crs=st_crs(l)) |> 
  st_set_crs(st_crs(l)) 

#plot
ggplot() + theme_bw() + 
  geom_stars(data=l, alpha=0.2) +
  geom_stars(data=l_affine_p1, alpha=0.2) +
  geom_stars(data=l_affine_p2, alpha=0.2) +
  geom_stars(data=l_affine_p3, alpha=0.8) +
  scale_fill_viridis_c()

The first step adds offset_x and offset_x, and the final step removes offset_x and offset_x, yet the rotated image is offset from the centroid, and the rotation is applied in both steps 2 & 3.

Is the "length 6 numeric vector" specified correctly here in st_geotransform and is there an approach to translate before rotating? Apologies in advance if this is more of a geometry than a stars specific issue but I can't seem to get st_geotransform to achieve this based on the documentation