如何反转st_transform，将变换后CRS坐标还原为真实地理坐标？

问题背景 #

我通过带gamma参数的proj4string实现了sf对象的旋转，在ggplot2和plotly中渲染效果正常，但在Shiny应用里，用户点击返回的是变换后CRS的坐标，而非原始的WGS84经纬度（如示例中的-79.8524, 36.634），需要反转该变换还原原始坐标，要求理论准确，近似精度即可。

用到的旋转CRS与示例数据 #

# 原始点数据
realpoint <- structure(list(lon = -79.8524, lat = 36.634), row.names = c(NA, -1L), class = "data.frame")
realpoint_sf <- structure(list(geometry = structure(list(structure(c(-79.8524, 36.634), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = -79.8524, ymin = 36.634, xmax = -79.8524, ymax = 36.634), class = "bbox"))), row.names = 1L, sf_column = "geometry", agr = structure(integer(0), class = "factor", levels = c("constant", "aggregate", "identity"), names = character(0)), class = c("sf", "data.frame"))

# 旋转CRS定义
crs_rot <- " + proj=omerc + lat_0=36.634073240 + lonc=-79.851657820 + datum=WGS84 + alpha=0 + gamma=67.000"

# 正向变换示例
realpoint_sf_rot <- sf::st_transform(realpoint_sf, crs = crs_rot)

逆变换解决方案 #

核心是将点击得到的变换后坐标标记为旋转CRS，再转换回原始WGS84坐标系：

# 假设userclick是Shiny中用户点击返回的变换后坐标数据框
userclick <- data.frame(sf::st_coordinates(realpoint_sf_rot))

# 1. 将点击坐标转为带旋转CRS的sf对象
clicked_sf_rot <- sf::st_as_sf(userclick, coords = c("X", "Y"), crs = crs_rot)

# 2. 转换回原始WGS84坐标系（EPSG:4326）
clicked_sf_original <- sf::st_transform(clicked_sf_rot, crs = 4326)

# 3. 提取原始经纬度
original_coords <- sf::st_coordinates(clicked_sf_original)

# 输出结果，近似原始的-79.8524, 36.634
print(original_coords)

说明 #

sf::st_transform支持双向坐标转换，只要明确源CRS和目标CRS即可完成逆变换。关键是要确保点击的坐标被正确关联到旋转后的CRS，这样转换回WGS84时才能得到准确的原始经纬度。

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