best practice for constructing compound CRSes from EPSGs?

[twest820]

st_crs(epsg) calls OGRSpatialReference::importFromEPSG() by way of OGRSpatialReference::SetFromUserInput() and therefore can handle either a projected EPSG or a vertical EPSG, but not both at the same time. There doesn't seem to be a way to call OGRSpatialReference::SetCompoundCS() from R.

Cleanest workaround I've found for this using R code is to create both components and glue their wkts together. Function below seems to work fine but comparative testing against layers generated by code compiled against GDAL 3.10.2 APIs shows differences in whether or not the geoid gets included in the vertical CRS's name. While this appears to be of little functional importance (and is maybe just inconsistent behavior between SetCompoundCS() and importFromEPSG()), I'm wondering if there's a better approach short of dropping into C++.

make_compound_crs = function(projectedEpsg, verticalEpsg) # EPSG numbers or "EPSG:<number>" strings
{
  projectedCrs = st_crs(projectedEpsg)
  projectedCrsName = stringr::str_extract(projectedCrs$wkt, "PROJCRS\\[\"([A-Za-z0-9_\\[\\]\\(\\){}<=>\\.,:;\\+\\- #%&'*^/\\@|°]+)", group = 1) # <quoted Latin text> except <doublequote symbol>, https://docs.ogc.org/is/18-010r11/18-010r11.pdf
  verticalCrs = st_crs(verticalEpsg)
  verticalCrsName = stringr::str_extract(verticalCrs$wkt, "VERTCRS\\[\"([A-Za-z0-9_\\[\\]\\(\\){}<=>\\.,:;\\+\\- #%&'*^/\\@|°]+)", group = 1)
  
  compoundWkt = paste0("COMPOUNDCRS[\"", projectedCrsName, " + ", verticalCrsName, "\",\n    ", projectedCrs$wkt, ",\n    ", verticalCrs$wkt, "]")
  return(st_crs(compoundWkt))
}

[edzer]

I see no problems with calling the C++ code; while at it, we may want to allow for setting the epoch.

[rsbivand]

A side comment on the proj list (I think) pointed to the imminent demise of NAD83 and NAVD88 https://www.ngs.noaa.gov/datums/newdatums/index.shtml as an expected shock to CRS and transformations. This also involves epochs.