This tutorial walks you through the steps necessary to improve the accuracy of the RPCs (Rational Polynomial Coefficients) provided with satellite imagery, by automatically collecting GCPs (Ground Control Points) and tie points. The new and improved RPCs will then be exported to a text file in a common ASCII format (format used by DigitalGlobe).
Being able to improve the quality and accuracy of your RPCs by collecting high quality GCPs and tie points and then exporting it to a common form1at is beneficial for:
- Using an improved geometric model in 3rd party software. For example, 3D Stereo software for DEM Editing or 3D Feature Extraction
- Distributing improved geometric models in common formats to clients
- Geomatica 2016 or later installed
- Python 101 (basic understanding of python coding concepts)
- How to access modules, classes, attributes and methods
- How to call a function
- How to use with statements
- How to create if-statements
- Some experience programming with Geomatica (recommended)
- Click Here to download the data required for this tutorial
To demonstrate this workflow, we will use Pleiades tri-stereo imagery over captured over Melboune, Australia. Each of the three tri-stereo images are distributed with Rational Polynomial Coefficients (RPCs), which are held in the RPC*.XML file for each scene. Through this tutorial, you will automatically collect sub-pixel GCPs, sub-pixel tie points and use them to improve the geometric model of the images. We will then calculate new RPCs based on the improved model and write it to a common ASCII (text) file.
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 file(s) and output file(s).
- link – Geomatica’s algorithm for creating a lightweight link file from raw satellite formats. This is required for photogrammetric operations.
- autogcp – Geomatica’s algorithm for automatically finding ground control points (GCPs) based on image matching techniques
- gcprefn – Geomatica’s algorithm for automatically detecting bad GCP matches and removing them from the model
- crproj – Geomatica’s algorithm for creating OrthoEngine projects, which are required for performing a block bundle adjustment (tie points)
- autotie – Geomatica’s algorithm for automatically collecting tie points between neighboring images in a project
- tprefn – Geomatica’s algorithm for automatically detecting bad tie points and removing them from the model
- oemodel – Geomatica’s algorithm for computing a block bundle adjustment and writing the updated geometric model to the image .pix files
- model2rpc – Geomatica’s algorithm for calculating new RPCs from the updated geometric models in the .pix files and writing the coefficients to an ASCII (text) file.
- Report, enableDefaultReport – are two functions/classes that allow dumping output notifications to a report file
- PCIException – is a class for handling errors specific to Geomatica’s algorithms
On line 14 we create a variable called working_dir, set it equal to the path where you extracted the tutorial data (root folder). On lines 15 to 21, we create a number of variables hold paths to various input files and folders and output files and folders, which are required by this workflow.
2. Check if output directories exist. If not, create them
We run a simple for-loop that checks to see if the ingest_gcp_dir and bundle_dir directories exist. If not, they are created.
On line 24, we create a for-loop that first checks the ingest_gcp_dir and then checks the bundle_dir. Line 25 is a condition with if not statement that checks to see if each directory does not exist. If it does exist, it will skip to the next directory to check. If it does not exist, it will run line 26 and make the directory.
Note: os.makedirs() will make all directories in the path that do not exist.
3. Search input directory for all valid images
In this step, we will create a simple generator that will search the input directory for all valid input files and prepare them for processing
On line 30 we define a function called get_batch, which takes in the in_dir variable and img_filter variable. The in_dir variable points to the input directory that it should recursively check for valid input files in. The img_filter variable is a string that is used to determine which files are valid. Lines 31 to 32 loop through the directory and sub directories (if exist) and search for all files that match the img_filter string, using fnmatch.
On line 33, the yield operator is used to return every item (valid input file) on-the-fly.
On line 35, we output each iteration from the get_batch generator and temporarily store it in the raw_image_batch variable.
4. Loop through raw images and create link .pix file for them
In this step, we will create a for-loop that will loop through each input image one at a time and perform all processing (link à GCP collection à GCP Refinement). However, at this point in the tutorial, we will only focus on the link part of the loop.
On line 37 we create the for-loop that all ingesting and GCP collection code will reside in. In the for-loop we will get each raw_img item that is yielded by raw_image_batch.
On line 38, we create a variable called base_name and set it equal to the filename of the input file. The os.path.splitext() function is used to retrieve the filename without the path or file extension. Line 39 defines a variable called pix_image, and use the function os.path.join() to add the output path to the base_name and append a suffix ‘_ingested’ and the .pix file extension.
On line 40 we run Geomatica algorithm, link() and pass the following input keyword arguments:
- fili=raw_img – raw_img (input image) variable is assigned to the fili keyword
- filo=pix_img – pix_img (output image) variable is assigned to the filo keyword
5. Automatically collect and refine GCPs
The remaining code in the for-loop is used to setup and run the Geomatica algorithms autogcp() and gcprefn() in order to collect GCPs and remove potential blunders.
On line 42 we define a variable called refined_gcp_seg that will hold an empty list. It will be used as an output parameter by the gcprefn() algorithm and populated with a value that tells us which .pix segment the refined GCPs are held in.
On line 47 we create a try: statement, which tells python to first try running the code inside and if an error occurs, handle it afterwards. This is known as EAFP (Easier to ask forgiveness than permission).
On lines 44 to 48 we call the autogcp() algorithm. We set it equal to gcp_seg, as the algorithm returns the output segment that holds the GCPs that were collected. Check out the autogcp() documentation to understand what each parameter does.
On lines 49 to 54 we call the gcprefn() algorithm. In the set of arguments passed when calling the gcprefn() algorithm, we set lasc=refined_gcp_seg variable, which will be populated with the .pix segment number that contains the refined GCPs. Check out the gcprefn() documentation to understand what each parameter does.
Finally, on lines 55 and 56, we handle any PCI exception that occurs by printing the exception and moving on.
On line 57 we exit the for-loop by tabbing back
6. Automatically collect and refine tie points and compute bundle adjustment
After exiting the for-loop used for importing and GCP collection, we run the necessary Geomatica algorithms for collection tie points and performing a block bundle adjustment.
On line 59 we create a try: statement, which tells python to first try running the code inside and if an error occurs, handle it afterwards. This is known as EAFP (Easier to ask forgiveness than permission).
On lines 60 to 65, we create an OrthoEngine project file and import the refined GCPs (refined_gcp_seg), by calling the Geomatica algorithm crproj(). Check out the crproj() documentation to understand what each parameter does.
Next we run the Geomatica algorithm autotie() and tprefn() on lines 66 to 70 and 71, respectively. The autotie() algorithm automatically collects tie points between adjacent images and tprefn() removes potential tie point blunders. On line 72, we call the Geomatica algorithm oemodel(), which performs the block bundle adjustment based on the GCPs and tie points in the project.
Finally, on lines 73 and 74, we handle any PCI exception that occurs by printing the exception and moving on.
7. Calculate new RPCs based on improved model and export to ASCII file
We call the get_batch() function we created near the top of this script to get all of the link files with the updated geometric models. Then calculate the RPCs from these updated models and write them to a new text file.
On line 77, we output each iteration from the get_batch() generator and temporarily store it in the updated_rpc_batch variable. Note that in this call, we are using different arguments (pix_img and ‘DIM*.pix’)
On line 79 we create the for-loop that loops through each of the three images that we updated the geometric models for and then calculates and exports the RPCs for each image using the Geomatica algorithm model2rpc().
On line 80, we create a variable called base_name and set it equal to the filename of the input file. The os.path.splitext() function is used to retrieve the filename without the path or file extension. Line 81 defines a variable called rpc_file, and use the function os.path.join() to add the output path to the base_name and append a suffix ‘_RPC’ and the .txt file extension.
Note: The output DIM*_RPC.txt files will be written to the output folder (C:\...\Improve_export_rpc\output\)