Transformación no finita detectada - lidR spTransform

Question

Me gustaría volver a proyectar el siguiente archivo laz utilizando lidR y spTransform de sp . El código EPSG no está definido, pero me parece que la cadena CRS está completa:

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"

Intenté proyectar al datum métrico global EPSG 3857 y obtuve el siguiente 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

He probado otra proyección (local), EPSG 6417 y obtuve el mismo 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

Bueno, muchas preguntas aquí

epsg es 0

Su archivo es un LAS 1.4. No almacena la proyección con un código epsg, sino con una cadena WKT. epsg() debería dar error (o al menos avisar) y decirte que tu archivo las no almacena código epsg sino cadena WKT. Pero en lugar de eso devuelve 0. Diseño mejorable...

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)
#> [...]

transformación no finita detectada

Este es un sp tema relacionado. No soy experto en proyección y transformación así que no puedo ayudarte. Usted debe hacer otra pregunta específica con la etiqueta enfocada. Mi conjetura es que usted está tratando de convertir cosas incompatibles. Aquí un ejemplo mínimo, independiente y reproducible con 3 puntos solamente:

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

Creado el 2020-10-01 por el paquete reprex (v0.3.0)

Uso de la pila C

Es probable que esto se deba a que tu nube de puntos es grande. No lo he reproducido porque sólo he cargado una fracción de los datos, pero estoy bastante seguro de que está relacionado con la cantidad de datos. Como se menciona en el post que enlazaste, spTransform es sólo una conveniencia para pequeños conjuntos de datos, pero no se recomienda su uso para realizar grandes transformaciones porque es inherentemente subóptima en términos de uso de memoria, ya que crea una copia de las coordenadas en sp formato. Utilización de las2las de LAStools es más aconsejable.

De todos modos, puedes intentarlo una vez que hayas resuelto tu transformación no finita. Tal vez esté relacionada y funcione.