PCI Geomatics -

# Overview

This tutorial outlines how to extract a digital surface model (DSM) from Pleiades 1A imagery. There are optional steps for converting the DSM to a digital terrain model (DTM). The data used is Triscopic stereo Pleiades imagery Melbourne, Australia.

Benefits

This is a useful workflow if you are working on imagery in an area where you do not have a DEM. Using high resolution optical imagery to extract a DEM will produce a high resolution DEM, usually higher than most common DEMs. You can then use your DEM to create more accurate orthorectified images.

# Prerequisites

• Geomatica 2015 or later installed
• Python 101 (basic understanding of python coding concepts)
• Some experience programming with Geomatica (recommended)

# Tutorial

In order to generate our DSM and DTM we will run the following algorithms. In this tutorial it is assumed that GCPs have already been collected for the imagery. However, if you are adapting this script for different imagery you may wish to add the GCP collection step. In this case, please see the Orthorectification Workflow tutorial which outlines GCP collection:

You will notice in this tutorial that many of the parameters for these algorithms are left as their defaults. You can click on any of the algorithm names in this tutorial to go to the Geomatica Help and see all of the algorithm’s parameters and default values.

## 1. Import necessary modules & setup inputs/outputs

The first step is to import the Geomatica modules and setting up the variables that will point to our input and output directories and files.

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```import os from pci.rfmodel import rfmodel from pci.epipolar import epipolar from pci.autodem import autodem from pci.dsm2dtm import dsm2dtm working_dir = r"I:\Data\Tutorial\Python\DEMExtraction\MelbourneDataset" epi_dir = os.path.join(working_dir, "Epipolars") dem_dir=os.path.join(working_dir, "DEM") geocoded_dsm=os.path.join(dem_dir, "geocoded_dsm.pix") DTM = os.path.join(dem_dir, "DTM.pix") ```

The os module, which we will use to access the operating system’s file structure, is imported on line 1. Lines 2 to 5 import the four required algorithms, as outlined above. On line 7, working_dir is a variable that contains the path to the directory where your tutorial data is stored. You will need to change this to the directory on your computer where your tutorial data is saved. Two output directories are set on lines 8 and 9; one for the epipolar images and one for the DEMs (DSM and DTM). Finally, two output images are set on lines 10 and 11. Geocoded_dem will be the output for autodem, while DTM will be the output for the optional dsm2dtm step.

## 2. Calculate math model from GCPs

The algorithm rfmodel is used to calculate the mathematical model required to create epipolars.

 ```13 14 15 16``` ```# Calculate the model from the previously collected GCPs for the stereo pair rfmodel(mfile=working_dir, dbgc=[3], numcoeff=[0]) ```

On line 14 the parameter mfile is set to the working directory. Each file within the directory will be processed, and the math model calculated. The parameter dbgc is the input GCP segment which will be segment 3. Keep in mind that if you are processing multiple files, the same parameters will be used for each file. The last parameter we need to set is numcoeff which specifies the number of coefficients that will be used to calculate the rational function. Since our input images do not have any RPC segments we will set this to 0. If your input imagery has an RPC segment you will need to change this value. Check the Geomatica Help to find out more information about this variable.

## 3. Generate the epipolar images and extract the DSM

The next step will be to generate the epipolar images which will be used to triangulate the 3D coordinates for the extracted DSM.

 ``` 18 19 20 21 22 23 24 25 26 27 ``` ```epipolar(mfile=working_dir, outdir=epi_dir) # Extract a DSM from stereo pairs autodem(mfile=epi_dir, outdir=dem_dir, filedem=geocoded_dsm, pxszout=[1,1]) # Convert the extracted DSM to DTM dsm2dtm(fili=geocoded_dsm, dbic=[1], ```

Working_dir is again used as the input file. The math models that were created in the previous step for each image will be used as an input mmseg for epipolar. However, we do not have to set this parameter, since it will automatically use the last RPC segment added to the file. The epipolar directory epi_dir was defined in the first step and will be set as the value for the outdir parameter.

Once the epipolars have been generated, we can extract the dsm from the stereo pairs. On line 23 the input mfile for autodem is set to the epi_dir, where the epipolar images were saved. The outdir is set the variable dem_dir. If you want to generate a geocoded DSM you need to set the filedem parameter. In this case we set this parameter to the variable geocoded_dsm that we defined in step 1. Finally, on line 26 the output pixel size is set to 1m by 1m.

Geocoded DSM

## 4. Optional: Convert the DSM to DTM

You will notice the DSM image above includes all objects on the earth’s surface such as buildings and trees. Many applications require a DTM which represents the bare ground surface without any objects. There are two ways to generate a DTM in Geomatica.

### DEM Editing Tool

The first option is to use the DEM Editing tool in Focus. We have an html tutorial and video tutorials which explain how to use this tool to remove the effects of buildings and trees from your DSM. Specifically, the video Part 2: Working with Complex Terrain outlines how remove the effects of buildings using the terrain filter. This tool allows you to specifically edit certain areas of the DSM. There are various filters available in this tool, which you can use to edit different types of terrain. Although, this is a semi-automatic process, you can create a better looking output.

Result from DEM Editing

### DSM2DTM algorithm

The second option is to run the Geomatica algorithm dsm2dtm. This tool automatically removes the surface features from your DSM. The following section of code can be added to the end of your script to run this algorithm.

 ```28 29 30 31 32 33``` ```# Convert the extracted DSM to DTM dsm2dtm(fili=geocoded_dsm, dbic=[1], filo=DTM, dboc=[1], objsize=[300]) ```

The input for this algorithm will be the geocoded_dsm that you created in autodem. On line 30 the channel that contains the DSM is selected for the variable dbic. The output file (filo) is set to the DTM variable that we defined earlier. Since DTM.pix is a new file the output channel (dboc) will be set to 1. Finally, on line 33 the objsize is set to 300. This parameter specified the size of the filters that should be used to remove surface features. If you are using your own dataset you will likely need to test out different values for this parameter to achieve the best result. You will notice in the image below that most of the buildings are removed. You can further edit this new DTM using the DEM editing tool mentioned above.

Result from DSM2DTM