3D Digital Surface Model with Python and Pylidar

from pylidar import lidarprocessor
from pylidar.toolbox import interpolation
from pylidar.toolbox import spatial
from functools import partial
import pyproj

import numpy as np

import matplotlib.pyplot as plt
import ipyvolume.pylab as p3

# To display interactive plot
# %matplotlib notebook

# To display static plot below the code block
%matplotlib inline

import cv2

from tqdm import tqdm

Download and import LiDAR data in Python

lidar_file_name = 'input1/LiDAR/USGS_LPC_WI_DodgeCo_2017_841666_LAS_2019.laz'

# Read all points of LiDAR data
data_all = spatial.readLidarPoints(lidar_file_name)

# Read only ground points of LiDAR data
data_ground = spatial.readLidarPoints(lidar_file_name,
        classification = lidarprocessor.CLASSIFICATION_GROUND)

EPSG: 26916=> NAD_1983_UTM_Zone_16N

# Input projection / coordinate referance system format
input_projection = pyproj.Proj(init='epsg:32136')
# Output projection / coordinate referance system format
output_projection = pyproj.Proj(init='epsg:4326')
transform_projection = partial(pyproj.transform, input_projection, output_projection)

# Convert foot to meter
def foot_to_meter(x):
    return x*0.3048

# Function to process
def data_conversion(p):
    xin, yin, zin = p
    xout, yout = transform_projection(foot_to_meter(xin), foot_to_meter(yin))
    return [xout, yout, foot_to_meter(zin)]

# Process all LiDAR data points
np_all = len(data_all)
x_all = np.zeros(np_all)
y_all = np.zeros(np_all)
z_all = np.zeros(np_all)

for x in tqdm(range(np_all)):
    x_all[x],y_all[x],z_all[x] = data_conversion(data_all[x])

# Process only Ground LiDAR data points
np_ground = len(data_ground)
x_g = np.zeros(np_ground)
y_g = np.zeros(np_ground)
z_g = np.zeros(np_ground)

for x in tqdm(range(np_ground)):
    x_g[x],y_g[x],z_g[x] = data_conversion(data_ground[x])

Now that we read LiDAR data into our notebook, now we need to create subset of our LiDAR data as we cannot process and plot all points because processing all LiDAR points would take huge time to process and render (to plot).

# Select randomly 1% of LiDAR data
percent_view = 0.01

num_pts_all = len(x_all)
mask_all = np.full(num_pts_all, False)
mask_all[:int(num_pts_all*percent_view)] = True
np.random.shuffle(mask_all)

x_mask_all = x_all[mask_all]
y_mask_all = y_all[mask_all]
z_mask_all = z_all[mask_all]

num_pts_g = len(x_g)
mask_ground = np.full(num_pts_g, False)
mask_ground[:int(num_pts_g*percent_view)] = True
np.random.shuffle(mask_ground)

x_mask_g = x_g[mask_ground]
y_mask_g = y_g[mask_ground]
z_mask_g = z_g[mask_ground]

Plot all LiDAR points

fig1 = plt.figure(figsize=(15,15))

ax = fig1.add_subplot(111, projection = '3d')
ax.axis('off')
plot = ax.scatter3D(xs = x_mask_all, ys = y_mask_all, zs = z_mask_all, c = z_mask_all, cmap = plt.cm.RdYlGn)
ax.set_zlim3d(250,400)
fig1.colorbar(plot)
plt.title("All Points")

All LiDAR Points

You remember, we read ground points separately. So let’s plot ground points separately.

fig2 = plt.figure(figsize = (15,15))

ax = fig2.add_subplot(111, projection = '3d')
ax.axis('off')
plot = ax.scatter3D(xs = x_mask_g, ys = y_mask_g, zs = z_mask_g, c = z_mask_g, cmap = plt.cm.RdYlGn)

ax.set_zlim3d(250,400)
fig2.colorbar(plot)
plt.title("Ground Points")

Ground LiDAR Points

binGeoSize (or in short bin size) is the geographic (i.e., in metres) size of each square bin (i.e., 1 m). Let’s assume a wall is created by multiple bricks. So 1 bin size is size of one brick.

BIN_SIZE =  1.71e-05

(xMin_all, yMax_all, ncols_all, nrows_all) = spatial.getGridInfoFromData(x_all, y_all,BIN_SIZE)
pxlCoords_all = spatial.getBlockCoordArrays(xMin_all, yMax_all, ncols_all, nrows_all, BIN_SIZE)

# All Surface model
asm = interpolation.interpGrid(x_all, y_all, z_all, pxlCoords_all, 'pynn')

Now let’s plot all surface model.

fig_surf_all = plt.figure(figsize = (15,15))

ax = fig_surf_all.add_subplot(111, projection = '3d')
ax.axis('off')

ax.plot_trisurf(x_mask_all, y_mask_all, z_mask_all)
ax.set_zlim3d(250,400)
plt.title("All Surface")

All Surface Model

Similar to all points let’s create surface model for ground. Now surface model of ground is also known as Digital Elevation Model (DEM) or Digital Terrain Model (DTM).

(xMin_g, yMax_g, ncols_g, nrows_g) = spatial.getGridInfoFromData(x_g, y_g, BIN_SIZE)
pxlCoords_g = spatial.getBlockCoordArrays(xMin_g, yMax_g, ncols_g, nrows_g, BIN_SIZE)

# Degital Elevation model (DEM)
dem = interpolation.interpGrid(x_g, y_g, z_g, pxlCoords_g, 'pynn')

Now let’s plot DEM in Python with matplotlib.

fig_surf_g = plt.figure(figsize = (15,15))

ax = fig_surf_g.add_subplot(111, projection = '3d')
ax.axis('off')
ax.plot_trisurf(x_mask_g, y_mask_g, z_mask_g)

ax.set_zlim3d(250,400)
plt.title("Ground Surface")

Ground Surface/ Digital Elevation Model (DEM)/ Digital Terrain Model (DTM)

We can subtract DEM from all surface model to get DSM.

diff_surf = asm-dem
diff_surf = np.nan_to_num(diff_surf)

temp_x = np.arange(diff_surf.shape[1])
temp_y = np.arange(diff_surf.shape[0])

temp_y = np.flip(temp_y,0)

X_diff, Y_diff = np.meshgrid(temp_x, temp_y)
Z_diff = diff_surf

Now let’s plot DSM in matplotlib.

# Plot DSM in Python

fig_dsm = plt.figure(figsize = (15,15))

ax = fig_dsm.add_subplot(111, projection = '3d')
ax.axis('off')

surf = ax.plot_surface(X_diff, Y_diff, Z_diff, rstride=1, cstride=1,linewidth=0)

ax.set_zlim3d(0,200)
plt.title("Digital Surface Model")

Digital Surface Model (DSM)

h = 1e-02
height = (Z_diff > h)

color_map = np.zeros([750, 925, 3], dtype=np.uint8)

for row_num in tqdm(range(len(height))):
    for col_num in range(len(height[row_num])):
        val = height[row_num][col_num]
        if (val==True):
            color_map[row_num, col_num, 0] = 255
            color_map[row_num, col_num, 1] = 0
            color_map[row_num, col_num, 2] = 0
        else:
            color_map[row_num, col_num, 0] = 0
            color_map[row_num, col_num, 1] = 255
            color_map[row_num, col_num, 2] = 0

# Plot Ground and non Ground
fig_classify = plt.figure(figsize = (15,15))

ax = fig_classify.add_subplot(111, projection = '3d')
ax.axis('off')

# To make RGB values in range of 0-1
drape_color = color_map/255
surf = ax.plot_surface(X_diff, Y_diff, Z_diff, rstride=1, cstride=1,facecolors=drape_color)

ax.set_zlim3d(0,200)
plt.title("Drape image over DSM")

Ground vs Non-Ground

Now we can overlay imagery over Digital Surface Model (DSM) for better visualization.

# Read clipped and transformed satellite imagery with openv
satellite_imagery = cv2.imread("input1/imagery/clipped_USGS_trans.tif")

# Overlay image over DSM
fig_drape_img = plt.figure(figsize = (15,15))

ax = fig_drape_img.add_subplot(111, projection = '3d')
ax.axis('off')

# To make RGB values in range of 0-1
drape_color=satellite_imagery / 255
surf = ax.plot_surface(X_diff, Y_diff, Z_diff, rstride=1, cstride=1,facecolors = drape_color)

ax.set_zlim3d(0,200)
plt.title("Drape image over DSM")

Drape image over DSM

Conclusion

• What is Digital Surface Model?
• What is Digital Elevation Model?
• How to read LiDAR data in pylidar in Python
• How to make subset of LiDAR points
• Plot LiDAR point cloud data in Python
• Create Surface Model in Python
• Create Digital Elevation Model in Python
• Create Digital Surface Model in Python
• Plot 3d Digital Surface model in Python