Rescaling and reoffsetting a point cloud with lidR?

Question

I'm trying to analyse a a 3 x 3 metres scan derived from terrestrial photogrammetry (i.e. local coordinates, with 0,0 in the centre of the point cloud) to obtain a DSM.

Importing the file in lidR, I get the following warnings:

Warning messages:
1: Invalid header: X scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 2.23339664936066e-09 
2: Invalid header: Y scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 2.33016455173492e-09 
3: Invalid header: Z scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 4.71837997436523e-10 

I specified the scale to 0.001 using lasrescale (suggested here too).

Then, I tried using lasreoffset to shift the point cloud to the correct geographical coordinates (here a similar implementation using LAStools), but I get this error:

las = lasreoffset(las, xoffset = 1306977.096, yoffset = 5661354.023, zoffset=265.8018)
#> Error: Incorrect xoffset: integer overflow.

When manually adding the offset to each point (e.g. las@data$X + 1306977.096, etc), the procedure works.

Finally, once I try to rasterize the point cloud using:

r = grid_canopy(las, res=0.01, dsmtin())

I get another error:

Error: Internal error in C_interpolate_delaunay: xy coordinates were not converted to integer. Scale factors are likely to be invalid.

It looks like I'm using the wrong tools for the purpose, or totally missing the point. Any suggestions?

Furthermore, can anybody help me understanding what does "integer overflow" refer to in this case?

Answer 1

This question is much more complex that it may appears.

The LAS format

First let's talk about the LAS format defined by the ASPRS. 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.

The LAS class in lidR

Second let's talk about the LAS class in lidR. It holds a point cloud in the LAS format and assumes that the object perfectly respects the LAS standard. For example, if you have a coordinate of 123456.789 with a scale factor of 0.01 and an offset of 0 this is not valid because (123456.789 - 0)/0.01 = 12345678.9 and this is not an integer, and thus this coordinate cannot be stored in a LAS file without information loss.

Rescaling and reoffseting

Now let's talk about las_rescale and las_reoffset. These functions intend to recompute valid coordinates according to the specification. Let's assume that your original las file contains 123456 with a scale of 0.0123 (which is not valid; it is a weird accuracy) and an offset of 0 (to simplify). At read time, the coordinate becomes 1234560.0123 + 0 = 1518.509. This is what is loaded in R. When rescaling, you want to transform the coordinates in such a way that the new values are valid with respect to the LAS specification. If you want a scale of 0.01, for example. We get back to the closest integer (1518.509 + 0)/0.01 = 151850.9 ~= 151851 and recompute the actual coordinate with this new scale factor 1518510.01 + 0 = 1518.51. 1518.51 is the closest value to 1518.509 with a scale of 0.01. And it is a valid value that can be converted to integer and stored in a LAS file. Reoffsetting is pretty much the same, but with the offset.

So you understand that rescaling and reoffseting are not about changing the coordinates but are about how the coordinates are stored in a LAS file.

Moving the coordinates

Now if you want to move manually the coordinates (las@data$X + 1306977.096) this is not that easy. Let's say you have X = 0 + 1306977.096 and a scale factor of 0.01 and an offset of 0. Now X is 1306977.096 which is not compatible with 0.01. You need to take care of the validity of the coordinates. That is easy to understand. More complex is the question of integer overflow.

Integer overflow is when an integer is greater than what can be stored in an integer. In a LAS file, coordinates are stored in 32 bits integer. Let's recompute 0 with a scale of 0.001 and an offset of 5661354.023. It becomes (0-5661354.023)/0.001 = -5661354023 which is greater than -2^31. This number is not an integer storable on 32 bits. You cannot store 0 with a scale factor of 0.001 and an offset of 5661354.023.

Specific issue with tin()

Most of the problem described here are only issues at write time. I mean, if you have 123456.789 with a scale factor of 0.1 this won't be written in a LAS file properly. But in R, all the functions work with the value 123456.789 and do not care of the scale factor and the offset.

All functions, but tin(). In lidR the triangulation is integer based. I won't explain why here (in short, for speed). The key is that your point cloud should always be convertible to integers since it is what is actually stored in the files. So lidR can take any coordinate, make an integer of it with the scale factor and the offset and compute the triangulation with integers.

In your case, you manually changed the coordinates; the LAS object is no longer valid and the conversion to integer failed because the header did not match with the coordinates, and the function tin() triggered an error.

In fact, it should have fell back to another decimal based (but slower) triangulation method. For an unknown reason, it failed. This might be reported as a bug.

Here we are! I told you that the question was complex.