3D Digital Surface Model with Python and Pylidar

Digital Elevation model in python

What does DSM stand for?

DEM vs DSM

Digital Surface Model (DSM) represents the top of the earth’s surface. It includes TREES, BUILDINGS and other objects that sit on the earth.

DSM is useful by creating 3D model for telecommunications, urban planning, aviation etc.

In this post I will show you how to create DSM and DEM with python.
Now what is DEM?

Digital Elevation Model

Digital Elevation Model (DEM) also known as Digital Terrain Model (DTM). It represents elevation of ground. It does not include other objects of earth like buildings, tree etc.

In this post I will cover below points:
·       Create Digital surface model (DSM) with Python
·       Plot DSM with Python
·       Create Digital Elevation model (DEM) with Python
·       Plot DEM with Python
·       Visualize Ground and Non-Ground DSM
·       Drape image over 3d DSM model

Install Packages

In this post I will use Python library named Pylidar to process LiDAR data. You can check Pylidardocumentation to know about Pylidar in detail.

To follow this article you have to install some packages. Listing commands to install those packages.

conda install -c rios pylidar
pip install pyproj==1.9.6
pip install matplotlib
pip install tqdm
pip install opencv-python

It is a good practice to make a fresh virtualenvironment while working with this kind of project.


Now let’s import all required packages

Import packages

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

You can download LiDAR data for your analysis area (Area of Interest) or just download LiDAR data I have used in this article.
You need to import LiDAR data after downloading. We will import LiDAR data into two files separately.
1.    Import all LiDAR data
2.    Import only ground points
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)

Change Projection System of LiDAR data in Python

To work with imagery data, we need to change Coordinate Reference System (CRS) of LiDAR data. The imagery and LiDAR data should have same projection.

In addition to that Pylidar and other libraries expect LiDAR data in meter not foot. So we need to convert LiDAR data from foot to meter.

Sharing details of projected coordinate system used in this article for your reference.

EPSG:32136=> NAD_1983_StatePlane_Tennessee_FIPS_ 4100
EPSG:4326=> WGS84

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])

Make Subset of LiDAR points

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")

plot lidar data in matplotlib
All LiDAR Points

Plot Only Ground 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")

plot 3d lidar data in python
Ground LiDAR Points

Create all Surface Model in Python

To create Surface Model (DSM) we need to convert LiDAR data point to raster surface.
Let’s first convert all LiDAR points to raster surface. To crate surface model in Pylidar setting bin size is required.

What is Bin size?

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')

Plot all Surface Model in Python

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")

plot surface model in matplotlib
All Surface Model

Create Digital Elevation Model in Python

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')

Plot Digital Elevation Model in Python

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")

plot DEM in python plot DTM in Python
Ground Surface/ Digital Elevation Model (DEM)/ Digital Terrain Model (DTM)

Create Digital Surface Model in Python

As you know Digital Surface Model (DSM) contains top of the earth’s surface like tree, buildings etc. and Digital Elevation Model (DEM) or Digital Terrain Model (DTM) represents elevation of ground.

Note:All surface model contains everything on earth (including ground)

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

Plot Digital Surface Model in Python

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")

plot DSM in python
Digital Surface Model (DSM)

Ground Filtering of LiDAR data in Python

After creating Digital Surface model (DSM) now let’s segment or classify ground and non-ground by different colour.

Considering height above 1e-02 meter as non-ground.
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 filtering of lidar data in Python
Ground vs Non-Ground

Drape image over DSM in Python

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

To drape overlay image over DSM you need to follow below steps:
3.    Transform projection of imagery same to LiDAR data
If you want to avoid above steps and just want to follow this tutorial, you can download clipped and transformed imagery used for this tutorial.

Note: We have changed projection of LiDAR data to EPSG: 4326 so projection of clipped and transformed imagery must be EPSG: 4326. Ensure this projection before draping imagery over DSM in python.
You can check my previous article to check and transform projection system in Python with Rasterio Library.

# 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")

overlay image over DSM in python
Drape image over DSM

Conclusion

In this post you learned
  • 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

If you have any question or suggestion regarding this topic see you in comment section. I will try my best to answer.

1 thought on “3D Digital Surface Model with Python and Pylidar”

Leave a Comment

Your email address will not be published. Required fields are marked *