Jump to content

[Solved] Custom LiDAR importer with python


sant0s81

Recommended Posts

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_import.thumb.jpg.ef4848746a1adf9d08271e6515596915.jpg

lidar_importer.hiplc

Edited by sant0s81
Link to comment
Share on other sites

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 by sant0s81
Link to comment
Share on other sites

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 :D

Link to comment
Share on other sites

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 by sant0s81
Link to comment
Share on other sites

@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.

Link to comment
Share on other sites

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 by Atom
Link to comment
Share on other sites

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 by sant0s81
Link to comment
Share on other sites

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 by sant0s81
Link to comment
Share on other sites

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 by pezetko
fixed inFile.x vs inFile.X
  • Like 2
Link to comment
Share on other sites

@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 :D

 

Edited by sant0s81
Link to comment
Share on other sites

@pezetko

works perfect, so good! :D

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.

Link to comment
Share on other sites

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

  • Like 1
  • Thanks 1
Link to comment
Share on other sites

@pezetko

That is freaking awesome!!!!! :D

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 :D

lidar_import2.thumb.jpg.813c8cf1b55056c48114354f6925670e.jpg

Edited by sant0s81
Link to comment
Share on other sites

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.

laspy_lidar.png.43ae43f4233ee8af21cdfbf2f74cf7c6.png

 

Edited by pezetko
Fixed memory issues with closing file in try block
  • Like 2
Link to comment
Share on other sites

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 by pezetko
with context works just fine
Link to comment
Share on other sites

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 :D

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 by sant0s81
Link to comment
Share on other sites

@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 by sant0s81
Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
×
  • Create New...