Non finite transformation detected - lidR spTransform

Question

I would like to re project the following laz file using lidR and spTransform from sp. The EPSG code is not defined but it seems to me that the CRS string is complete:

tmp_laz <- readLAS("./USGS_LPC_CA_FEMA_R9_Russian_2017_10SDJ993053_LAS_2018.laz")
epsg(tmp_laz)
#>[1] 0
tmp_laz@proj4string@projargs
#> [1] "+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"

I tried to project to the global metric datum EPSG 3857 and got the following error non finite transformation detected:

tmp_laz_2 <- spTransform(tmp_laz, sp::CRS("+init=epsg:3857") )
#> non finite transformation detected:
#>   X       Y        
#> [1,] 6275252 2094773 Inf Inf
#> [2,] 6275249 2094771 Inf Inf
#> [3,] 6275246 2094769 Inf Inf
#> [4,] 6275243 2094768 Inf Inf
#> [5,] 6275251 2094768 Inf Inf
#> [6,] 6275251 2094775 Inf Inf
#> [7,] 6275248 2094774 Inf Inf
#> [8,] 6275245 2094772 Inf Inf
#> [9,] 6275242 2094770 Inf Inf
#> [10,] 6275239 2094768 Inf Inf
#> ... goes on till reaching print limit
#> Erreur : C stack usage  100062272 is too close to the limit
#> De plus : 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
#> 3: In sp::proj4string(x) : CRS object has comment, which is lost in output
#> 4: In sp::spTransform(spts, CRSobj) :
#>   12346945 projected point(s) not finite

I tried another (local) projection, EPSG 6417 and got the same error:

tmp_laz_6417 <- spTransform(tmp_laz, sp::CRS("+init=epsg:6417") )

#> non finite transformation detected:
#>   X       Y        
#> [1,] 6275252 2094773 Inf Inf
#> [2,] 6275249 2094771 Inf Inf
#> [3,] 6275246 2094769 Inf Inf
#> [4,] 6275243 2094768 Inf Inf
#> [5,] 6275251 2094768 Inf Inf
#> [6,] 6275251 2094775 Inf Inf
#> [7,] 6275248 2094774 Inf Inf
#> [8,] 6275245 2094772 Inf Inf
#> [9,] 6275242 2094770 Inf Inf
#> [10,] 6275239 2094768 Inf Inf
#> ... goes on till reaching print limit

#> Erreur : C stack usage  100062272 is too close to the limit
#> De plus : Warning messages:
#>   1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
#>   Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition
#> 2: In sp::proj4string(x) : CRS object has comment, which is lost in output
#> 3: In sp::spTransform(spts, CRSobj) :
#>   12346945 projected point(s) not finite

Answer 1

Well, many questions here

epsg is 0

Your file is a LAS 1.4. It does not store the projection with an epsg code but with a WKT string. epsg() should error (or at least warn) and tell you that your las file does not store epsg code but WKT string. But it returns 0 instead. Improvable design...

tmp_las@header
#> File signature:           LASF 
#> File source ID:           0 
#> Global encoding:
#>  - GPS Time Type: Standard GPS Time 
#>  - Synthetic Return Numbers: no 
#>  - Well Know Text: CRS is WKT 
#>  - Aggregate Model: false 
#> [...]
#> Variable length records: 
#>    Variable length record 1 of 3 
#>        Description: OGC WKT Coordinate System 
#>        WKT OGC COORDINATE SYSTEM: COMPD_CS["NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geo [...] (truncated)
#> [...]

non finite transformation detected

This one is an sp related issue. I'm not an expert in projection and transformation so I can't help you. You should ask another specific question with focused tag. My guess is that you are trying to convert incompatible things. Here a minimal, standalone, reproducible example with 3 points only:

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

^{Created on 2020-10-01 by the reprex package (v0.3.0)}

C stack usage

This one is likely because your point cloud is big. I didn't reproduce because I only loaded a fraction of the data but I'm pretty sure it is related to the amount of data. As mentioned in the post you linked, spTransform is only a convenience for small datasets but it is not recommended to use to perform big transformation because it is inherently sub optimal in term of memory usage as it creates a copy of the coordinates in sp format. Using las2las from LAStools is more advised.

You can give it a try anyway once you solved your non finite transformation. Maybe it is related and it will work.