Satellite Photogrammetry¶

from pci.fimport import fimport

fimport(fili=r"C:\Kompsat\MSC_080407140610_09045_24220801PMSN14_1G.eph",
        filo=r"C:\Kompsat\Ingest\kompsat2_image.pix", dblayout="TILED512")
print "Complete"

from pci.link import link

link(fili=r"C:\geoeye\po_311159_0000000\po_311159_bgrn_0000000.tif", filo=r"C:\geoeye\Output\link_GeoEye.pix")

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(fili=raw_ms_l8, filo=pix_ms_l8)

# 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(fili=raw_pan_l8, filo=pix_pan_l8)

# 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(fili=raw_therm_l8, filo=pix_therm_l8)

from pci.link import link
from pci.autogcp import autogcp
from pci.gcprefn import gcprefn

# 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(fili=raw_geoeye_image, filo=link_pix_file)

# Run GCP Collection and store new GCP segment in link pix file
reference_image = r"c:\geoeye\Ref\S-55-40_lr_2000.tif"
dem = r"C:\geoeye\po_311159_0000000\dem\gmted2010.pix"
gcp_seg = autogcp(mfile=link_pix_file, dbic=[3],filref=reference_image,
                  dbic_ref=[3], filedem=dem, searchr=[50])

# autogcp function returns the new gcp segment (saved to variable gcp_seg), which can then be directly passed to the gcp refinement function
# Setup required parameters for automatic GCP refinement
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(file=link_pix_file, dbgc=gcp_seg, modeltyp="RF", order=order, reject=reject)
# the refined GCPs are now stored in the link pix file and can be used for further photogrammetric processing

import fnmatch, os, sys
from pci.link import link
from pci.autogcp import autogcp
from pci.gcprefn import gcprefn

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 and then adds .pix extension
        ingest_link_file = os.path.join(ingest_link_dir, ingest_link_filename)
        link(fili=raw_file, filo=ingest_link_file)

        # Perform automatic GCP Collection
        ref = r"c:\geoeye\Ref\S-55-40_lr_2000.tif"
        dem = r"C:\geoeye\dem\gmted2010.pix"
        gcp_seg = autogcp(mfile=ingest_link_file, dbic=[1],filref=ref,
                  dbic_ref=[3], filedem=dem, searchr=[50])

        # Refine the GCPs based on RMSE
        gcprefn(file=ingest_link_file, dbgc=gcp_seg, modeltyp="RF", order=[0], reject=[5,0.5,0.5])

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

import fnmatch, os, sys
from pci.autogcp 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 = os.path.join(r"D:\Kompsat-2\ingested", os.path.splitext(os.path.basename(raw_file))[0] + ".pix")
        link(fili=raw_file, filo=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

            # Run autogcp to collect GCPs
            gcp_match_chan = [1]
            ref = r"D:\Kompsat-2\Reference\SPOT\o2_5m.pix"
            dem = 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 = autogcp(mfile=ingest_file, dbic=[1],filref=ref,
                              filedem=dem, searchr=search_radius, numgcps=numgcps)
            increase_radius += 50 # adds 50 pixels to the search radius for the next iteration (If required)

        # Run automatic GCP refinement (remove bad gcps)
        gcprefn(file=ingest_file, dbgc=[-1], modeltyp="RF", order=[0], reject=[5,1,1])

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

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(fili=raw_kompsat2_file, filo=imported_k2_image, dblayout="TILED512")

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

imageinv(fili=imported_k2_image, filo=inverse_image, dbic=[3])

# copy the inverted image channel into the imported image that will be used for GCP collection
pcimod(file=imported_k2_image, pciop="ADD", pcival=[0,0,0,1])

iii(fili=inverse_image, filo=imported_k2_image, dbic=[1], dboc=[5])

# Run FFTMVEC to automatically collect GCPs from road vectors
roadvec = r"C:\Kompsat-2\Reference\roads\toronto_vector.pix"
gcp_seg = [] # the fftmvec lasc parameter will populate this variable with the output gcp segment number

fftmvec(file=imported_k2_image, dbic=[1,2,4], dbinvc=[5], filv=roadvec, dbvs=[2], vecwidth=[15], lasc=gcp_seg)

# Refine GCPs - remove bad GCPs until RMSE falls below 1 pixel in x and y direction

gcprefn(file=imported_k2_image, dbgc=gcp_seg, modeltyp="RF", order=[0], reject=[5,1,1])

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

from pci.rfmodel import rfmodel
from pci.ortho import ortho

# Calculate the mathematical model from GCPs (rational function model) required to orthorectify the imagery
ingest_dir = r"D:\Python_Cookbook\Kompsat-2\ingested\*.pix"

rfmodel(mfile=ingest_dir, dbgc=[9], numcoeff=[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"

ortho(mfile=ingest_dir, filo=ortho_dir, ftype="tif", bxpxsz="4", bypxsz="4", filedem=dem)

import os
from pci.autotie import autotie
from pci.tprefn import tprefn
from pci.crproj import crproj
from pci.cpmmseg import cpmmseg
from pci.ortho import ortho

# Create OrthoEngine Project, which is required for the bundle adjustment
image_dir = r"C:\kompsat\ingest"
work_dir = r"C:\kompsat\working\TiePoints"
proj = os.path.join(work_dir, "import_oemodel.prj")
projection = "UTM 17 T D000" # D000 is WGS84

crproj(mfile=image_dir, dbgc=[4], oeproj=proj, modeltyp=u"RFI1", mapunits=projection)

# Run Automatic Tie Point Collection
tp_oeproj = os.path.join(work_dir,"output_oetp.prj")
match_chan = [4] # NIR channel often yields best matching results
dem = r"C:\dem\gmted2010.pix"

autotie(mfile=tp_oeproj, oeprojo=tp_oeproj, dbic=match_chan, filedem=dem)

# Run automatic Tie Point Refinement

tprefn(oeproji=tp_oeproj, reject=[5, 1.0, 1.0])

# Copy Math Model out to file
# The second parameter in cpmmseg is FILE, which provides the name of the input image file for which the math model will be copied.
# If FILE is not specified, all images contained within the project are copied.
cpmmseg(oeproj=tp_oeproj)

# Orthorectify images
ortho_dir = r"C:\Kompsat\orthos"
format = "TIF"

ortho(mfile=image_dir, filo=ortho_dir, ftype=format, bxpxsz=[4], bypxsz=[4], filedem=dem)

from pci.rfmodel import rfmodel
from pci.epipolar import epipolar
from pci.autodem import autodem
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(mfile=ingest_dir, dbgc=gcp_seg, numcoeff=[0])

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

epipolar(mfile=ingest_dir, dbic=gcp_seg, selmthd="max", sampling=[1], outdir=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]

autodem(mfile=epi_dir, outdir=out_dir, filedem=geocoded_dsm, pxszout=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(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)

import fnmatch, os, sys
from pci.link import *
from pci.autogcp import *
from pci.gcprefn import *
from pci.autotie import *
from pci.tprefn import *
from pci.crproj import *
from pci.cpmmseg import *
from pci.ortho import *

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

# Create the link file directory if it does not yet exist
if not os.path.exists(ingest_link_dir):
    os.makedirs(ingest_link_dir)

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_filename = os.path.splitext(os.path.basename(raw_file))[0] + ".pix" # get the filename of the raw image (without extension)
        ingest_link_file = os.path.join(ingest_link_dir, ingest_link_filename)

        link(fili=raw_file, filo=ingest_link_file)

        # Collect GCPs with autogcp
        ref_image = r"C:\geoeye\Ref\S-55-40_lr_2000.tif"

        gcp_seg = autogcp(mfile=ingest_link_file, dbic=[1], filref=ref_image,filedem=dem, searchr=[50])

        # Refine GCPs with gcprefn
        gcprefn(file=ingest_link_file, dbgc=gcp_seg, modeltyp="RF", order=[0], reject=[5,2.0,2.0])

else:
    print "No image files found in " + raw_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)

proj = os.path.join(work_dir, "import_oemodel.prj")
projection = "UTM 55 G D000" # D000 is WGS84
refined_gcps = [-1] # setting the dbgc parameter to -1 will ensure that the last GCP segment in a file will be imported

crproj(mfile=ingest_link_dir, dbgc=[-1], oeproj=proj, modeltyp="RFI1", mapunits=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

autotie(mfile=proj, oeprojo=tp_oeproj, dbic=tp_match_chan, filedem=dem)

# Run automatic Tie Point Refinement
tprefn(oeproji=tp_oeproj, reject=[5, 0.5, 0.5])

# Copy Math Model out to file
# The second parameter in cpmmseg is FILE, which provides the name of the input image file for which the math model will be copied.
# If FILE is not specified, all images contained within the project are copied.
cpmmseg(oeproj=tp_oeproj)

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

ortho(mfile=ingest_link_dir, filo=ortho_dir, ftype="TIF", bxpxsz=[2], bypxsz=[2], filedem=dem)

import fnmatch, os, sys, shutil
from pci.link import *
from pci.autogcp import *
from pci.gcprefn import *
from pci.autotie import *
from pci.tprefn import *
from pci.crproj import *
from pci.cpmmseg import *
from pci.epipolar import *
from pci.autodem import *
from pci.dsm2dtm import *
from pci.ortho import *
from pci.exceptions import *

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

# Create the link file directory if it does not yet exist
if not os.path.exists(ingest_link_dir):
    os.makedirs(ingest_link_dir)

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_filename = os.path.splitext(os.path.basename(raw_file))[0] + ".pix" # get the filename of the raw image (without extension)
        ingest_link_file = os.path.join(ingest_link_dir, ingest_link_filename)

        link(fili=raw_file, filo=ingest_link_file)

        # Collect GCPs with autogcp
        ref_image = r"C:\geoeye\Ref\S-55-40_lr_2000.tif"

        gcp_seg = autogcp(mfile=ingest_link_file, dbic=[1], filref=ref_image,filedem=dem, searchr=[50])

        # Refine GCPs with gcprefn
        gcprefn(file=ingest_link_file, dbgc=gcp_seg, modeltyp="RF", order=[0], reject=[5,2.0,2.0])

else:
    print "No image files found in " + raw_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)

proj = os.path.join(work_dir, "import_oemodel.prj")
projection = "UTM 55 G D000" # D000 is WGS84
refined_gcps = [-1] # setting the dbgc parameter to -1 will ensure that the last GCP segment in a file will be imported

crproj(mfile=ingest_link_dir, dbgc=[-1], oeproj=proj, modeltyp="RFI1", mapunits=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

autotie(mfile=proj, oeprojo=tp_oeproj, dbic=tp_match_chan, filedem=dem)

# Run automatic Tie Point Refinement
tprefn(oeproji=tp_oeproj, reject=[5, 0.5, 0.5])

# Copy Math Model out to file
# The second parameter in cpmmseg is FILE, which provides the name of the input image file for which the math model will be copied.
# If FILE is not specified, all images contained within the project are copied.
cpmmseg(oeproj=tp_oeproj)

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

epipolar(mfile=ingest_link_dir, dbic=[3], selmthd="max", sampling=[1], outdir=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))

autodem(mfile=epi_dir, outdir=epi_dem_dir, filedem=geocoded_dsm, pxszout=[6,6])

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

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 images
ortho_dir = r"C:\geoeye\orthos"
if not os.path.exists(ortho_dir):
    os.makedirs(ortho_dir)

ortho(mfile=ingest_link_dir, filo=ortho_dir, ftype="TIF", bxpxsz=[2], bypxsz=[2], filedem=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)