**What does DSM stand for?**

**how to create DSM and DEM with python**.

**Digital Elevation Model**

**Drape image over 3d DSM 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**

**Python library named Pylidar to process LiDAR data**. You can check Pylidardocumentation to know about

**Pylidar**in detail.

**conda install -c rios pylidar**

**pip install pyproj==1.9.6**

**pip install matplotlib**

**pip install tqdm**

**pip install opencv-python**

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

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**

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

**All surface model contains everything on earth (including ground)**

__Note:__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**

**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 vs Non-Ground |

**Drape image over DSM in Python**

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

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

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

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