Digital Aerial Photogrammetry¶

This section provides some basic recipes for performing photogrammetric operations with digital aerial imagery, using Geomatica’s python library.

In this section you will find recipes for ingesting aerial imagery, extracting Digital Elevation Models (DEMs) from multiple aerial images in a block and orthorectification.

Create .pix Link File from Raw Aerial Images¶

This recipe demonstrates how to create a .pix link file from raw aerial imagery. This creates the .pix file format for an image, which links back to the .tif file rather than copying the actual pixels.

The pix format is highly recommended for photogrammetric processing with Geomatica functions.
from pci.link import link

raw_airphoto = r"c:\airphoto_project\raw_images\00025.tif"
pix_link_file = r"c:\airpoto_project\ingested\00025.pix"

link(fili=raw_airphoto, filo=pix_link_file)

Discover all airphotos in a directory and create .pix link file - Batch Processing¶

There are two examples below: one for discovering all valid raw airphotos in a single folder directory and the second for searching a multi-folder directory. The first example is simpler but only works if all raw image files are in a single folder.

Once all of the raw airphotos are discovered you can then run the Geomatica algorithm in a loop to process each image. In this case, link will be run to create a linked .pix file for each file in the directory.

import os, fnmatch, sys, glob
from pci.link import link

# Example 1: Valid if all raw images are in a single folder
input_files_list = glob.glob(r"C:\aerial\photos\*.pix")

# Example 2: Valid if all raw images are in a single folder or multiple folders in the directory
raw_airphoto_dir =r"C:\aerial\photos"
input_files_list = []

# r=root, d=directories, f = files
for r,d,f in os.walk(raw_airphoto_dir):
    for file in fnmatch.filter(f, "*.tif"):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r,file))

# Run link in a loop to process each image in the list - input_files_list
ingest_dir = r"C:\aerial\Ingested"
for raw_image in input_files_list:
    input_file = raw_image
    link_pix = os.path.join(ingest_dir, os.path.basename(os.path.splitext(raw_image)[0]) + ".pix")

    link(fili=input_file, filo=link_pix)

Setup Digital Airphoto Project for Photogrammetric Processing¶

This recipe demonstrates how to setup an Orthoengine project for airphoto photogrammetric processing. Running this script will result in a digital airphoto project that is ready for DEM extraction and orthorectification.

Note: this script requires the user to provide the camera calibration parameters in a formatted camera calibration xml file.

Template Camera Calibration XML file

import os
import fnmatch
from pci.link import *
from pci.crproj import *
from pci.camimport import *
from pci.eoimport import *

raw_airphoto_dir = r"D:\AerialOrtho\photos"
input_files_list = []

for r, d, f in os.walk(raw_airphoto_dir):
    for file in fnmatch.filter(f, "*.tif"):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r, file))

# Create Link file
ingest_dir = r"D:\AerialOrtho\Ingest"
for raw_image in input_files_list:
    input_file = raw_image
    link_pix = os.path.join(ingest_dir, os.path.basename(os.path.splitext(raw_image)[0]) + ".pix")

    link(fili=input_file, filo=link_pix)

# Create OrthoEngine project file
mfile = os.path.join(ingest_dir, "*.pix")
oeproj = os.path.join(ingest_dir, "oe_initial.prj")

crproj(mfile=mfile, oeproj=oeproj, modeltyp="APDU", mapunits="UTM 17 D122")

# Import Camera Parameters
camfile = os.path.join(raw_airphoto_dir, "camera_calib.xml")
prepared_oe_proj = os.path.join(ingest_dir, "oe_configured.prj")

camimport(oeproji=oeproj, tfile=camfile, oeprojo=prepared_oe_proj)

# Import Exterior Orientation File
eofile = os.path.join(raw_airphoto_dir, "eo.txt")

eoimport(oeproji=prepared_oe_proj, tfile=eofile, eoformat=[3], mapunits="UTM 17 D122")

Automatic DEM Extraction from Aerial Imagery - Batch Processing¶

This recipe demonstrates how to extract a DEM from multiple Aerial images, using an OrthoEngine project that was previously setup for digital aerial image processing (See above recipe).

import os
import fnmatch
from pci.cpmmseg import cpmmseg
from pci.epipolar import epipolar
from pci.autodem import autodem
from pci.dsm2dtm import dsm2dtm

# Create math model from OE to ingested images
ingest_airphoto_dir = r"D:\AerialOrtho\Ingest"
output_dir = r"D:\AerialOrtho\AutoDEMBatch"
prepared_oe_proj = os.path.join(ingest_airphoto_dir, "oe_configured.prj")

ingest_image_list = []
for r, d, f in os.walk(ingest_airphoto_dir):
    for file in fnmatch.filter(f, "*.tif"):
        print "found valid input file: ", file
        ingest_image_list.append(os.path.join(r, file))

cpmmseg(oeproj=prepared_oe_proj)

# Generate Epipolar Images from Stereo pairs
epi_outdir = os.path.join(output_dir, 'epipolars')
if not os.path.exists(epi_outdir):
    os.makedirs(epi_outdir)

epipolar(mfile=ingest_airphoto_dir, dbic=[1, 2, 3], selmthd="max", outdir=epi_outdir, memsize=[1024])

# Extract Digital Surface Model (DSM) from multiple stereo airphoto images and mosaic them together
geocoded_dsm = os.path.join(output_dir, "airphoto_dsm.pix")

autodem(mfile=epi_outdir, dbic=[1], demdet='high', tertype='hilly', datatype='32R', extinter=[2], demedit='yes',
        outdir=output_dir, filedem=geocoded_dsm, mapunits='UTM 17 D122', pxszout=[0.9, 0.9])

# convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
dem_file = os.path.join(output_dir, "airphoto_dem.pix")
feature_size = [250]
gradient = [30]
remove_bumps = [9, 5, 9, 5]
fill_pits = [7, 5]
smooth_dem = [9]

dsm2dtm(fili=geocoded_dsm, dbic=[1], filo=dem_file, dboc=[1], objsize=feature_size, gthresh=gradient,
        tertype="hilly", bumpflt=remove_bumps, pitflt=fill_pits, medflt=smooth_dem)

Automatic Orthorectification of Aerial Imagery¶

This recipe provides a working example of how to orthorectify a single aerial image, using an OrthoEngine project that was previously setup for digital aerial image processing (See above recipe).

from pci.cpmmseg import cpmmseg
from pci.ortho import ortho

#Create math model from OE to ingested images
ingest_airphoto = r"C:\Aerial\Ingested\A17563-001.pix"
prepared_oe_proj = r"C:\Aerial\Ingested\oe_configured.prj"

# Copies the math model from the project file into a math model segment in FILE
cpmmseg(oeproj=prepared_oe_proj, file=ingest_airphoto)

#Orthorectify imagery with a user provided DEM
ortho_airphoto = r"C:\Aerial\Orthos\A17563-001.pix"
dem = r"C:\Aerial\DEM\1962_Saskatchewan_DEM.pix"

ortho(mfile=ingest_airphoto, filo=ortho_airphoto, ftype="PIX", bxpxsz="0.2", bypxsz="0.2", filedem=dem)

Automatic Orthorectification of Aerial Imagery - Batch Processing¶

This recipe provides a working example of how to orthorectify multiple aerial images, using an OrthoEngine project that was previously setup for digital aerial image processing (See above recipe).

import os, fnmatch
from pci.cpmmseg import cpmmseg
from pci.ortho import ortho

ingest_airphoto_dir = r"C:\Aerial\Ingested"
prepared_oe_proj = os.path.join(ingest_airphoto_dir, "oe_configured.prj")

# Copy math model from prepared orthoengine project file to all ingested images
# If the FILE parameter is not specified, all images contained within the project are copied.
cpmmseg(oeproj=prepared_oe_proj)

# Orthorectify imagery with a user provided DEM
ortho_dir = r"C:\Aerial\Orthos"
dem = r"C:\Aerial\dtm\PCI_Hamilton14_DTM_45cm.pix"

ortho(mfile=ingest_airphoto_dir, filo=ortho_dir, ftype="PIX", bxpxsz="0.2", bypxsz="0.2", filedem=dem)

Digital Aerial Orthorectification with user provided DEM - Full Workflow¶

This is a complete orthorectification workflow for digital aerial imagery (not including DEM extraction). This script combines multiple concepts from the previous recipes in this section.

import os
import fnmatch
from pci.link import *
from pci.crproj import *
from pci.camimport import *
from pci.eoimport import *
from pci.cpmmseg import cpmmseg
from pci.ortho import ortho

raw_airphoto_dir = r"D:\AerialOrtho\photos"

input_files_list = []

for r, d, f in os.walk(raw_airphoto_dir):
    for file in fnmatch.filter(f, "*.tif"):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r, file))

# Create Link file
ingest_dir = r"D:\AerialOrtho\ingest"
link_image_list = []
for raw_image in input_files_list:
    link_pix = os.path.join(ingest_dir, os.path.basename(os.path.splitext(raw_image)[0]) + ".pix")

    link(fili=raw_image, filo=link_pix)
    link_image_list.append(raw_image)

# Create OrthoEngine project file
mfile = os.path.join(ingest_dir, "*.pix")
oeproj = os.path.join(ingest_dir, "oe_initial.prj")

crproj(mfile=mfile, oeproj=oeproj, modeltyp="APDU", mapunits="UTM 17 D122")

# Import Camera Parameters
camfile = r"D:\AerialOrtho\camera\camera_calib.xml"
prepared_oe_proj = os.path.join(ingest_dir, "oe_configured.prj")

camimport(oeproji=oeproj, tfile=camfile, oeprojo=prepared_oe_proj)

# Import Exterior Orientation File
eofile = r"D:\AerialOrtho\eo\eo.txt"

eoimport(oeproji=prepared_oe_proj, tfile=eofile, eoformat=[3], mapunits="UTM 17 D122")

# Copy math model from prepared orthoengine project file to all ingested images
cpmmseg(oeproj=prepared_oe_proj)

# Orthorectify imagery with a user provided DEM
ortho_dir = r"D:\AerialOrtho\orthos"
dem = r"D:\AerialOrtho\dtm\PCI_Hamilton6_DTM_45cm.pix"

ortho(mfile=ingest_dir, filo=ortho_dir, ftype="PIX", bxpxsz="0.2", bypxsz="0.2", filedem=dem)

Digital Aerial Orthorectification with DEM Extraction - Full Workflow¶

This is a complete orthorectification workflow for digital aerial imagery (including DEM extraction). This script combines multiple concepts from the previous recipes in this section.

import os
import fnmatch
from pci.link import *
from pci.crproj import *
from pci.camimport import *
from pci.eoimport import *
from pci.cpmmseg import cpmmseg
from pci.epipolar import epipolar
from pci.autodem import autodem
from pci.dsm2dtm import dsm2dtm
from pci.ortho import ortho

raw_airphoto_dir = r"D:\AerialOrtho\photos"

input_files_list = []

for r, d, f in os.walk(raw_airphoto_dir):
    for file in fnmatch.filter(f, "*.tif"):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r, file))

# Create Link file
ingest_dir = r"D:\AerialOrtho\Ingest"
link_image_list = []
for raw_image in input_files_list:
    link_pix = os.path.join(ingest_dir, os.path.basename(os.path.splitext(raw_image)[0]) + ".pix")

    link(fili=raw_image, filo=link_pix)
    link_image_list.append(raw_image)

# Create OrthoEngine project file
mfile = os.path.join(ingest_dir, "*.pix")
oeproj = os.path.join(ingest_dir, "oe_initial.prj")

crproj(mfile=mfile, oeproj=oeproj, modeltyp="APDU", mapunits="UTM 17 D122")

# Import Camera Parameters
camfile = os.path.join(raw_airphoto_dir, "camera_calib.xml")
prepared_oe_proj = os.path.join(ingest_dir, "oe_configured.prj")

camimport(oeproji=oeproj, tfile=camfile, oeprojo=prepared_oe_proj)

# Import Exterior Orientation File
eofile = os.path.join(raw_airphoto_dir, "eo.txt")

eoimport(oeproji=prepared_oe_proj, tfile=eofile, eoformat=[3], mapunits="UTM 17 D122")

# Copy math model from prepared orthoengine project file to all ingested images
cpmmseg(oeproj=prepared_oe_proj)

#Generate Epipolar Images from Stereo pairs
epi_outdir = os.path.join(ingest_dir, 'epipolars')
if not os.path.exists(epi_outdir):
    os.makedirs(epi_outdir)

epipolar(mfile=ingest_dir, dbic=[3], selmthd="max", sampling=[1], outdir=epi_outdir)

#Extract Digital Surface Model (DSM) from multiple stereo airphoto images and mosaic them together
outdir = r"D:\AerialOrtho\OrthoDEM"
if not os.path.exists(outdir):
    os.makedirs(outdir)

geocoded_dsm = os.path.join(outdir, "airphoto_dsm.pix")

autodem(mfile=epi_outdir, outdir=outdir, filedem=geocoded_dsm, pxszout=[0.9, 0.9])

#convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
dem_file = os.path.join(outdir, "airphoto_dem.pix")

dsm2dtm(fili=geocoded_dsm, dbic=[1], filo=dem_file, dboc=[1], objsize=[150], gthresh=[35], tertype="hilly",
        bumpflt=[9, 5, 9, 5], pitflt=[7, 5], medflt=[9])

#Orthorectify imagery with a user provided DEM

ortho(mfile=ingest_dir, filo=outdir, ftype="PIX", bxpxsz="0.2", bypxsz="0.2", filedem=dem_file)

Create Camera Calibration (.xml) file¶

This recipe demonstrates how to quickly create a camera calibration file in the XML format required by Geomatica’s automated aerial processing algorithms.

In most cases it will be easier to acquire a template of the XML file (See link below) and edit it manually. However, this script provides an example of a logical process for creating a GUI that accepts the user inputs required for creating the camera calibration xml file in a user-friendly fashion; it limits the input fields based on the sensor type.