[GIS] Visualizing LAS with matplotlib

Question

I'm trying to create a 3D scatterplot using an las file that contains xyz and return values. I want to keep all returns but would like to subset the cloud. I've learned how to subset using laspy so no issue there. The trouble comes when I'm trying to load the las to matplotlib.

How do I do I load the las to matplotlib? Is there a better way to view it? I'm doing it this way so I have scales around the outside of the las and 3d view capability.

The code below, loads two points that are at the extents of the file, upper left and lower right. I want to load those and all the points in between.

if I try to use the line

x,y,z = points

I get the error

Traceback (most recent call last):
  File "<stdin>", line 1 in <module>
ValueError: too many values to unpack

Code:

from liblas import file
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
las_file = r"path/to/las"
f = file.File(las_file, mode='r')
num_points = int(f.__len__())
PointsXYZ = np.empty(shape=(num_points, 3))
counter = 0
for p in f:
    newrow = [p.x, p.y, p.z]
    PointsXYZ[counter] = newrow
    counter += 1    

ax = fig.add_subplot(111, projection= '3d')

ax.set_axis_bgcolor('black')
ax.xaxis.label.set_color('white')
ax.yaxis.label.set_color('white')
ax.zaxis.label.set_color('white')
ax.tick_params(axis='x',colors='white')
ax.tick_params(axis='y',colors='white')
ax.tick_params(axis='z',colors='white')

listmin = np.amin(PointsXYZ, axis=0)
listmax = np.amax(PointsXYZ, axis=0)
zmin = listmin[2]
xmin = listmin[1]
ymin = listmin[0]
zmax = listmax[2]
xmax = listmax[1]
ymax = listmax[0]

xs = (counter, xmin, xmax)
ys = (counter, ymin, ymax)
zs = (counter, zmin, zmax)

ax.scatter(xs, ys, zs, c='r', marker='o')
plt.show

Answer 1

This is a working solution with laspy:

import numpy as np
import laspy
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# reading las file and copy points
input_las = laspy.file.File("test.las", mode="r")
point_records = input_las.points.copy()

# getting scaling and offset parameters
las_scaleX = input_las.header.scale[0]
las_offsetX = input_las.header.offset[0]
las_scaleY = input_las.header.scale[1]
las_offsetY = input_las.header.offset[1]
las_scaleZ = input_las.header.scale[2]
las_offsetZ = input_las.header.offset[2]

# calculating coordinates
p_X = np.array((point_records['point']['X'] * las_scaleX) + las_offsetX)
p_Y = np.array((point_records['point']['Y'] * las_scaleY) + las_offsetY)
p_Z = np.array((point_records['point']['Z'] * las_scaleZ) + las_offsetZ)

# plotting points
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(p_X, p_Y, p_Z, c='r', marker='o')
plt.show()