How to add points to LAS copy using liblas?

Question

I am trying to make a copy of a LAS file while inserting new points using python liblas module. The code below creates the copy and adds the points, but for some reason the added points are not in the correct location. Also, it seems that x, y, z values are being stored int instead of float values. Can anyone explain why I am having this problem? Thanks in advance for your help.

import liblas
from liblas import file as lasfile

las_file = r"path here"
out_las = r"path here"

# Read las file
f = lasfile.File(las_file, mode='r')
h = f.header

# Open file to write to, using original header
fout = lasfile.File(out_las, mode='w', header=h)

counter = 0
for p in f:
    if counter != 0:
        # Code to Interpolate
        # x, y, z for new point
        # from existing points
        vD = (new_x, new_y, new_z)
        # vD is tuple of floats
        new_pt = liblas.point.Point()          # Create Point object
        new_pt.x = vD[0]                       # Set x property
        new_pt.y = vD[1]                       # Set y property
        new_pt.z = vD[2]                       # Set z property
        fout.write(new_pt)                     # Write new_pt to fout
    # Write p to fout, retaining original properties
    fout.write(p)
    counter += 1
fout.close()

Answer 1

I found my answer in the liblas.header documentation. It states the following:

"The scale factors in [x, y, z] for the point data. libLAS uses the scale factors plus the :obj:liblas.header.Header.offset values to store the point coordinates as integers in the file.

The scale factor fields contain a double floating point value that is used to scale the corresponding X, Y, and Z long values within the point records. The corresponding X, Y, and Z scale factor must be multiplied by the X, Y, or Z point record value to get the actual X, Y, or Z coordinate. For example, if the X, Y, and Z coordinates are intended to have two decimal point values, then each scale factor will contain the number 0.01

Coordinates are calculated using the following formula(s):

• x = (x_int * x_scale) + x_offset
• y = (y_int * y_scale) + y_offset
• z = (z_int * z_scale) + z_offset

.. note:: libLAS uses this header value when reading/writing raw point data to the file. If you change it in the middle of writing data, expect the unexpected."

That being said, I modified my original code to look like this:

# Get header
h = f.header
# Get scale factors
x_scale, y_scale, z_scale = h.scale[0], h.scale[1], h.scale[2]
# Get offset values
x_offset, y_offset, z_offset = h.offset[0], h.offset[1], h.offset[2]

for p in f:
    if counter != 0:
        # Code to Interpolate
        # x, y, z for new point
        # from existing points
        vD = (new_x, new_y, new_z)
        # vD is tuple of floats
        x, y, z = vD[0], vD[1], vD[2]
        new_pt = liblas.point.Point()           # Create Point object
        x_int = (x - x_offset) / x_scale        # De-scale x value
        y_int = (y - y_offset) / y_scale        # De-scale y value
        z_int = (z - z_offset) / z_scale        # De-scale z value                
        new_pt.x = x_int                        # Set x property
        new_pt.y = y_int                        # Set y property
        new_pt.z = z_int                        # Set z property
        fout.write(new_pt)                      # Write new_pt to fout
    # Write p to fout, retaining original properties
    fout.write(p)
    counter += 1
fout.close()