Reproject using WKT for just the horizontal CRS

Question

I'm trying to reproject LiDAR point clouds from California zone 2 (ftUS) using spTransform from sp R package to a metric CRS like EPSG 3857 (initial post here) and got non finite transformation detected: error. Here bellow a reproducible example.

library(sp)
X = c(6275235.88, 6275229.48, 6275163.35)
Y = c(2094771.54, 2094815.15, 2094777.31)
WKT = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",39.83333333333334],PARAMETER[\"standard_parallel_2\",38.33333333333334],PARAMETER[\"latitude_of_origin\",37.66666666666666],PARAMETER[\"central_meridian\",-122],PARAMETER[\"false_easting\",6561666.667],PARAMETER[\"false_northing\",1640416.667],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]],VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AUTHORITY[\"EPSG\",\"6360\"]]]"
XY = cbind(X,Y)
XY = sp::SpatialPoints(XY, sp::CRS(SRS_string = WKT))
sp::spTransform(XY, sp::CRS("+init=epsg:3857"))
#> Warning in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): 3 projected point(s)
#> not finite
#> non finite transformation detected:
#>            X       Y        
#> [1,] 6275236 2094772 Inf Inf
#> [2,] 6275229 2094815 Inf Inf
#> [3,] 6275163 2094777 Inf Inf
#> Error in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): failure in points 1:2:3

Basically, in this project I'm dealing with various CRS across the world and I don't need absolute Z accuracy but I really want to get rid of US feet survey CRS !! The EPSG 3857 is convenient because it is a worldwide metric datum but it is not the cause of the error. I received the same error using local datum like EPSG 6417 and 6418. As suggested by @mkennedy in this post, the Z component of the CRS transformation could be the cause of the error and I am wandering how to use spTransform just for the horizontal CRS.

In other words, is it possible to ask to spTransform to focus on the horizontal component of the targeted CRS? I know this is related to the proj library but I'm looking for a solution within R environment. I'm open to solutions based on other R package like sf

Edit: the problem is reproducible with a Z coordinate or with sf::st_transform(). But sf doesn't fail. It returns NAs instead

X = c(6275235.88, 6275229.48, 6275163.35)
Y = c(2094771.54, 2094815.15, 2094777.31)
Z = c(10,15,20)
WKT = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",39.83333333333334],PARAMETER[\"standard_parallel_2\",38.33333333333334],PARAMETER[\"latitude_of_origin\",37.66666666666666],PARAMETER[\"central_meridian\",-122],PARAMETER[\"false_easting\",6561666.667],PARAMETER[\"false_northing\",1640416.667],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]],VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AUTHORITY[\"EPSG\",\"6360\"]]]"
XYZ = cbind(X,Y,Z)
XYZ = sp::SpatialPoints(XYZ, sp::CRS(SRS_string = WKT))
sp::spTransform(XYZ, sp::CRS("+init=epsg:3857"))

XYZ = sf::st_as_sf(XYZ)
sf::st_transform(XYZ, sf::st_crs("EPSG:3857"))

Answer 1

You can turn this into a 2D transform by removing the VERT_CS component of the WKT.

WKT2 = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",
PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",
GEOGCS[\"NAD83(2011)\",
DATUM[\"NAD83_National_Spatial_Reference_System_2011\",
SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],
AUTHORITY[\"EPSG\",\"1116\"]],
PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],
UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],
AUTHORITY[\"EPSG\",\"6318\"]],
PROJECTION[\"Lambert_Conformal_Conic_2SP\"],
PARAMETER[\"standard_parallel_1\",39.83333333333334],
PARAMETER[\"standard_parallel_2\",38.33333333333334],
PARAMETER[\"latitude_of_origin\",37.66666666666666],
PARAMETER[\"central_meridian\",-122],
PARAMETER[\"false_easting\",6561666.667],
PARAMETER[\"false_northing\",1640416.667],
UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],
AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]]]"

and then I only get the warning messages and what looks like a sensible transformation of the X and Y coordinates:

> XYZ = sp::SpatialPoints(XYZ, sp::CRS(SRS_string = WKT2))
Warning message:
In showSRID(SRS_string, format = "PROJ", multiline = "NO") :
  Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition
> sp::spTransform(XYZ, sp::CRS("+init=epsg:3857"))
SpatialPoints:
             X       Y       Z
[1,] -13693038 4708772 2.0e+07
[2,] -13693041 4708789 1.5e+01
[3,] -13693066 4708774 2.0e+01
Coordinate Reference System (CRS) arguments: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m
+nadgrids=@null +wktext +no_defs 
Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum WGS_1984 in CRS definition

the heights being unchanged (so multiply by meter-per-USft to fix).

By sensible I mean if I transform them to lat-long I get something that seems to be around California somewhere...

> sp::spTransform(XYZ, sp::CRS("+init=epsg:4326"))
SpatialPoints:
             X        Y       Z
[1,] -123.0067 38.90989 2.0e+07
[2,] -123.0067 38.91001 1.5e+01
[3,] -123.0069 38.90990 2.0e+01

I don't know why a VERT_CS messes this up, possibly worth a bug report, but I'm not sure who to...

Versions:

> library(sf)
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

Also I don't know if the vertical coordinate and horizontal components are independent. In other words if transforming (X,Y,Z) to (X', Y', Z') can be written X'=f(X,Y), Y'=g(X,Y), Z'=h(Z). If not then dropping the vertical coordinate system might result in errors because of a difference in the definition of "vertical" in different systems.

To investigate if this is a bug, I used R to convert your 3d CRS to a Proj string, and feed into the proj command line tools to see what it did:

# cs2cs +proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333 +lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016 +ellps=GRS80 +units=us-ft +vunits=us-ft +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs
6275235.88 2094771.54 100
-13693037.82    4708772.00 30.48

And this looks correct, even down to the change of height from 100 us-ft to 30.48 metres. Feeding this proj string into R still returns empty points on transform:

> WKTc = st_crs(WKT)
> projstring = as(WKTc,"CRS")
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
> projstring
CRS arguments:
 +proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333
+lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016
+ellps=GRS80 +units=us-ft +vunits=us-ft +no_defs 
> pts = st_as_sf(data.frame(x=6275235.88,y= 2094771.54,z=100), coords=1:3, crs=projstring)
> st_transform(pts, "+init=epsg:3857")
Simple feature collection with 1 feature and 0 fields (with 1 geometry empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
projected CRS:  WGS 84 / Pseudo-Mercator
     geometry
1 POINT EMPTY

which seems to be because R has converted the proj string back to a WKT. Given that (I think) proj internally only deals with proj strings, it might be a bug in the way R is passing the coordinates to Proj, because it seems capable of generating the correct proj string (or my cs2cs command line wouldn't work).