"non finite transformation detected" using sp R

Question

This is a part of a broader question already posted (Non finite transformation detected - lidR spTransform) but I am creating a new question for the sp related part of that.

I'm trying to reproject points using spTransform and get non finite transformation detected error. Below is a reproducible example with which I was helped by @JRR:

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

I also got the same error when projecting to other CRS like EPSG 6417. I know this must be because I'm asking something mathematically impossible but I don't see what as it looks like all the info are available to perform the reprojection.

Answer 1

The problem was indeed related to the vertical CRS. If you are looking for a quick (and a little dirty if you are interested about Z absolute accuracy after re-projection) you can drop the VERT_CS[] (thanks @Spacedman, see Reproject using WKT for just the horizontal CRS). In the code below, the solution integrates a regex to remove the problematic VERT_CS part (thanks to this post)

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\"]]]"

# remove VERT_CS part of WKT
WKT2 = gsub("VERT_CS\\[.*?\\]+,", "", WKT)
XY = cbind(X,Y)
XY = sp::SpatialPoints(XY, sp::CRS(SRS_string = WKT2))
sp::spTransform(XY, sp::CRS("+init=epsg:3857"))

Answer 2

Few months later, I would like to report that this is no longer an issue and the trick is not needed anymore... I guess this was a bugg