Rotate a raster an arbitrary angle in R

Question

I'd like to rotate a raster by an arbitrary angle (say... 32 degrees) around the centre of the raster. Using any of the R packages (terra, raster, stars etc). How can I approach a general solution to rotating rasters? Example below:

library(stars)
library(terra)
library(raster)

# stars
l <-  st_downsample(st_as_stars(L7_ETMs), 9) 
plot(l)

# terra
l <- rast(l)
plot(l)

#raster
l <- raster(l)
plot(l)

I've tried affine transformations in stars, but I can't get the rotation around the centre to work correctly. gdalutils offers gdalwarp but again, I can't seem to get the translation of the affine to function correctly.

edit:

I'd like to rotate the cell values while keeping coordinate system constant.

The following example in stars uses affine transformation to rotate around the lower left corner (while maintaining the CRS), but I can't figure out how to correctly use translate in the affine transformation to offset the rotation around the centre of the raster.

enter image description here

library(stars)
library(tidyverse)
library(ggplot2)

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

# 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)[2] |> as.numeric()

# get centroid
bbox <- st_bbox(l)
sfc_bbox <- st_as_sfc(bbox)
centroid <- st_centroid(sfc_bbox) |> st_coordinates()

# set angle for rotation
theta_deg <- 45

# convert to radians:
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 = centroid[1] - trans_x
offset_y = centroid[2] + trans_y

# Use the offsets in the affine transformation
e = trans_x #+ offset_x - (a * centroid[1] + b * centroid[2])
f = trans_y #+ offset_y - (c * centroid[1] + d * centroid[2])


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

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

# convert to crs
l_affine <- l_affine |> st_transform(crs=st_crs(31985)) |> 
  st_set_crs(31985) 

#plot
ggplot() + theme_bw() + 
  geom_stars(data=l, alpha=0.8) +
  geom_stars(data=l_affine, alpha=0.4) +
  scale_fill_viridis_c()