What does DSM stand for?
Digital Elevation Model
Must Read:
- How to Download LiDAR Data free
- Beginner sGuide to LiDAR: Light Detection and Ranging
- How to Create Boundary shape file from LiDAR data
- Clip Raster with a Shape file in Python
- Download High Resolution Satellite Imagery free online
Install 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
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
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")
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")
Ground LiDAR Points |
Create all Surface Model in Python
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")
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")
Ground Surface/ Digital Elevation Model (DEM)/ Digital Terrain Model (DTM) |
Create Digital Surface Model in Python
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")
Digital Surface Model (DSM) |
Ground Filtering of LiDAR data in Python
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 |
Drape image over DSM in Python
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
Finally came to the article, which I was searching for so long. Thanks for sharing this article.
If you are commenting out the matplotlib inline and matplotlib notebook lines use plt.show() under each plt.title() line to show the images.
whoah this blog is wonderful i really like studying your posts.
Keep up the great work! You already know, a lot of people are looking around for this info, you
could aid them greatly.
Thanks
Hi! Someone in my Myspace group shared this site with us so I came to give it a look.
I’m definitely loving the information. I’m bookmarking and will be tweeting this to my followers!
Outstanding blog and fantastic style and design.
It’s hard to come by educated people for this subject, but you
sound like you know what you’re talking about! Thanks