**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 .. code-block:: python 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. .. code-block:: python 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" Create .pix Link File for Satellite Images from New Sensors ----------------------------------------------------------- This recipe demonstrates how to create a .pix link file from an input image (often raw satellite image from a newer sensor). A link file creates a .pix file, but rather than copying the imagery, it only links back to the original image file. However, the file will have all the same file support as a regular .pix file. Importing the file (fimport) creates a copy of the input image in .pix format, which uses more disk space and takes more time to process. The link file is a nice alternative to quickly import data into the .pix format while saving on diskspace and processing time. It is important to keep in mind that the relative path between the link file and the raw imagery must be maintained. Otherwise the link will be broken. .. code-block:: python from pci.link import link raw_ms_geoeye = r"C:\geoeye\po_311159_0000000\po_311159_bgrn_0000000.tif" link_image = r"C:\geoeye\Output\link_GeoEye.pix" link( raw_ms_geoeye, link_image, [] ) 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. .. code-block:: python 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. .. code-block:: python 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. .. code-block:: python 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. .. code-block:: python #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 .. code-block:: python 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 .. code-block:: python #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. .. code-block:: python #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 .. code-block:: python 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 .. code-block:: python 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) .. code-block:: python 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. .. code-block:: python 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. .. code-block:: python 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)