[Solved] Custom LiDAR importer with python

[sant0s81]

Morning everyone,

a week ago I was starting to build a custom LiDAR Importer with the help of the guy who wrote a module to work with LiDAR data.
Actually it works, but there are some problems in speed, constant crashing etc.
Also the installation was quite a pain in the a$$ since I had to install lots of extra modules.

So I was searching the data base of python to find a more simple module and I think I found a really simple and quick one: https://pypi.org/project/laspy/

I also installed the daily buld with python 3 support since lots of the LiDAR modules are written in python 3 and I like the freedom of using Pip to install the modules.
I tested the example script on the page to load and write *.*las and its incredibly fast.
Even the different classifications of the LiDAR files can be filtered like seen in the screenshot.
The grey points are the original LiDAR, the roofs are filtered classifications with the Laspy module.

But right now I have to first write out the points and import again.

Is there a way to filter the points without saving but instantly show the result in Houdini?

I checked the manual for Laspy, but I could not find any usefull information.

Actually, instead of writing out the point positions, I have to fetch the xyz and use that to add points in Houdini I guess?

Any ideas, how that could be done?

Thanks alot
sant0s

 

edit: https://pythonhosted.org/laspy/util.html#laspy.util.Spec

There are some infos about the Util - also something with xml(). Could that help to fetch the data and use them in Houdini?

lidar_importer.hiplc

Edited July 11, 2020 by sant0s81

[sant0s81]

I just realised what Numpy actually is

Maybe with Numpy its possible, to get the xyz positions of the points and use them in Houdini?

Listening to a tutorial right now and damn, thats some crazy math for Einstein's.

But since Numpy also works with 3D arrays, it could be the way to go?

Edit:

I print(outFile.points) and I got a very promising output:

[((46300342, 6900720, 72207, 345, 73, 5, -4, 0, 587, 96018581.83408736),)outFile.close()
 ((46300252, 6900579, 72213, 291, 73, 5, -4, 0, 587, 96018581.79534031),)outFile.close()
 ((46308576, 6900400, 72220, 422, 73, 5, -7, 0, 587, 96018583.62939823),)outFile.close()
 ...outFile.close()
 ((46399155, 6999461, 82537, 178, 89, 5, 28, 0, 595, 96029297.37466702),)outFile.close()
 ((46399166, 6999254, 82147,  24, 90, 5, 28, 0, 595, 96029297.37466702),)outFile.close()
 ((46399155, 6999438, 82589,  95, 89, 5, 28, 0, 595, 96029297.37466954),)]outFile.close()
>>> outFile.close()

So something interesting is going on.
I wonder, why I get only 6 lines of the the point positions (if that are point positions).

Edited July 11, 2020 by sant0s81

[sant0s81]

Getting closer

from laspy.file import File
import numpy as np

inFile = File('Z:\lidar\GK_463_69.las', mode='r')

I = inFile.Classification == 2

outFile = File('Z:\lidar\output13.las', mode='w', header=inFile.header)
outFile.points = inFile.points[I]
X = inFile.X
Y = inFile.Y
Z = inFile.Z
print(X)
print(Y)
print(Z)


outFile.close()

 

With that I get my XYZ coordinates. Output looks like that:

[46300840 46300805 46300768 ... 46399155 46399166 46399155]outFile.close()
[6900092 6900066 6900038 ... 6999461 6999254 6999438]outFile.close()
[72201 72205 72217 ... 82537 82147 82589]outFile.close()

 

Are these '...' actually representing the rest of the coordinates but to not have millions of numbers, its shortening the print output?

But with that I now could feed something in Houdini.

Probably very similar to the CSV module reading GPS date etc.

If someone has ideas, my head is burning already

[Atom]

You could check out my Moon DEM Python reader. I believe those original NASA files were in a similar format.

 

[sant0s81]

Update:

I am trying to modify the CSV reader to read the coordinates, but something there is not working - brain overheating.

from laspy.file import File
import numpy as np

inFile = File('Z:\lidar\GK_463_69.las', mode='r')

I = inFile.Classification == 2

outFile = File('Z:\lidar\output13.las', mode='w', header=inFile.header)
outFile.points = inFile.points[I]

x = inFile.X
y = inFile.Y
z = inFile.Z

#print(X)
#print(Y)
#print(Z)

with outFile as f:
    reader = pt_reader(f)
    
    for i in pt_reader:
        if i == 0:
            continue
        pt = geo.createPoint()
        x = float(row[x])
        y = float(row[y])        
        z = float(row[z])

outFile.close()

 

Edited July 11, 2020 by sant0s81

[sant0s81]

@Atom

Thanks for the link - just trying to understand your code there

But do you have an idea, whats wrong with my last with-for script?
With CSV it works similar - but I guess I am mixing something since the 'with outFile' gives me already the error.

[Atom]

It looks like your are trying to read from an out file that was opened with mode 'w'. Wouldn't you need to read from the in file?

reader = pt_reader(inFile)

I think I would use standard Python file operations. Do your LIDAR files meet the requirements of the laspy API?

Quote

limited to reading LAS version 1.0-1.3 files.

 

Edited July 11, 2020 by Atom

[sant0s81]

39 minutes ago, Atom said:

It looks like your are trying to read from an out file that was opened with mode 'w'. Wouldn't you need to read from the in file?

reader = pt_reader(inFile)

 

Actually I dont need the outFile since I dont want to write out a new *.*las - and I am not sure, but do I actually need the 'with'?

Wouldnt something like that already come closer:

import hou
import re

from laspy.file import File
import numpy as np

inFile = File('Z:\lidar\GK_463_69.las', mode='r')

I = inFile.Classification == 2

#outFile = File('Z:\lidar\output1.las', mode='w', header=inFile.header)
inFile.points = inFile.points[I]

x = inFile.X
y = inFile.Y
z = inFile.Z

#print(X)
#print(Y)
#print(Z)

for i in inFile:        
    if i == 0:
        continue
    pt = geo.createPoint()
    x = float(row[x])
    y = float(row[y])        
    z = float(row[z])
        
    pt.setPosition(hou.Vector3(x, y, z))

 

Edited July 11, 2020 by sant0s81

[sant0s81]

Btw, here is the scene file and the wetransfer to the *.*las file (Its from slovenia geo department, so available for public, but the website seems to be down at the moment and I dont find any other example file)

Just in case someone wanna test around

https://wetransfer.com/downloads/4f62b5945e9108a9dea3fdd9b965a24820200711140423/c0518a1995a0223d571f55ec142e893c20200711140439/1921cf

and the link to the LiDAR module:
https://pypi.org/project/laspy/

lidar_importer_v1.hiplc

Edited July 11, 2020 by sant0s81

[pezetko]

Hi, pretty neat library! Thank you for the tip.

There is no need for csv, you can do a lot with laspy and numpy itself.

Attached example scene to load data from las file. Seems that Lidar Import SOP ignores scale and offset.

 

To make it work (18.0.499, Python 2.7 branch) I cloned the https://github.com/laspy/laspy repository. Then copied content of laspy folder to $HOME/Houdini18.0/python2.7libs/laspy so I have $HOME/houdini18.0/python2.7libs/laspy/__init__.py (and rest of the library) and it's possible to load it into Houdini with import laspy in Python shell. (Numpy is already included with Houdini)

I used example file from repository: https://github.com/laspy/laspy/blob/master/laspytest/data/simple.las

import logging
from laspy.file import File
import numpy as np

node = hou.pwd()
geo = node.geometry()

file_path = geo.attribValue("file_path")
inFile = File(file_path, mode='r')

try:
	# --- load point position
	coords = np.vstack((inFile.X, inFile.Y, inFile.Z)).transpose()
	scale = np.array(inFile.header.scale)
	offset = np.array(inFile.header.offset) # there is no offset in simple.las example from laspy library
	# offset = np.array([1000, 20000, 100000])  # just for testing that offset works

	# geo.setPointFloatAttribValues("P", np.concatenate(coords))    # same as Lidar Import SOP - seems that it ignores scale (and offset?)
	geo.setPointFloatAttribValues("P", np.concatenate(coords*scale+offset))

	# --- load color
	color = np.vstack((inFile.red, inFile.green, inFile.blue)).transpose()
	geo.addAttrib(hou.attribType.Point, "Cd", (1.0,1.0,1.0), False, False)  # add color atttribute
	geo.setPointFloatAttribValues("Cd", np.concatenate(color / 255.0)) # transform from 1-255 to 0.0-1.0 range)

except Exception:
	logging.exception("Processing lidar file failed")
finally:
	inFile.close()

 

pz_load_las_with_python.hipnc

Edited July 15, 2020 by pezetko
fixed inFile.x vs inFile.X

[sant0s81]

@pezetko

awesome - gonna try it in some minutes.

Did you try

 inFile.Classification == x

?

Thats the idea of the custom LiDAR importer.
Than you can for example build a Heightfield only from the ground points(in most cases its classification 2) and than filter out only the vegetation (normaly class 3,4 and 5), buildings (class 6), water (class 9) and build HF masks from that.

But thanks alot, cannot wait to test, that will save so much time instead of going over a GIS application

 

Edited July 11, 2020 by sant0s81

[sant0s81]

@pezetko

works perfect, so good!

But somehow I dont get the classification anymore.
I was trying something to change your script to that, but it doesnt work and leads to a crash:

file_path = hou.evalParm("lidar_file")
I = inFile.Classification == 2
inFile = File(file_path, mode='r')
inFile = inFile.points[I]

Do you have an idea why that doesnt work?
When exporting with outFile it works.

[sant0s81]

And here a link to the classifications:

https://geodetics.com/lidar-point-clouds/

[pezetko]

The classification works fine. But with

inFile = inFile.points[I]

you are overwriting inFile object with multi-dimensional arrays of all attributes so you no longer get them by .x/.y/.z or other named properties.

I uploaded a modified scene, where you can set classification and then it returns only a subset of points that match the correct classification.

inFile.points[I]

Returns subset of points in an array where I is True

inFile.Classification == 2

Returns array of True/False values to determine which points are classified with id 2.

 

Another approach would be adding Classification as an attribute to all points and then use vex expressions, groups, or other partitioning mechanisms to separate points.

pz_load_las_with_python_classification.hipnc

[sant0s81]

@pezetko

That is freaking awesome!!!!!

Thank you so much for your help.
I hope SideFX gonna read that post and implement a better LiDAR Importer in on of there future releases since you can do so much cool stuff with it now.

Muito obrigado!

Here a screenshot of the seperated classifications. Still smiling all over my face

Edited July 11, 2020 by sant0s81

[pezetko]

You are welcome.

Submitting RFEs doesn't hurt. C++ implementation is still much faster than laspy, but it's possible to do an easy and quick specific modification for las import until LIDAR Import SOP gets some improvements without going the C++ route.

I attached an example with adding classification as an attribute (added closing the las file to avoid memory issues).

I added an updated file: pz_load_las_with_python_classification_attribute_memory_fix.hipnc as it looks like context management protocol is not implemented in laspy (with File(): block will not close the file at exit) so I switched to try/except/finally instead. It will not error on the node, so watch python console for exception logging.

 

Edited July 15, 2020 by pezetko
Fixed memory issues with closing file in try block

[pezetko]

If you have any memory issue just check latest update version that using try/except block https://forums.odforce.net/applications/core/interface/file/attachment.php?id=56259

I'm not sure if laspy does implement context management protocol (so using "with File() as file_pointer:" statement may not close the file at exit) Just verified that and with File() as fp: statement works as it should.

Edited July 11, 2020 by pezetko
with context works just fine

[bunker]

just in case, did you try this?
https://www.sidefx.com/docs/houdini/nodes/sop/lidarimport.html
edit: worked for me with H18.0.327, even on my 11 years old macbook pro, it takes under a sec to load 

Edited July 11, 2020 by bunker

[sant0s81]

2 hours ago, bunker said:

just in case, did you try this?
https://www.sidefx.com/docs/houdini/nodes/sop/lidarimport.html
edit: worked for me with H18.0.327, even on my 11 years old macbook pro, it takes under a sec to load 

Hey @bunker

yes - the build in LiDAR Import is cool.
But it has no settings to filter the classifications and that was the goal, at least for me.

What do you mean loading under a second to load?
Also the wetransfer I uploaded some posts above?
Damn, you must have a NASA super computer

@pezetko

gonna go for a nap now, than I try.
Thx again for that lessons in python. Learned so much today!
Already used the script in an HDA to use in UE4 - gonna show some stuff tomorrow!

Have a good night guys,
sant0s
 

Edited July 12, 2020 by sant0s81

[sant0s81]

@pezetko

Just playing around with your most actuall version. Working like expected.

Maybe I am missing something, but is there a drop down menu for the classifications I dont see?

# Add code to modify contents of geo.
# Use drop down menu to select examples.

But anyway - that tool just blows it and makes the workflow between Houdini -> UE4 so smooth.

Cheers,
sant0s

 

Edit:

18 hours ago, pezetko said:

To make it work (18.0.499, Python 2.7 branch) I cloned the https://github.com/laspy/laspy repository. Then copied content of laspy folder to $HOME/Houdini18.0/python2.7libs/laspy so I have $HOME/houdini18.0/python2.7libs/laspy/__init__.py (and rest of the library) and it's possible to load it into Houdini with import laspy in Python shell. (Numpy is already included with Houdini)

I also just installed Laspy to get it work with Python 2.7 - but instead of copying the content of the folder into python2.7libs, I had to copy it into $HOME\Houdini 18.0.499\python27\lib - dont know why, but just in case someone gets error that Laspy module is not available.

Edited July 12, 2020 by sant0s81