Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

Building Your First Geomatica Python Workflow (Orthorectification)

PCI Geomatics -

Overview


 

This tutorial outlines how to run multiple Geomatica algorithms in python to complete a workflow. We will work through an orthorectification workflow using Kompsat 3 imagery.

Benefits


 

Automating a workflow, such as orthorectification can save you time. Instead of individually running an algorithm, waiting and then running another algorithm, all of the algorithms will run in sequential order.

Prerequisites


 

  • Geomatica 2015 or later installed
  • Python 101 (basic understanding of python coding concepts)
    • try, except, finally statements
    • How to access modules and classes
    • How to call a function
  • Some experience programming with Geomatica (recommended)

Data Package


 

  • Click Here to download the data required for this tutorial

Tutorial


 

To demonstrate how to create a script which will automate an entire Geomatica workflow, we will orthorectify a Kompsat 3 image. The workflow that we will automate is as follows:

Import Kompsat 3 imagery
(fimport)

Collect and refine GCPs (autogcp & gcprefn)

Calculate model and orthorectify (satmodel & ortho)

 

1. Import necessary modules & setup inputs/outputs


 

The first step is to import the Geomatica modules we will be calling in this tutorial and setting up the variables that will point to our input and output directories and files.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
from pci.fimport import fimport
from pci.autogcp import autogcp
from pci.gcprefn import gcprefn
from pci.satmodel import satmodel
from pci.ortho import ortho
import os

# Set directories and pathnames
working_dir = r'I:\Data\Tutorial\Python\OrthoWorkflow'
image_header = os.path.join(working_dir, 'raw', 'K3_20131003174133_07365_20171339_L1R_Aux.xml-MS')
pix_image = os.path.join(working_dir, 'K3_Import.pix')
ref_image = os.path.join(working_dir, 'ref', 'LC80160282015300LGN00_CLIP.pix')
dem = os.path.join(working_dir, 'ref', 'Ottawa-NCR_LCC_DEM_10m.pix')
output_dir = os.path.join(working_dir, 'Ortho')
if not os.path.isdir(output_dir):
    os.mkdir(output_dir)

# Import Landsat 7 image
fimport(fili=image_header,
        filo=pix_image)

 

Lines 1 to 5 import the five required algorithms, as outlined above. The os module, which we will use to access the operating system’s file structure, is imported on line 6. On line 9, working_dir is a variable that contains the path to the directory where your tutorial data is stored. You will need to change this to the directory on your computer where your tutorial data is saved. On lines 10 to 13 each of the variables contain a required file. The image_header is the raw Kompsat 3 xml file, while pix_image will be the output image from the ingest step. Both ref_image and dem will be used as inputs for autogcp.

The output directory is then defined on line 14. In this variable definition, os.path.join() is used to create a new pathname string that combines working_dir with a new folder name. On line 15 os.path.isdir() is used to check if output_dir is currently a directory. If it is a directory, line 16 will be skipped. If it is not a directory, line 16 will be executed. The os.mkdir() function will create the new output_dir directory.

2. Import Kompsat 3 image


 

The first step in the workflow is to import the Kompsat 3 image to create a new pix file. This pix file will be used in any further algorithms.

18
19
20
# Import Landsat 7 image
fimport(fili=image_header,
        filo=pix_image)

 

In the code above fimport is used to import the raw imagery (image_header) into a new pix file (pix_image). All of the Kompsat 3 bands will be imported including the georeferencing and orbit segment.

3. Collect and Refine GCPs


 

In order to orthorectify imagery, GCPs must be collected. Autogcp automatically collects GCPs for your imagery and then returns a GCP segment number. You can then use this segment number in gcprefn to refine the newly created GCPs.

 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
#  GCP Collection and Refinement
autogcp_seg = autogcp(mfile=pix_image,
                      dbic=[1, 2, 3, 4],
                      filref=ref_image,
                      filedem=dem,
                      searchr=[50])

# Refine the GCPs created by autogcp2

refngcp_seg = []

gcprefn(file=pix_image,
        dbgc=autogcp_seg,
        orbit=[3],
        modeltyp=u'SAT',
        reject=[5, 0.5, 0.5],
        lasc=refngcp_seg)

 

First the GCPs need to be collected. However, we need to make sure that we know which new segment is created by autogcp. Therefore we need make the call to autogcp a variable (autogcp_seg). We can then use autogcp_seg in gcprefn. On line 23 the variable autogcp_seg is set and we call the algorithm autogcp. On lines 23-26 the input parameters are set. This is when the two variables that you set in step 1, ref_image and dem, are required. In this example many of the variables are left as default. We only set the search radius (searchr). You can check the help for an explanation of all the input parameters.

Next, we will refine the GCPs that were created in autogcp. The gcprefn algorithm uses a computer model to remove ground control points that have high errors. In our case we will use Toutin’s Satellite Orbital mode (SAT). The newly edited GCPs are then saved to a new GCP segment. You will need to save the segment number of this new GCP segment to a list. Create the variable refngcpSeg and define it as an empty list []. The output segment autogcp_seg will be used as the input segment on line 34. Gcprefn will then refine the GCPs from autogcp_seg, removing GCPs with large errors. The final GCPs are then saved to a new segment using the lasc variable. The empty list (refngcp_seg) that you specified earlier will be populated with the new segment number.

In this example two other parameters are listed reject. The reject parameter specifies the method used for rejecting (refining) the GCPs. There are five different rejection modes and in this example mode 5, RMS Error, will be used. GCP points with the highest residuals will be removed until xRMS and yRMS are both 0.5.

4. Calculate Model and Orthorectify


 

The GCPs are now ready to be used in model calculation and orthorectification.

 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
print refngcp_seg

#  Calculate Rigorous Model and Orthorectify
satmodel(mfile=pix_image,
         dbgc=refngcp_seg,
         orbit=[3])

# Create an orthorectified image
ortho(mfile=pix_image,
      filo=os.path.join(output_dir, 'output_ortho.pix'),
      filedem=dem)

 

On lines 44 to 46 the satellite model is calculated using satmodel. This algorithm computes the math model of the image using Toutin's model and stores the result as a binary segment in the image file.

Now that GCPs have been collected and the model has been calculated we can orthorectify our Kompsat 3 image. On lines 49 and 50 the input and output images are specified. The filedem parameter is also set to the variable dem that we defined originally. In this example all of the other ortho parameters are left as their defaults. When you are running all of these algorithms on your own datasets you may have to adjust the parameter values to better represent your specific dataset. A list of parameters for each of these algorithms is available in the Geomatica Python help.

Have more questions? Submit a request

Comments

Powered by Zendesk