1

I am trying to reproject a las file in lidR (doawnload file here) from epsg:2232 to epsg:32613, and get the following error message after about 3 seconds of processing:

las <- readLAS("LD31201353.las")
projection(las) <- sp::CRS("EPSG:2232")
spTransform(las, CRSobj = "+init=epsg:32613")
#> Error in fast_quantization(X, scalex, offsetx) : 
#>  Non quantizable value outside the range of representable values of type 'int'

Has anyone encountered this, and might you have suggestions for a fix (or what I'm doing wrong)?

The las file natively does not have its 'WKT OGC COORDINATE SYSTEM' field or projargs slot populated, but I get the error with and without assigning them.

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 
#> Project ID - GUID:        00000001-0012-07e3-0000-000000000000 
#> Version:                  1.4
#> System identifier:        Merrick 
#> Generating software:      MARS 
#> File creation d/y:        188/2018
#> header size:              375 
#> Offset to point data:     377 
#> Num. var. length record:  0 
#> Point data format:        6 
#> Point data record length: 30 
#> Num. of point records:    1965664 
#> Num. of points by return: 1634171 285040 42901 3469 82 1 0 0 0 0 0 0 0 0 0 
#> Scale factor X Y Z:       0.001 0.001 0.001 
#> Offset X Y Z:             3122231 1354500 11639 
#> min X Y Z:                3121464 1353000 11008.84 
#> max X Y Z:                3123000 1356000 11782.48 
#> Variable length records: 
#>    Variable length record 1 of 1 
#>        Description: WKT Information 
#>        WKT OGC COORDINATE SYSTEM: PROJCRS["NAD83 / Colorado Central (ftUS)",
#>     BASEGEOGCRS["NAD83",
#>    [...] (truncated)
6
  • Anyway you can share your data or a subset of it that exhibits this problem? Or can you show your complete workflow as well as summary output and plots of your las object? What package are you using for LAS objects?
    – Spacedman
    Jan 28, 2021 at 21:39
  • @Spacedman - sorry, new to this, I don't see an upload file option. 11MB is the size of a 5x5m sample, anyway.
    – G. Moore
    Jan 28, 2021 at 22:08
  • @Kahdir Sahbaz - thank you for the formatting!
    – G. Moore
    Jan 28, 2021 at 22:08
  • @Spacedman: Win 10 Home, rgdal_1.5-21, lidR_3.1.1, raster_3.4-5, sp_1.4-5
    – G. Moore
    Jan 28, 2021 at 22:19
  • You'll need a third-party to share data files - dropbox, onedrive, box, or some other file storage. I've added "lidr" to the tags since some people have searches triggered on tags.
    – Spacedman
    Jan 28, 2021 at 23:22

1 Answer 1

0

So, first, we need a reproducible example

library(lidR)
data = data.frame(X = c(3176138, 3173598),Y = c(1434968, 1432401), Z = 0)
header = rlas::header_create(data)
header[["X offset"]] = 3172500
header[["Y offset"]] = 1432500
header[["X scale factor"]] = 0.001
header[["Y scale factor"]] = 0.001
las = LAS(data, header)
projection(las) = 2232

spTransform(las, CRSobj = "+init=epsg:32613")
#> Error in fast_quantization(X, scalex, offsetx) : 
#>  Non quantizable value outside the range of representable values of type 'int'

The problem comes from the fact that in a las file, coordinates are stored as integers. Because in lidR, LAS object are strictly (as far as possible) compliant to ASPRS specification of the LAS format it sometimes add a layer of complexity and triggers some error difficult to understand.

I already provided a comprehensive explanation of the LAS format and the importance of the scale and offset factors here

In a las file coordinates are stored as integers with a scale factor and an offset to compute back decimal positioning. For example, a X coordinate might be 123456 with a scale factor 0.01 and offset of 100000 for an actual coordinate of 123456*0.01+100000 = 101234.56. The scale factor is the accuracy of your point (0.01 is a centimeter accuracy if the units are meters) and the offset enable memory and accuracy storage optimizations. This is the meaning of scale and offset.

In your case the scale factors are 0.001 and offsets are 3172500 and 1432500. When re-projecting from 2232 to 326013 the new coordinates range in [509611.1, 510389.3, 4318892, 4319670]

spTransform(as.spatial(las), CRSobj = "+init=epsg:32613")
#> class       : SpatialPointsDataFrame 
#> features    : 2 
#> extent      : 509611.1, 510389.3, 4318892, 4319670  (xmin, xmax, ymin, ymax)
#> crs         : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
#> variables   : 1
#> names       : Z 
#> min values  : 0 

The problem is that new scale and offset factors are unknown, thus lidR re-use the old ones to quantize the value. Let do the math for Y

(4319670- 1432500)/0.001 = 2886392000

which is more than 2^31 and thus cannot be stored as integer with a scale factor of 0.001 and an offset of 1432500.

Now we understand the problem there is no solution because there is a bug in the function. The offsets are actually expected to be recomputed in spTransform but are not. This will be fixed in lidR v3.1.2 (already fixed on github).

In addition, for consistency across multiple tiles is is possible to set the offsets manually (in v >= 3.1.2)

spTransform(las, CRSobj = "+init=epsg:32613", xoffset = 509600, yoffset = 4319600)
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.