Satellite Photogrammetry

This section provides some basic recipes for photogrammetric processing of satellite imagery using Geomatica’s Py library. Topics include: Ingest raw satellite imagery, automatic GCP collection, automatic tie point collection, automatic refinement of GCPs and Tie Points, DSM extraction, DSM-to-DEM filtering, orthorectification, batch processing, advanced conditions and more...

Ingest Older Raw Satellite Images

This recipe demonstrates how to ingest an image from a Formosat-2 sensor. The use of specific CD* functions to ingest data is only required for older sensors. However, running the CD* functions for older sensors is required before the images can be viewed or processed

from pci.cdformosat import cdformosat

in_file = r"C:\Formosat\FS2_G025_MS_L1A_20071024_021110.tif"
out_file = r"C:\Formosat\Ingest\Formosat.pix"
linkfile = u"NO"
channels = [1,2,3,4]

cdformosat(in_file, out_file, linkfile, channels, u"")

Import Satellite Images from New Sensors

This recipe demonstrates how to import a satellite image from a newer satellite sensor into a .pix file. This is NOT required for viewing the imagery and performing many types of processes on the data, but is often useful for photogrammetric processing.

from pci.fimport import fimport

raw_kompsat2_image = r"C:\Kompsat\MSC_080407140610_09045_24220801PMSN14_1G.eph"
imported_image = r"C:\Kompsat\Ingest\kompsat2_image.pix"

fimport(raw_kompsat2_image, imported_image, [], "", "TILED512")
print "Complete"

Import satellite images from auxiliary files

Some raw satellite images are loaded with auxiliary text or xml files, such as, Landsat-8 and Kompsat-3 images. This concept can sometimes confuse scripters when it comes to importing the imagery through code.

This recipe demonstrates how to ingest panchromatic, multispectral and thermal imagery from landsat-8 imagery, using the auxiliary text file.

from pci.fimport import fimport

#import landsat-8 multispectral imagery
raw_ms_l8 = r"C:\LC80150292013111LGN01\LC80150292013111LGN01_MTL.txt-MS" #notice the -MS at the end
pix_ms_l8 = r"C:\ingest\landsat8_ms.pix"

fimport(raw_ms_l8, pix_ms_l8, [], "", "TILED512")


#import landsat-8 panchromatic imagery
raw_pan_l8 = r"C:\LC80150292013111LGN01\LC80150292013111LGN01_MTL.txt-PAN" #notice the -PAN at the end
pix_pan_l8 = r"C:\ingest\landsat8_pan.pix"

fimport(raw_pan_l8, pix_pan_l8, [], "", "TILED512")


#import landsat-8 thermal imagery
raw_therm_l8 = r"C:\LC80150292013111LGN01\LC80150292013111LGN01_MTL.txt-Thermal" #notice the -Thermal at the end
pix_therm_l8 = r"C:\ingest\landsat8_therm.pix"

fimport(raw_therm_l8, pix_therm_l8, [], "", "TILED512")

Automatic GCP Collection and Refinement

This recipe demonstrates how to automatically collect GCPs on a raw satellite image and refine them based on statistical accuracy. In this case, root mean square error (RMSE) will be used to remove blunder GCPs.

from pci.link import link
from pci.autogcp2 import *
from pci.gcprefn import *

#Create Link File to hold new GCP segment(s)
raw_geoeye_image = r"C:\geoeye\po_311159_0000000\po_311159_bgrn_0000000.tif"
link_pix_file = r"C:\geoeye\Output\po_311159_bgrn_0000000.pix"

link(raw_geoeye_image, link_pix_file, [] )

#Run GCP Collection and store new GCP segment in link pix file
raw_match_chan = [3]       #When available NIR channel is recommended
reference_image = r"c:\geoeye\Ref\S-55-40_lr_2000.tif"
ref_match_chan = [3]
dem = r"C:\geoeye\po_311159_0000000\dem\gmted2010.pix"
gcp_search_radius = [50]

gcp_seg = autogcp2(link_pix_file, raw_match_chan, u"", [], [], reference_image, ref_match_chan, dem,
                [], [], u"", u"", [], u"", u"", gcp_search_radius, u"", [], u"", [])
#autogcp2 function returns the new gcp segment, which can then be directly passed to the gcp refinement function

#Setup required parameters for automatic GCP refinement
model_type = u"RF"    #rational function model
order = [0]           #a polynomial order of 0 is best for most DigitalGlobe sensors (due to very accurate RPCs)
reject = [5,0.5,0.5]  #Continue to reject worse GCPs until RMSE falls below 0.5 pixel in x and y directions

gcprefn(link_pix_file, gcp_seg, [], [], model_type, order, reject, u"NO", [], [])

#accurate GCPs should now be stored in the link pix file and can be used for further photogrammetric processing

Automatic GCP Collection and Refinement - Batch Processing

This recipe demonstrates how to automatically discover all valid images in a folder and subfolders, import the images as a pci link file, and collect and refine GCPs on all discovered images.

import fnmatch, os, sys
from pci.link import link
from pci.autogcp2 import *
from pci.gcprefn import *

raw_dir = r"c:\geoeye\RAWDATA"
file_filter = '*_bgrn_*.tif' #Looks for the Multispectral image only
input_files_list = []

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

if len(input_files_list) > 0: #handles condition if the input list is empty
    for raw_file in input_files_list: #iterates through each image file found.
        #Create .pix link file to store the GCP segement that holds the GCPs
        ingest_link_dir = r"C:\geoeye\output\ingest"
        ingest_link_filename = os.path.splitext(os.path.basename(raw_file))[0] + ".pix" #gets the filename of the raw image (without extension)
        ingest_link_file = ingest_link_dir + os.sep + ingest_link_filename

        link(raw_file, ingest_link_file, [] )

        #Perform automatic GCP Collection
        dbic = [1]      # Channel from input raw file that is used for automatic GCP Collection
        filref = r"c:\geoeye\Ref\S-55-40_lr_2000.tif"
        filedem = r"C:\geoeye\dem\gmted2010.pix"
        searchr = [50]      # Search radius. Default is to search in pixels

        gcp_seg = autogcp2(ingest_link_file, dbic, u"", [], [], filref, [], filedem, [],
                           [], u"", u"", [], u"", u"", searchr, u"", [], u"", [])

        #Refine the GCPs based on RMSE
        model_type = u"RF"
        order = [0]
        reject = [5,0.5,0.5]

        gcprefn(ingest_link_file, gcp_seg, [], [], model_type, order, reject, u"NO", [], [])

else:
    print "No image files found in " + root_dir + " that match your filter criteria."
    sys.exit()

Iterate through GCP collection/refinement until success - Batch Processing

This recipe provides a basic demonstration of how one might create an advanced workflow for intelligently running multiple scenarios with different GCP collection and refinement settings until acceptable results are achieved. Specifically, this recipe runs a maximum of 10 iterations of the GCP collection, each time increasing the search radius until more than 10 GCPs are found or the search radius exceeds 500 pixels.

#Needs work. Condition for when it fails to find anything.

import fnmatch, os, sys
from pci.autogcp2 import *
from pci.gcprefn import *
from pci.link import *

root_dir = r"D:\Kompsat-2\Raw_Scenes_1G"
file_filter = '*G_1G.eph' #Looks for the multispectral Kompsat2 image
input_files_list = []

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

if len(input_files_list) > 0:
    for raw_file in input_files_list:
        ingest_file = r"D:\Kompsat-2\ingested" + os.sep + os.path.splitext(os.path.basename(raw_file))[0] + ".pix"
        link(raw_file, ingest_file, [])
        increase_radius = 0
        numgcps = [0] #It is required to initially set the number of GCPs to zero to satisfy the condition

        while ((int(numgcps[0]) > 9) == True) !=  ((increase_radius < 500) == True): #Mimicking an XOR condition with two expressions normalized as boolean
            #Setup required parameters for automatic GCP Collection
            gcp_match_chan = [1]
            filref = r"D:\Kompsat-2\Reference\SPOT\o2_5m.pix"
            filedem = r"D:\Kompsat-2\DEM\torontodem.pix"
            search_radius = [50 + increase_radius] #Search radius. Increase the search radius by 50 pixels each iteration
            numgcps = [] #AUTOGCP2 will return the number of gcps it found in this list

            gcp_seg = autogcp2(ingest_file, dbic, u"", [], [], filref, [], filedem, [],
                     [], u"", u"METER", [], u"", u"", search_radius, u"", [], u"", numgcps)

            increase_radius += 50 #adds 50 pixels to the search radius for the next iteration (If required)

        #Run automatic GCP refinement (remove bad gcps)
        model_type = u"RF"
        order = [0]
        reject = [5,1,1]

        gcprefn(ingest_file, gcp_seg, [], [], model_type, order, reject, u"NO", [], [])

else:
    print "No image files found in " + root_dir + " that match your filter criteria."
    sys.exit()

Automatic GCP Collection and refinement from road vector layer

This recipe demonstrates how to automatically collect GCPs from a road vector layer on a raw input image

import os
from pci.fftmvec import *
from pci.imageinv import *
from pci.fimport import *
from pci.iii import *
from pci.pcimod import *
from pci.gcprefn import *

#Import the raw Kompsat-2 image to .pix file
raw_kompsat2_file = r"C:\Kompsat-2\Raw_Scenes_1G\MSC_070611154402_04650_23191328_1G\MSC_070611154402_04650_23191328M4P17R_1G.eph"
imported_k2_image = r"C:\Kompsat-2\working\kompsat2.pix"

fimport(raw_kompsat2_file, imported_k2_image, [], "", "TILED512")

#Invert the Near-Infrared channel - Helps FFTMVEC to successfully collect GCPs
inverse_image = r"C:\Kompsat-2\working\inv_k2.pix"
inv_channel = [3] #Near-Infrared channel

imageinv(imported_k2_image, inverse_image, inv_channel, [], [])

#copy the inverted image channel into the ingested image we are going to collect GCPs for
pcimod(imported_k2_image, "ADD", [0,0,0,1])
iii(inverse_image, imported_k2_image, [1], [5], [], [])

#Run FFTMVEC to automatically collect GCPs from road vectors
roadvec = r"C:\Kompsat-2\Reference\roads\toronto_vector.pix"
gcp_seg = [] #fftmvec provides the output gcp segment number to this parameter
fftmvec(imported_k2_image, [1,2,4], [5], [], roadvec, [2], [], [15], "", gcp_seg)

#Refine GCPs - remove bad GCPs until RMSE falls below 1 pixel in x and y direction
dbgc = gcp_seg
modeltyp = u"RF"
order = [0]
reject = [5,1,1]

gcprefn(imported_k2_image, dbgc, [], [], modeltyp, order, reject, u"NO", [], [])

#Clean up unnecessary files and channels
pcimod(imported_k2_image, "DEL", [5])
os.remove(inverse_image)

Subset an OrthoEngine Project

This recipe demonstrates how to split an existing OrthoEngine project into smaller subsets for improved processing

This can be useful when performing a bundle adjustment on big projects. Breaking big projects into smaller subsets and performing GCP and tie point collection independently on each subset can improve bundle adjustment performance and in some cases accuracy. The projects can then be merged together

#UNDER CONSTRUCTION

Merge multiple OrthoEngine projects with tie points

This recipe demonstrates how to merge multiple OrthoEngine projects. It also performs automatic tie point collection between the images at the connection points between the OrthoEngine project being merged to improve coregistration between them.

This can be useful when performing a bundle adjustment on big projects. Breaking big projects into smaller subsets and performing GCP and tie point collection independently on each subset can improve bundle adjustment performance and in some cases, accuracy.

#UNDER CONSTRUCTION

Satellite Orthorectification with GCPs previously collected - Batch Processing

This recipe demonstrates how to use GCPs that were previously collected to improve the accuracy of a math model (RPCs), which is then used to orthorectify the imagery

from pci.rfmodel import *
from pci.ortho2 import *

#Calculate the mathematical model from GCPs (rational function model) required to orthorectify the imagery
ingest_dir = r"D:\Python_Cookbook\Kompsat-2\ingested\*.pix"
gcp_seg = [9] #use GCPs from gcp segment 9 in the ingest files

rfmodel( ingest_dir, gcp_seg, [], [1], "" )

#Orthorectify the imagery with the new and improved math model (RPC segment)
ortho_dir = r"D:\Python_Cookbook\Kompsat-2\orthos"
dem = r"D:\Python_Cookbook\Kompsat-2\DEM\torontodem.pix"
mapunits = ""

ortho2( ingest_dir, [], [], [], "", ortho_dir, "tif", "", [], "", "", "", "", [],
       "", mapunits, "4", "4", dem, [], [], "", "", [], "", [], "" )

Automatic Bundle Adjustment and Orthorectification (GCPs previously collected)

This recipe demonstrates how to automatically collect tie points from overlapping regions of images and perform a block bundle adjustment to ensure relative alignment between images

import os
from pci.autotie2 import *
from pci.tprefn import *
from pci.crproj import *
from pci.cpmmseg import *
from pci.ortho2 import *

#Create OrthoEngine Project, which is required for the bundle adjustment
image_dir = r"C:\ompsat\ingest"
work_dir = r"C:\kompsat\working\TiePoints"
dbgc = [4] #Import GCPs from GCP segment 4 (see GCP collection recipes)
oeproj = work_dir + os.sep + "import_oemodel.prj"
projection = "UTM 17 T D000" #D000 is WGS84

crproj(image_dir, dbgc, oeproj, u"RFI1", projection)

#Run Automatic Tie Point Collection
tp_oeproj = work_dir + os.sep + "output_oetp.prj"
match_chan = [4] #NIR channel is often yeilds best matching results
dem = r"C:\dem\gmted2010.pix"

autotie2( oeproj, tp_oeproj, "", "", match_chan, dem, [], [], "", "", [], "", "", "", [], "", [], "" )

#Run automatic Tie Point Refinement
reject = [5, 1.0, 1.0]

tprefn(tp_oeproj, reject, "")

#Copy Math Model out to file
cpmmseg(tp_oeproj, "", [])

#orthorectify images
ortho_dir = r"C:\Kompsat\orthos"
res = ["4", "4"]
format = "TIF"

ortho2( image_dir, [], [], [], "", ortho_dir, format, "", [], "", "", "", "", [], "", "",
        res[0], res[1], dem, [], [], "", "", [], "", [], "" )

Automatic DEM Extraction with Model Refinement(GCPs previously collected)

This recipe demonstrates how to automatically extract a Digital Surface Model (DSM) from satellite imagery and automatically filter out the surface features to create a bare-earth Digital Terrain Model (DTM)

from pci.rfmodel import rfmodel
from pci.epipolar2 import epipolar2
from pci.autodem2 import autodem2
from pci.dsm2dtm import dsm2dtm

#Update the geometric accuracy of the input stereo pair based on the GCPs (previously collected)
ingest_dir = r"C:\geoeye\ingest"
gcp_seg = [4]

rfmodel(ingest_dir, gcp_seg, [], [0], "")

#Generate the epipolar images, which are used to triangulate the 3D coordinates for the DEM
pairing_method = "max"
epi_dir = r"C:\geoeye\epipolars"

epipolar2(ingest_dir, gcp_seg, [], [], "", "max", [], [1], epi_dir, [], "", "", [])

#Automatically extract a Digital Surface Model from satellite stereo pairs
out_dir = r"C:\geoeye\working\dsms"
geocoded_dsm = r"c:\geoeye\dsm_new\geoeye_dsm.pix"
dem_res = [ 6, 6 ]

autodem2(epi_dir, [], [], [], [], [], "", "", "", [], "", "", "", "", out_dir,
         geocoded_dsm, "", "", [], [], dem_res, "")

#convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
dem_file = r"C:\geoeye\dsm_new\geoeye_dem.pix"
feature_size = [150]
gradient = [35]
remove_bumps = [9,5,9,5]
fill_pits = [7,5]
smooth_dem = [9]

dsm2dtm( geocoded_dsm, [1], dem_file, [1], [], [], "", feature_size, gradient, "hilly",
         remove_bumps, fill_pits, smooth_dem, [] )

Satellite Photogrammetric Batch Workflow with user provided DEM - Full Workflow

This recipe is a complete batch satellite photogrammetric workflow that performs the following: Discovers and ingests multiple raw satellite images, automatically collects GCPs from a reference image, refines GCPs based on their statistical accuracy, automatically collects tie points and refines them based on statistical accuracy, performs a block bundle adjustment and finally, orthorectifies the imagery.

import fnmatch, os, sys
from pci.link import *
from pci.autogcp2 import *
from pci.gcprefn import *
from pci.autotie2 import *
from pci.tprefn import *
from pci.crproj import *
from pci.cpmmseg import *
from pci.ortho2 import *
from pci.exceptions import *

raw_dir = r"C:\geoeye\RAWDATA"
file_filter = '*_bgrn_*.tif' #Looks for the Multispectral image only

input_files_list = []
for r,d,f in os.walk(raw_dir):
    for file in fnmatch.filter(f, file_filter):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r,file))

if len(input_files_list) > 0:
    for raw_file in input_files_list:

        #Create link files of raw images so that they can store GCP segments
        ingest_link_dir = r"C:\geoeye\ingest"
        ingest_link_filename = os.path.splitext(os.path.basename(raw_file))[0] + ".pix" #gets the filename of the raw image (without extension)
        ingest_link_file = ingest_link_dir + os.sep + ingest_link_filename

        if not os.path.exists(ingest_link_dir):
            os.makedirs(ingest_link_dir)

        try:
            link( raw_file, ingest_link_file, [] )
        except PCIException, e:
            print e
        except Exception, e:
            print e

        #Setup required parameters for automatic GCP Collection
        in_channel = [1]      # Channel from input raw file that is used for automatic GCP Collection
        ref_image = r"C:\geoeye\Ref\S-55-40_lr_2000.tif"
        dem = r"C:\geoeye\dem\gmted2010.pix"
        searchr = [50]      # Search radius. Default is to search in pixels

        try:
            gcp_seg = autogcp2(ingest_link_file, in_channel, u"", [], [], ref_image, [], dem, [],
                            [], u"", u"", [], u"", u"", searchr, u"", [], u"", [])
        except PCIException, e:
            print e
        except Exception, e:
            print e


        #Setup required parameters for automatic GCP refinement
        model_type = u"RF" #Rational Function Model
        order = [0]
        reject = [5,2.0,2.0]     # Continue to reject worse GCPs until RMSE falls below 0.5 pixel in x and y directions

        try:
            gcprefn(ingest_link_file, gcp_seg, [], [], model_type, order, reject, u"NO", [], [])
        except PCIException, e:
            print e
        except Exception, e:
            print e

else:
    print "No image files found in " + root_dir + " that match your filter criteria."
    sys.exit()

#Create OrthoEngine Project, which is required for the bundle adjustment
work_dir = r"C:\geoeye\working"
if not os.path.exists(work_dir):
    os.makedirs(work_dir)

dbgc = [gcp_seg[0] + 1]
oeproj = work_dir + os.sep + "import_oemodel.prj"
projection = "UTM 55 G D000" #D000 is WGS84

crproj(ingest_link_dir, dbgc, oeproj, "RFI1", projection)

#Run Automatic Tie Point Collection
tp_oeproj = work_dir + os.sep + "output_oetp.prj"
tp_match_chan = [4] #NIR channel often yields best matching results

autotie2( oeproj, tp_oeproj, "", "", tp_match_chan, dem, [], [], "", "", [], "", "", "", [], "", [], "" )

#Run automatic Tie Point Refinement
reject = [5, 0.5, 0.5]

tprefn(tp_oeproj, reject, "")

#Copy Math Model out to file
cpmmseg(tp_oeproj, "", [])

#orthorectify images
ortho_dir = r"C:\geoeye\orthos"
if not os.path.exists(ortho_dir):
    os.makedirs(ortho_dir)

res = ["2", "2"]

ortho2( ingest_link_dir, [], [], [], "", ortho_dir, "TIF", "", [], "", "", "", "", [], "", "",
        res[0], res[1], dem, [], [], "", "", [], "", [], "" )

Satellite Photogrammetric Batch Workflow with DEM Extraction - Full Workflow

This recipe is a complete batch satellite photogrammetric workflow that performs the following: Discovers and ingests multiple raw satellite images, automatically collects GCPs from a reference image, refines GCPs based on their statistical accuracy, automatically collects tie points and refines them based on statistical accuracy, performs a block bundle adjustment, performs DEM extraction and finally, orthorectifies the imagery.

import fnmatch, os, sys, shutil
from pci.link import *
from pci.autogcp2 import *
from pci.gcprefn import *
from pci.autotie2 import *
from pci.tprefn import *
from pci.crproj import *
from pci.cpmmseg import *
from pci.epipolar2 import *
from pci.autodem2 import *
from pci.dsm2dtm import *
from pci.ortho2 import *
from pci.exceptions import *

raw_dir = r"C:\geoeye\RAWDATA"
file_filter = '*_bgrn_*.tif' #Looks for the Multispectral image only

input_files_list = []
for r,d,f in os.walk(raw_dir):
    for file in fnmatch.filter(f, file_filter):
        print "found valid input file: ", file
        input_files_list.append(os.path.join(r,file))

if len(input_files_list) > 0:
    for raw_file in input_files_list:

        #Create link files of raw images so that they can store GCP segments
        ingest_link_dir = r"C:\geoeye\ingest"
        ingest_link_filename = os.path.splitext(os.path.basename(raw_file))[0] + ".pix" #gets the filename of the raw image (without extension)
        ingest_link_file = ingest_link_dir + os.sep + ingest_link_filename

        if not os.path.exists(ingest_link_dir):
            os.makedirs(ingest_link_dir)

        try:
            link( raw_file, ingest_link_file, [] )
        except PCIException, e:
            print e
        except Exception, e:
            print e

        #Setup required parameters for automatic GCP Collection
        in_channel = [1]      # Channel from input raw file that is used for automatic GCP Collection
        ref_image = r"C:\geoeye\Ref\S-55-40_lr_2000.tif"
        dem = r"C:\geoeye\dem\gmted2010.pix"
        searchr = [50]      # Search radius. Default is to search in pixels

        try:
            lasc = autogcp2(ingest_link_file, in_channel, u"", [], [], ref_image, [], dem, [],
                            [], u"", u"", [], u"", u"", searchr, u"", [], u"", [])
        except PCIException, e:
            print e
        except Exception, e:
            print e


        #Setup required parameters for automatic GCP refinement
        gcp_seg = lasc
        modeltyp = u"RF"
        order = [0]      # For rational function models (RF), you need to specify a polynomial order to use for calculating error
        reject = [5,2.0,2.0]     # Continue to reject worse GCPs until RMSE falls below 0.5 pixel in x and y directions

        try:
            gcprefn(ingest_link_file, gcp_seg, [], [], modeltyp, order, reject, u"NO", [], [])
        except PCIException, e:
            print e
        except Exception, e:
            print e

else:
    print "No image files found in " + root_dir + " that match your filter criteria."
    sys.exit()

#Create OrthoEngine Project, which is required for the bundle adjustment
work_dir = r"C:\geoeye\Working\bundle"
if not os.path.exists(work_dir):
    os.makedirs(work_dir)

gcp_seg = [lasc[0] + 1]
oeproj = work_dir + os.sep + "import_oemodel.prj"
projection = "UTM 55 G D000" #D000 is WGS84

crproj(ingest_link_dir, gcp_seg, oeproj, "RFI1", projection)

#Run Automatic Tie Point Collection
tp_oeproj = work_dir + os.sep + "output_oetp.prj"
tp_match_chan = [4] #NIR channel is often yields best matching results

autotie2( oeproj, tp_oeproj, "", "", tp_match_chan, dem, [], [], "", "", [], "", "", "", [], "", [], "" )

#Run automatic Tie Point Refinement
reject = [5, 0.5, 0.5]

tprefn(tp_oeproj, reject, "")

#Copy Math Model out to file
cpmmseg(tp_oeproj, "", [])

#Generate the epipolar images, which are used to triangulate the 3D coordinates for the DEM
pairing_method = "max"
epi_dir = r"C:\geoeye\epipolars"
if not os.path.exists(epi_dir):
    os.makedirs(epi_dir)

epipolar2(ingest_link_dir, [3], [], [], "", "max", [], [1], epi_dir, [], "", "", [])

#Automatically extract a Digital Surface Model from satellite stereo pairs
epi_dem_dir = r"C:\geoeye\epi_dems"
if not os.path.exists(epi_dem_dir):
    os.makedirs(epi_dem_dir)

geocoded_dsm = r"C:\geoeye\dsm_new\geoeye_dsm.pix"
if not os.path.exists(os.path.dirname(geocoded_dsm)):
    os.makedirs(os.path.dirname(geocoded_dsm))

dem_res = [ 6, 6 ]

autodem2(epi_dir, [], [], [], [], [], "", "", "", [], "", "", "", "", epi_dem_dir, geocoded_dsm, "",
         "", [], [], dem_res, "")

#convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
dem_file = r"C:\geoeye\dsm_new\geoeye_dem.pix"
feature_size = [150]
gradient = [35]
remove_bumps = [9,5,9,5]
fill_pits = [7,5]
smooth_dem = [9]

dsm2dtm( geocoded_dsm, [1], dem_file, [1], [], [], "", feature_size, gradient, "hilly", remove_bumps, fill_pits, smooth_dem, [] )

#orthorectify images
ortho_dir = r"C:\geoeye\orthos"
if not os.path.exists(ortho_dir):
    os.makedirs(ortho_dir)

res = ["2", "2"]

ortho2( ingest_link_dir, [], [], [], "", ortho_dir, "TIF", "", [], "", "", "", "", [], "", "",
        res[0], res[1], dem_file, [], [], "", "", [], "", [], "" )

#Clean up unnecessary files and directories
del_dirs_list = [epi_dem_dir, epi_dir, work_dir, ingest_link_dir]
for del_dir in del_dirs_list:
    shutil.rmtree(del_dir)