如何在R中围绕栅格中心旋转任意角度？（支持terra/stars等包）

问题描述 #

需在R语言中实现围绕栅格中心旋转任意角度（如32度），支持stars、terra等包，要求保持坐标系不变，仅旋转像元值。此前尝试stars仿射变换时因平移参数计算错误，导致旋转围绕左下角而非中心；gdalwarp也无法正确处理平移逻辑。

基于stars包的修正实现 #

以下代码修正了仿射变换的平移参数计算，确保旋转围绕栅格中心进行：

library(stars)
library(ggplot2)

# 加载并简化Landsat示例数据
l <- st_downsample(st_as_stars(L7_ETMs), 9) 
original_crs <- st_crs(l)

# 获取栅格基础参数
scale_x <- st_dimensions(l)$x$delta
scale_y <- st_dimensions(l)$y$delta
bbox <- st_bbox(l)

# 计算栅格中心地理坐标
centroid <- st_centroid(st_as_sfc(bbox)) |> st_coordinates()
cx <- centroid[1]
cy <- centroid[2]

# 设置旋转角度（可替换为任意数值）
theta_deg <- 32
theta_rad <- theta_deg * (pi / 180)

# 定义旋转+缩放的仿射参数
a <- scale_x * cos(theta_rad)
b <- scale_x * sin(theta_rad)
c <- -scale_y * sin(theta_rad)
d <- scale_y * cos(theta_rad)

# 计算中心对应的行列号（1-based）
col_c <- (cx - st_dimensions(l)$x$start) / scale_x + 1
row_c <- (cy - st_dimensions(l)$y$start) / scale_y + 1

# 计算平移参数：确保旋转后中心坐标与原中心一致
e <- cx - a * col_c - b * row_c
f <- cy - c * col_c - d * row_c

# 应用仿射变换并保持原坐标系
l_rotated <- l
st_geotransform(l_rotated) <- c(e, b, a, f, d, c)
st_crs(l_rotated) <- original_crs

# 可视化对比原栅格与旋转后栅格
ggplot() + 
  theme_bw() + 
  geom_stars(data = l, alpha = 0.8) +
  geom_stars(data = l_rotated, alpha = 0.4) +
  scale_fill_viridis_c() +
  labs(title = "原栅格（深色）与中心旋转32度栅格（浅色）")

基于terra包的实现 #

terra包通过自定义仿射变换矩阵也能实现相同效果，代码更简洁：

library(terra)

# 加载并简化数据
l <- rast(L7_ETMs)
l <- aggregate(l, 9) # 下采样减少计算量
original_crs <- crs(l)

# 获取栅格中心地理坐标与分辨率
cx <- xFromCol(l, ncol(l)/2 + 0.5) # 像素中心对应的x坐标
cy <- yFromRow(l, nrow(l)/2 + 0.5) # 像素中心对应的y坐标
res <- res(l)

# 设置旋转角度
theta_deg <- 32
theta_rad <- theta_deg * (pi / 180)

# 构建旋转+缩放的仿射矩阵参数
a <- res[1] * cos(theta_rad)
b <- res[1] * sin(theta_rad)
c <- -res[2] * sin(theta_rad)
d <- res[2] * cos(theta_rad)

# 计算中心对应的行列号（像素中心）
col_c <- ncol(l)/2 + 0.5
row_c <- nrow(l)/2 + 0.5

# 计算平移参数，保证中心位置不变
e <- cx - a * col_c - b * row_c
f <- cy - c * col_c - d * row_c

# 应用仿射变换并保持坐标系
l_rotated <- l
geom(l_rotated) <- c(a, b, c, d, e, f)
crs(l_rotated) <- original_crs

# 可视化对比
plot(l, alpha = 0.8, main = "原栅格与旋转后栅格")
plot(l_rotated, alpha = 0.4, add = TRUE)

核心逻辑说明 #

围绕中心旋转的关键在于平移参数的反向推导：

1. 先根据旋转角度计算缩放与旋转参数（a,b,c,d）；
2. 确定栅格中心对应的行列位置；
3. 通过原中心地理坐标反推平移参数（e,f），确保旋转后中心的地理坐标与原中心完全重合。

内容的提问来源于stack exchange，提问作者Thomas Moore