**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.
    
    .. code-block:: python
        
        import os
        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( raw_airphoto, pix_link_file, [] )

        
Discover all airphotos in a directory and create .pix link file - **Batch Processing**
--------------------------------------------------------------------------------------
    This recipe demonstrates how to automatically discover all valid raw airphotos and then create a .pix link
    file for them so that they can be run through Geomatica's photogrammetric processing chain.

    .. code-block:: python
    
        import os, fnmatch, sys
        from pci.link import link


        raw_airphoto_dir =r"C:\aerial\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
        for raw_image in input_files_list:
            input_file = raw_image
            ingest_dir = r"C:\aerial\Ingested"
            link_pix = ingest_dir + os.sep + os.path.basename(os.path.splitext(raw_image)[0]) + ".pix"

            link(input_file, 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 formated camera 
    calibration xml file 
    
    :download:`Template Camera Calibration XML file <downloads/camera_calib.xml>`
    
    .. code-block:: python

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

        raw_airphoto_dir =r"C:\Aerial\photos\Subset"

        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
        for raw_image in input_files_list:
            input_file = raw_image
            ingest_dir = r"C:\Aerial\Ingested\Subset"
            link_pix = ingest_dir + os.sep + os.path.basename(os.path.splitext(raw_image)[0]) + ".pix"

            link(input_file, link_pix, [] )

        #Create OrthoEngine project file
        mfile = ingest_dir + os.sep + "*.pix"
        oeproj = ingest_dir + os.sep + "oe_initial.prj"
        mapunits = "UTM 17 D122"

        crproj( mfile, [], oeproj, "APDU", mapunits )

        #Import Camera Parameters
        camfile = raw_airphoto_dir + os.sep + "camera_calib.xml"
        prepared_oe_proj = ingest_dir + os.sep + "oe_configured.prj"

        camimport( oeproj, camfile, prepared_oe_proj )

        #Import Exterior Orientation File
        eofile = raw_airphoto_dir + os.sep + "eo.txt"
        mapunits = "UTM 17 D122" #D122 is NAD83

        eoimport( prepared_oe_proj, eofile, [3], mapunits, "", [], [], "" )
        
        
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)

    .. code-block:: python
        
        import os, sys, fnmatch
        from pci.cpmmseg import cpmmseg
        from pci.epipolar2 import epipolar2
        from pci.autodem2 import autodem2
        from pci.dsm2dtm import dsm2dtm

        #Create math model from OE to ingested images
        ingest_airphoto_dir = r"C:\Aerial\Ingested"
        prepared_oe_proj = ingest_airphoto_dir + os.sep + "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))

        for ingest_image in ingest_image_list:
            oeproj = final_oe_proj
            file = ingest_image

            cpmmseg(oeproj, file, [])

        #Generate Epipolar Images from Stereo pairs
        epi_outdir = ingest_airphoto_dir + os.sep + 'epipolars'
        if not os.path.exists(epi_outdir):
            os.makedirs(epi_outdir)

        epipolar2( ingest_airphoto_dir, [1,2,3], [], [], u"", u"opt", [],
                   [1], epi_outdir, [0], u"", u"", [1024] )

        #Extract Digital Surface Model (DSM) from multiple stereo airphoto images and mosaic them together
        dsm_outdir = r"C:\Aerial\Ingested"
        geocoded_dsm = r"C:\Aerial\Ingested\airphoto_dsm.pix"
        pxszout = [0.9, 0.9]

        autodem2( epi_outdir, [1], [], [], [], [], 'high', 'hilly', '32r', [2], 'yes',
                  'no', 'no', '', dsm_outdir, geocoded_dsm, '', 'UTM 17 D122', [], [], pxszout, 'blend' )

        #convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
        dem_file = r"C:\Aerial\Ingested\airphoto_dem.pix"
        feature_size = [250]
        gradient = [30]
        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, [] ) 


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)
    
    .. code-block:: python 
        
        from pci.cpmmseg import cpmmseg
        from pci.ortho2 import ortho2

        #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"

        cpmmseg(prepared_oe_proj, ingest_airphoto, [])

        #Orthorectify imagery with a user provided DEM
        ortho_airphoto = r"C:\Aerial\Orthos\A17563-001.pix"
        x_res = "0.2" #x resolution set to 20cm
        y_res = "0.2" #y resolution set to 20cm
        dem = r"C:\Aerial\DEM\1962_Saskatchewan_DEM.pix"

        ortho2(ingest_airphoto, [], [], [], "", ortho_airphoto, "PIX", "", [], "", "",
                "", "", [], "", "", x_res, y_res, 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)
    
    .. code-block:: python

        import os, sys, fnmatch
        from pci.cpmmseg import cpmmseg
        from pci.ortho2 import ortho2

        ingest_airphoto_dir = r"C:\test\ShawnPy\Data\Aerial\Ingested"
        prepared_oe_proj = ingest_airphoto_dir + os.sep + "oe_configured.prj"

        #Copy math model from prepared orthoengine project file to all ingested images
        for r,d,f in os.walk(ingest_airphoto_dir):
            for file in fnmatch.filter(f, "*.pix"):
                print "found valid input file: ", file
                ingest_image = os.path.join(r,file)

                cpmmseg(prepared_oe_proj, ingest_image, [])

        #Orthorectify imagery with a user provided DEM
        ortho_dir = r"C:\test\ShawnPy\Data\Aerial\Orthos"
        x_res = "0.2" #x resolution set to 20cm
        y_res = "0.2" #y resolution set to 20cm
        dem = r"C:\test\ShawnPy\Data\Aerial\dtm\PCI_Hamilton14_DTM_45cm.pix"

        ortho2(ingest_airphoto_dir, [], [], [], "", ortho_dir, "PIX", "", [], "", "",
                "", "", [], "", "", x_res, y_res, 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.

    .. code-block:: python

        import os, sys, fnmatch
        from pci.link import *
        from pci.crproj import *
        from pci.camimport import *
        from pci.eoimport import *
        from pci.cpmmseg import cpmmseg
        from pci.ortho2 import ortho2

        raw_airphoto_dir =r"C:\Aerial\photos\Subset"

        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
        for raw_image in input_files_list:
            input_file = raw_image
            ingest_dir = r"C:\Aerial\Ingested\Subset"
            link_pix = ingest_dir + os.sep + os.path.basename(os.path.splitext(raw_image)[0]) + ".pix"

            link( input_file, link_pix, [] )

        #Create OrthoEngine project file
        mfile = ingest_dir + os.sep + "*.pix"
        oeproj = ingest_dir + os.sep + "oe_initial.prj"
        mapunits = "UTM 17 D122"

        crproj( mfile, [], oeproj, "APDU", mapunits )

        #Import Camera Parameters
        camfile = raw_airphoto_dir + os.sep + "camera_calib.xml"
        prepared_oe_proj = ingest_dir + os.sep + "oe_configured.prj"

        camimport( oeproj, camfile, prepared_oe_proj )

        #Import Exterior Orientation File
        eofile = raw_airphoto_dir + os.sep + "eo.txt"
        mapunits = "UTM 17 D122" #D122 is NAD83

        eoimport( prepared_oe_proj, eofile, [3], mapunits, "", [], [], "" )

        #Copy math model from prepared orthoengine project file to all ingested images
        for r,d,f in os.walk(ingest_dir):
            for file in fnmatch.filter(f, "*.pix"):
                print "found valid input file: ", file
                ingest_image = os.path.join(r,file)

                cpmmseg(prepared_oe_proj, ingest_image, [])

        #Orthorectify imagery with a user provided DEM
        ortho_dir = r"C:\Aerial\Orthos\subset"
        x_res = "0.2" #x resolution set to 20cm
        y_res = "0.2" #y resolution set to 20cm
        dem = r"C:\Aerial\dtm\PCI_Hamilton14_DTM_45cm.pix"

        ortho2(ingest_dir, [], [], [], "", ortho_dir, "PIX", "", [], "", "",
                "", "", [], "", "", x_res, y_res, 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.

    .. code-block:: python

        import os, sys, fnmatch
        from pci.link import *
        from pci.crproj import *
        from pci.camimport import *
        from pci.eoimport import *
        from pci.cpmmseg import cpmmseg
        from pci.autodem2 import autodem2
        from pci.epipolar2 import epipolar2
        from pci.dsm2dtm import dsm2dtm
        from pci.ortho2 import ortho2

        raw_airphoto_dir =r"C:\Aerial\photos\Subset"

        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
        for raw_image in input_files_list:
            input_file = raw_image
            ingest_dir = r"C:\Aerial\Ingested\Subset"
            link_pix = ingest_dir + os.sep + os.path.basename(os.path.splitext(raw_image)[0]) + ".pix"

            link( input_file, link_pix, [] )

        #Create OrthoEngine project file
        mfile = ingest_dir + os.sep + "*.pix"
        oeproj = ingest_dir + os.sep + "oe_initial.prj"
        mapunits = "UTM 17 D122"

        crproj( mfile, [], oeproj, "APDU", mapunits )

        #Import Camera Parameters
        camfile = raw_airphoto_dir + os.sep + "camera_calib.xml"
        prepared_oe_proj = ingest_dir + os.sep + "oe_configured.prj"

        camimport( oeproj, camfile, prepared_oe_proj )

        #Import Exterior Orientation File
        eofile = raw_airphoto_dir + os.sep + "eo.txt"
        mapunits = "UTM 17 D122" #D122 is NAD83

        eoimport( prepared_oe_proj, eofile, [3], mapunits, "", [], [], "" )

        #Create math model from OE to ingested images
        ingest_image_list = []
        for r,d,f in os.walk(ingest_dir):
            for file in fnmatch.filter(f, "*.pix"):
                print "found valid input file: ", file
                ingest_image_list.append(os.path.join(r,file))

        for ingest_image in ingest_image_list:
            cpmmseg(prepared_oe_proj, ingest_image, [])

        #Generate Epipolar Images from Stereo pairs
        epi_outdir = ingest_dir  + os.sep + 'epipolars'

        epipolar2( ingest_dir, [1,2,3], [], [], u"", u"opt", [],
                   [1], epi_outdir, [0], u"", u"", [1024] )

        #Extract Digital Surface Model (DSM) from multiple stereo airphoto images and mosaic them together
        dsm_outdir = r"C:\Aerial\Ingested\Subset\DEM"
        geocoded_sm = r"C:\Aerial\Ingested\Subset\DEM\airphoto_dsm.pix"
        pxszout = [0.9, 0.9]

        autodem2( epi_outdir, [1], [], [], [], [], 'high', 'hilly', '32r', [2], 'yes',
                  'no', 'no', '', dsm_outdir, geocoded_sm, '', 'UTM 17 D122', [], [], pxszout, 'blend' )

        #convert DSM to DEM (Automatically remove surface features such as trees and buildings from DSM)
        dem_file = r"C:\Aerial\Ingested\Subset\DEM\airphoto_dem.pix"
        feature_size = [250]
        gradient = [30]
        remove_bumps = [9,5,9,5]
        fill_pits = [7,5]
        smooth_dem = [9]

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

        #Orthorectify imagery with a user provided DEM
        ortho_dir = r"C:\Aerial\Orthos\subset"
        x_res = "0.2" #x resolution set to 20cm
        y_res = "0.2" #y resolution set to 20cm

        ortho2(ingest_dir, [], [], [], "", ortho_dir, "PIX", "", [], "", "",
                "", "", [], "", "", x_res, y_res, 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.
    
    :download:`Template Camera Calibration XML file <downloads/camera_calib.xml>`
    
    .. code-block:: python
            
        import os

        #Set camera calibration location
        cam_dir = raw_input("Please provide the directory for the Camera Calibration file: ")
        cam_file = cam_dir + os.sep + "camera_calib.xml"
        target = open(cam_file, 'w+')

        #Get Camera Calibration Values
        x = 0
        while x == 0:
            cam_sensor = raw_input("What Camera Sensor are you using (UltraCam/DMC/Other)? ")
            
            if cam_sensor.lower() == "ultracam":
                print "Progressing to UltraCam setup mode..."
                x = 1
            elif cam_sensor.lower() == "dmc":
                print "Progressing to DMC setup mode..."
                x = 1
            elif cam_sensor.lower() == "other":
                print "Progressing to generic camera setup mode..."
                x = 1
            else:    
                print "Your answer could not be verified. Please respond with: ultracam, dmc or other"
                x = 0
                
        if cam_sensor.lower() == "ultracam":
            focal_length = raw_input("Please provide the focal length in mm: ")
            sym_ppx =  raw_input("Please provide the principal Point offset x: ")
            sym_ppy =  raw_input("Please provide the principal Point offset y: ")
            chip_width = raw_input("Please provide the chip size width in mm: ")
            chip_length = raw_input("Please provide the chip size length in mm: ")
        elif cam_sensor.lower() == "dmc":
            focal_length = raw_input("Please provide the focal length in mm: ")
            chip_width = raw_input("Please provide the chip size width in mm: ")
            chip_length = raw_input("Please provide the chip size length in mm: ")
        else:
            focal_length = raw_input("Please provide the focal length in mm: ")
            sym_ppx =  raw_input("Please provide the principal Point offset x: ")
            sym_ppy =  raw_input("Please provide the principal Point offset y: ")
            chip_width = raw_input("Please provide the chip size width in mm: ")
            chip_length = raw_input("Please provide the chip size length in mm: ")
            
            radial_distortion = raw_input("Do you have radial distortion values? (y/n)")
            if radial_distortion == "y":
                rldc_r0 = raw_input("please provide RLDC R0 value: ")
                rldc_r1 = raw_input("please provide RLDC R1 value: ")
                rldc_r2 = raw_input("please provide RLDC R2 value: ")
                rldc_r3 = raw_input("please provide RLDC R3 value: ")
                rldc_r4 = raw_input("please provide RLDC R4 value: ")
                rldc_r5 = raw_input("please provide RLDC R5 value: ")
                rldc_r6 = raw_input("please provide RLDC R6 value: ")
                rldc_r7 = raw_input("please provide RLDC R7 value: ")
            else:
                print "Setting radial coefficients to 0 (zero)"
                rldc_r0 = 0
                rldc_r1 = 0
                rldc_r2 = 0
                rldc_r3 = 0
                rldc_r4 = 0
                rldc_r5 = 0
                rldc_r6 = 0
                rldc_r7 = 0    
                
            decentering_distortion = raw_input("Do you have decentering distortion values? (y/n)")
            if decentering_distortion == "y":
                ddc_p1 = raw_input("please provide DDC P1 value: ")
                ddc_p2 = raw_input("please provide DDC P2 value: ")
                ddc_p3 = raw_input("please provide DDC P3 value: ")
                ddc_p4 = raw_input("please provide DDC P4 value: ")
            else:
                print "Setting decentering coefficients to 0 (zero)"
                ddc_p1 = 0
                ddc_p2 = 0
                ddc_p3 = 0
                ddc_p4 = 0
                
        #Create Camera Calibration XML file

        target.write("<CAMERA NAME=\"" + cam_sensor.upper() + "\">\n")
        target.write("<LENS_SERIAL/>\n")
        target.write("<CALIB_SOURCE/>\n")
        target.write("<CALIB_DATE/>\n")
        target.write("<FL>" + str(focal_length) + "</FL>\n")
        target.write("<SYMMETRY PPX=\"" + str(sym_ppx) + "\" PPY=\"" + str(sym_ppy) + "\"/>\n")
        target.write("<RLDC name=\"R0\">" + str(rldc_r0) + "</RLDC>\n")
        target.write("<RLDC name=\"R1\">" + str(rldc_r1) + "</RLDC>\n")
        target.write("<RLDC name=\"R2\">" + str(rldc_r2) + "</RLDC>\n")
        target.write("<RLDC name=\"R3\">" + str(rldc_r3) + "</RLDC>\n")
        target.write("<RLDC name=\"R4\">" + str(rldc_r4) + "</RLDC>\n")
        target.write("<RLDC name=\"R5\">" + str(rldc_r5) + "</RLDC>\n")
        target.write("<RLDC name=\"R6\">" + str(rldc_r6) + "</RLDC>\n")
        target.write("<RLDC name=\"R7\">" + str(rldc_r7) + "</RLDC>\n")
        target.write("<DDC name=\"P1\">" + str(ddc_p1) + "</DDC>\n")
        target.write("<DDC name=\"P2\">" + str(ddc_p2) + "</DDC>\n")
        target.write("<DDC name=\"P3\">" + str(ddc_p3) + "</DDC>\n")
        target.write("<DDC name=\"P4\">" + str(ddc_p4) + "</DDC>\n")
        target.write("<DIGITAL>\n")
        target.write("<CHIP_SIZE WIDTH=\"" + str(chip_width) + "\" HEIGHT=\"" + str(chip_length) + "\"/>\n")
        target.write("<YSCALE>1.000</YSCALE>\n")
        target.write("</DIGITAL>\n")
        target.write("</CAMERA>")

        print ""
        print "completed writing Camera Calibration file"
        print ""
        print "Camera calibration file located here: " + cam_file