reverse st_transform for real-world coordinates

Question

I have a process that rotates specific sf objects, and this works well for rendering, both in ggplot2 and plotly.

However, in shiny apps we need to take a user's click and do something useful with it, but the coordinates are in the transformed CRS, not the original, I'm uncertain how to reverse the transformation.

We use a proj4string and its gamma= argument to affect the rotation:

realpoint
#        lon    lat
# 1 -79.8524 36.634
realpoint_sf
# Simple feature collection with 1 feature and 0 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -79.8524 ymin: 36.634 xmax: -79.8524 ymax: 36.634
# Geodetic CRS:  WGS 84
#                  geometry
# 1 POINT (-79.8524 36.634)
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)
realpoint_sf_rot
# Simple feature collection with 1 feature and 0 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -33.41711 ymin: 57.92568 xmax: -33.41711 ymax: 57.92568
# Projected CRS: +proj=omerc +lat_0=36.634073240 +lonc=-79.851657820 +datum=WGS84 +alpha=0 +gamma=67.000
#                     geometry
# 1 POINT (-33.41711 57.92568)

### in a shiny app showing a ggplot including the above point,
### a user "click" returns data effectively like this:
userclick <- data.frame(sf::st_coordinates(realpoint_sf_rot))

How can I reverse the process to attain the original 79.85, 36.63 coordinates? I don't need absolute precision, "close" will be good enough, though I would prefer as accurate as theoretically reasonable.

---

Raw data:

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"))

Answer 1

You can create an sf object from the user input in your rotated CRS and transform it back to the original CRS like this:

library(sf)
userclick_sf_rot <- st_as_sf(userclick, coords=c('X', 'Y'), crs=crs_rot)
userclick_sf <- st_transform(userclick_sf_rot, crs=st_crs(realpoint_sf))

st_coordinates(userclick_sf)
#          X      Y
# 1 -79.8524 36.634
st_coordinates(realpoint_sf)
#          X      Y
# 1 -79.8524 36.634

Answer 2

The following might also be used to apply the inverse operation with respect to what is specified in a proj4string:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
original_point <- st_point(c(0, 0)) |> st_sfc(crs = 4326)
transformed_point <- st_transform(
  original_point, 
  crs = "+proj=omerc +lat_0=36.634073240 +lonc=-79.851657820 +datum=WGS84 +alpha=0 +gamma=67.000"
)

# back to the original point
st_transform(
  x = transformed_point, 
  pipeline = "+proj=pipeline +step +proj=omerc +lat_0=36.634073240 +lonc=-79.851657820 +datum=WGS84 +alpha=0 +gamma=67.000", 
  reverse = TRUE
)
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 4.440892e-16 ymin: 0 xmax: 4.440892e-16 ymax: 0
#> CRS:           NA
#> POINT (4.440892e-16 0)

^{Created on 2024-01-17 with reprex v2.0.2}

Unfortunately, we cannot recover the CRS. More details about pipeline operator can be found here.