4

I have an input LAS file that has lots of points, that also includes RGB and other info besides the XYZ. My goal is to select only the XYZ of the ground points, preserving the geographic location, but discard non-ground points and discard RGB and other info from the ground points I am keeping. Here is what I have tried:

orig_las = laspy.read('original.las')
ground_points = orig_las.points[orig_las.points.classification == 2]
new_las = laspy.LasData(orig_las.header)
new_las.X = ground_points.X
new_las.Y = ground_points.Y
new_las.Z = ground_points.Z
new_las.write('new.las'))

However it fails trying to set new_las.X with ValueError: could not broadcast input array from shape (42,) into shape (17934,) I tried calling new_las.update_header() before writing XYZ but it didn't seem to help.

I also tried creating the LasData like this:

new_las = laspy.create(point_format=orig_las.header.point_format, file_version=orig_las.header.version)

This succeeds in creating a file with my 42 ground points, however it has lost the information about where these points are on the earth. (i.e. when I drag them both into QGIS, they do not overlap).

I notice the original has a vlr that has the CRS in it, so by doing this I am able to copy over the vlrs from the original:

new_las.header.vlrs = orig_las.header.vlrs

But they still do not appear in the same location in QGIS. I notice that the header.offsets and header.scales are still different, as well as the points.offsets and points.scales... I have not tried copying these over but I can't imagine that copying over all these values from one LasData to the other is the best solution. I figure I must be missing an easier way.

I did see the examples to do chunked writing to copy points using write_points, but the problem is I don't want to copy all the data from the source points; I am looking to reduce the dimensions to just the XYZ data to save space (removing RGB and other info).

1 Answer 1

2

To filter your file and write it without losing metadata doing the following should be enough

ground_las = orig_las[orig_las.classification == 2]
ground_las.write("ground.las")

In the code you posted, I think you are missing resetting the point count to 0.

from copy import copy


ground_points = orig_las.points[orig_las.points.classification == 2]
header = copy(orig_las.header)
header.point_count = 0
new_las = laspy.LasData(header)
new_las.X = ground_points.X
new_las.Y = ground_points.Y
new_las.Z = ground_points.Z
new_las.write('new.las')

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.