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 NA
s 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"))
sp
orsf
. ANd mysp::CRS
doesn't have anSRS_string
parameter. What version of everything (packages, gdal, proj, as seen bylibrary(sf)
output) are you on? Edit these into your question?SRS_string
is new with recent version ofsp
and needs proj6, gdal3 and is part of the migration to WTK representation of CRS in R (see also many doc written by Roger Bivand). It should be the preferred way to specify CRS. I tried with Z coordinates as well but it fails too.