The Geomatica InSAR package can be used to generate topographic products to characterize digital surface models (DSMs) or deformation products which identify subtle changes in surface elevation due to land subsidence. This InSAR tutorial is intended to provide two simple workflows to allow users to extract topographic or deformation products from a single interferometric pair.
This is an introductory tutorial to provide information on running InSAR algorithms and a simple workflow. For the full deformation workflow, multiple pairs would be processed and the deformation products would be stacked in chronological order for temporal analysis using the INSSTACK algorithm. For more information on INSTACK see: Deformation: Further Processing - ORTHO, INSSTACK
A full end-to-end workflow document is available from: InSAR Processing - Deformation and Topographic Workflows. This link includes a series of python scripts and a user manual which outline the full InSAR workflow when working with a stack of images. The manual outlines both the deformation and topographic workflows and includes additional information on running each of the Python scripts.
The InSAR User Manual provides a much more detailed explanation of the InSAR tools available, parameter options, helpful hints and troubleshooting steps for InSAR processing.
The InSAR workflow must be run in either Python or EASI. In this tutorial, these algorithms are run in Python. More information on creating python scripts using Geomatica libraries is available from the Developer Zone. If you are just starting out in Python, the Getting Started page is quite helpful. There are also a number of Geomatica Developer tutorials.
All of Geomatica algorithm parameters are outlined in the Geomatica Help. You can check the help for detailed information on each parameter.
This tutorial will outline the following steps of the InSAR workflow.
|INSTOPO – No External DEM||INSTOPO – External DEM|
User Manual – Pg 10
Coarse DEM (for optional adjustment)
Reference and Dependent Datasets: Prior to selecting your data sets, it is important to know what type of interferometric product you wish to produce. In general, it is a good idea to use data sets with the smallest time difference between acquisitions. This will minimize the effects of temporal decorrelation. If possible, you should avoid acquisitions that occur during or shortly after significant weather events such as rain or snow falls.
If you plan to generate a topographic product such as a DEM, you should select data pairs with a perpendicular baseline in the range of hundreds of meters. If you plan to generate subsidence (deformation) products you should select data sets with the smallest possible perpendicular baseline.
DEM: If the ELEVATION_DATUM metadata tag is set, the appropriate elevation model is used. If the elevation model is not specified, the DEM is assumed to be in meters above mean sea level (MSL). The accuracy of the topographic phase removal step is directly related to the accuracy of the elevations extracted from the DEM. The vertical accuracy of the DEM should be better than ± 20m. The spacing of the DEM should be within a factor of 10 spacing of the SAR data.
For this tutorial two Radarsat-2 datasets from Phoenix, Arizona are used. The datasets, DEM file, INSINFO report and full python scripts for each workflow (topographic/deformation) can be downloaded from the following link – https://pcigeomatics.sharefile.com/d-s4fb7ae7e8fc4cd3b
Note that in the data download there are five scripts:
InSAR_Topo_Sentinel1TOPS.py, InSAR_Defo_Sentinel1TOPS.py, InSAR_Topo.py, InSAR_Defo.py and INSINFO.py.
In this tutorial the scripts InSAR_Topo.py and InSAR_Defo.py are outlined. The “Sentinel1TOPS” scripts specifically process Sentinel-1 TOPS data. The only change between the TOPS scripts and the other scripts is the create_proc_window helper function.
* Note that currently Sentinel-1 TOPS data may not produce good results with the topographic workflow as the perpendicular baseline between Sentinel-1 TOPS pairs is not high enough (in the hundreds of meter range).
The INSINFO script can be run on any type of InSAR data. As outlined below, this algorithm is run separately from the other InSAR algorithms.
Generate Interferometric Parameter Report - INSINFO
User Manual – Pg. 21
This algorithm is not included in the Python script but can be run separately (in Python/EASI/Algorithm Librarian). This module can be used to report the interferometric parameters for SAR data sets. The INSINFO module can read data either directly from the vendor supplied format or from the PCIDSK format after data ingest (see SARINGEST). If raw vendor data is specified as input; the viewing geometry parameters such as slant ranges, incident angles and baseline information (first; middle and last) refer to the entire image dataset. If the output from SARINGEST is used as input, the slant range, incident angle range and baseline information correspond to the image subset (set by dbiw).
As noted above, when choosing reference and dependent datasets you want to note the time difference and perpendicular baseline information. This information is available in the INSINFO report. The report for the datasets used in this tutorial is available from the data download.
from pci import algo
full_ref_insar_image = r"D:\Data\Sensors\RADARSAT2\Insar\2008_05_04\Phoenix_RS2_UltraFine12_20080504_HH_SLC\product.xml"
full_dep_insar_image = r"D:\Data\Sensors\RADARSAT2\Insar\2008_05_28\Phoenix_RS2_UltraFine12_20080528_HH_SLC\product.xml"
info_report = r"D:\Data\Sensors\RADARSAT2\Insar\RS2_insinfo.txt"
algo.insinfo(filref=full_ref_insar_image, mfile=full_dep_insar_image, tfile=info_report)
Import modules and setup file variables
You will first need to import the Geomatica library, Geomatica Python API and Python OS module.
The file and directory variables then need to be set. The user also needs to provide upper left and lower right geocoded coordinates for a subset window (ul_coord, lr_coord). You can open the SAR image in Focus to determine the coordinates for your subset area.
A check has been added to ensure that the output directory (outdir) exists and if it does not then it will be created.
from pci import algo
from pci.api import datasource as ds
# FILE DISCOVERY AND SETUP
start = time.time()
full_ref_insar_image = r"C:\Phoenix_RS2_UltraFine12_20080504_HH_SLC\product.xml"
full_dep_insar_image = r"C:\Phoenix_RS2_UltraFine12_20080528_HH_SLC\product.xml"
external_dem = r"C:\dem.pix"
external_dem_projection = "LONG/LAT D000" # Projection information for the provided external DEM
low_coherence_mask = r"" # optional parameter - used if you want to mask out low coherence areas (i.e. waterbodies)
mask_segment =  # optional parameter - must be set if low_coherence_mask is set
outdir = r"C:\Defo_Output"
# PROCESSING WINDOW
# IT IS STRONGLY RECOMMENDED TO SELECT AN INPUT WINDOW FOR EXPERIMENTATION.
# PROCESSING THE ENTIRE DATASET CAN TAKE QUITE SOME TIME. THIS PROCESSING WINDOW WILL BE USED IN SARINGEST
ul_coord = [391639.993, 3690726.202]
lr_coord = [399629.004, 3687211.037]
# CHECK FOR OUTPUT DIRECTORY - CREATE IT IF REQUIRED
if not os.path.exists(outdir):
Helper Function – create_proc_window
A helper function is defined at the beginning of the script - create_proc_window. This function converts the user provided UL and LR geocoded coordinates to raster (pixel/line) coordinates and calculates the processing window required for the SARINGEST – dbiw parameter.
Note that the method for calculating the processing window for Sentinel-1 TOPS imagery differs slightly from other sensor data (ex. Radarsat-2). For Sentinel-1 TOPS data the processing window is calculated differently depending on whether the data is ascending or descending. If you are processing Sentinel-1 TOPS data you can check the corresponding scripts from the data download.
# HELPER FUNCTION TO CONVERT GEOCODED PROCESSING WINDOW TO PIXEL/LINE VALUES FOR EACH IMAGE
def create_proc_window(image, ul, lr):
This function converts the provided UL and LR geocoded coordinates to raster coordinates and calculates the
processing window required in saringest - dbiw
:param image: Raw SAR Image
:param ul: Upper Left Coordinate
:param lr: Lower Right Coordinate
:return: a list of the processing window information - UL (raster), LR (raster), offset x, offset y.
with ds.open_dataset(image) as dataset:
in_geocode = dataset.geocoding
ul_raster = in_geocode.map_to_raster(ul, ul)
lr_raster = in_geocode.map_to_raster(lr, lr)
offset = [lr_raster - ul_raster, lr_raster - ul_raster]
proc_wind = [int(ul_raster), int(ul_raster), abs(int(offset)), abs(int(offset))]
Ingest SAR Images – SARINGEST
User Manual – Pg. 21
SARINGEST imports a SAR dataset into the PCIDSK (PCI proprietary format) image database format. You can import the entire image, or you can specify a rectangular subset window of data to import. You can also specify the data calibration type.
It is very strongly recommended to use a processing window when running the InSAR process. As previously mentioned, you will set the geographic coordinates of your subset window. They are then converted to the correct processing window for the SARINGEST algorithm.
The processing window is calculated for both the reference and dependent images and both are imported into pix format.
# INGEST SAR IMAGES - PROCESSING WINDOW FROM ABOVE IS USED
step = 1 # INSERTS A PREFIX TO THE FILES CREATED BY THE WORKFLOW
window1 = create_proc_window(full_ref_insar_image, ul_coord, lr_coord)
window2 = create_proc_window(full_dep_insar_image, ul_coord, lr_coord)
ref_insar_image = os.path.join(outdir, "%s_subset_SAR1.pix" % step)
algo.saringest(fili=full_ref_insar_image, calibtyp='SIGMA', filo=ref_insar_image, dbiw=window1)
dep_insar_image = os.path.join(outdir, "%s_subset_SAR2.pix" % step)
algo.saringest(fili=full_dep_insar_image, calibtyp='SIGMA', filo=dep_insar_image, dbiw=window2)
step = step + 1
Co-register Dependent to Reference – INSCOREG
User Manual – Pg. 22
INSCOREG automatically resamples and co-registers the dependent file to the reference file. The module automatically acquires control points, removes outliers, and resamples the dependent file to match one-to-one with the reference file. The metadata of the resampled file is updated to reflect the modified geolocation information.
# INSCOREG - COREGISTER THE DEPENDENT IMAGE TO THE REFERENCE IMAGE AND CREATE A NEW COREGISTETED OUTPUT
inscoreg_output = os.path.join(outdir, "%s_INSCOREG.pix" % step)
algo.inscoreg(filref=ref_insar_image, fili=dep_insar_image, filo=inscoreg_output)
step = step + 1
Raw Interferogram – INSRAW
User Manual – Pg. 23
INSRAW produces the raw interferogram using the specified reference file and resampled dependent file generated by the INSCOREG module. The user has the option to filter the interferogram using an adaptive Lee filter by specifying a filter size (odd) in pixels for the FLSZ parameter. The Lee adaptive filter preserves sharp edges and generally provides a “smoother looking” interferogram. If a filter size is specified it must be five or bigger. In this tutorial we will be using a filter size of 15.
# INSRAW - GENERATE THE RAW INTERFEROGRAM USING THE REFERENCE IMAGE AND THE COREGEISTERED DEPENDENT IMAGE
insraw_output = os.path.join(outdir, "%s_INSRAW.pix" % step)
algo.insraw(filref=ref_insar_image, fili=inscoreg_output, filo=insraw_output, flsz=)
step = step + 1
You can open the interferograms in Focus. For an interferometric file containing complex values (designated in Focus as C32R for complex 32 bit reals, or as C16U for complex 16 bit unsigned integers) the default in Focus is to display the layer as intensity. The data interpretation of the complex layer can be changed by right clicking on the file in the MAPS tab and selecting Data Interpretation and selecting one of the following options Intensity (default), Amplitude (which represents coherence), Phase in Radians, Phase in Degrees, Decibel, Real Part, Imaginary Part. The display will now interpret the complex layer using the selected option. More information on displaying interferograms in Focus is available from the User Manual – Pg. 29.
Flat Earth and Topographic Correction – INSTOPO
User Manual – Pg. 33
The INSTOPO module extracts the reference and dependent file orbit data and calculates the required phase adjustments to correct for systematic flat earth effects projected to the specified earth model (which is generally the WGS-84 ellipsoid). Whether or not you specify a DEM when running this algorithm determines if a deformation or topographic interferogram is created. The InSARContent metadata tag of the output file will is updated to either DeformationInterferogram or TopographicInterferogram.
No DEM required - When no DEM is specified, only the systematic flat earth phase is removed and the output contains fringes due to topographic effects. These fringes can be converted to elevation using the INSUNWRAP module. Residual fringes (if any) due to orbital error can be corrected using INSADJUST.
# INSTOPO - COMPENSATE FOR FLAT EARTH. NOTE THE ABSENCE OF A DEM
# INDICATES WE WISH TO KEEP THE TOPOGRAPHIC EFFECTS TO GENERATE A DEM
instopo_output = os.path.join(outdir, "%s_Instopo.pix" % step)
step = step + 1
DEM is required - The DEM elevations are automatically adjusted to compensate for the nominal processing elevation which is extracted from the metadata. The modified elevations are projected to the reference file slant range prior to calculating the topographic phase correction.
# INSTOPO - COMPENSATE FOR FLAT EARTH AND TOPOGRAPHIC EFFECTS.
# THE INCLUSION OF A DEM INDICATES THAT THE DEFORMATION WORKFLOW IS BEING RUN
instopo_output = os.path.join(outdir, "%s_INSTOPO.pix" % step)
algo.instopo(fili=insraw_output, filedem=external_dem, filo=instopo_output)
step = step + 1
Orbital Correction and Residual Fringe Removal – INSADJUST
User Manual – Pg. 35
The INSADJUST module automatically removes low frequency residual phase caused by orbital drift and/or baseline errors. The module iteratively estimates both the horizontal and vertical spectral residuals and applies the correction based upon the viewing geometry for the interferometric pair.
# INSADJUST - COMPENSATE FOR ORBITAL EFFECTS ON THE INTERFEROGRAM
insadjust_output = os.path.join(outdir, "%s_Insadjust.pix" % step)
algo.insadjust(fili=instopo_output, dbic_ocs=, filo=insadjust_output, maxiter=)
step = step + 1
Modified Goldstein Filter – INSMODGOLD
INSMODGOLD applies a sliding 3x3 spatial boxcar filter to pre-filter the data prior to estimating the fringe frequency. The pre-filtered spatial data is converted to the spectral domain using the FFT size defined by FFTSIZE. Applying larger window sizes will reduce the interferometric noise but may destroy phase continuity in areas with a high number of fringes.
INSMODGOLD is normally applied after the topographic phase correction (INSTOPO) or after the orbit adjustment (INDSADJUST) to enhance interferometric fringes and speed up the following phase unwrapping (INSUNWRAP).
# INSMODGOLD - APPLY A MODIFIED GOLDSTEIN FILTER TO ENHANCE INFEROMETRIC FRINGES
insmodgold_output = os.path.join(outdir, "%s_Insmodgold.pix" % step)
algo.insmodgold(fili=insadjust_output, fftsize=, filo=insmodgold_output)
step = step + 1
Phase Unwrapping – INSUNWRAP
User Manual – Pg. 38
INSUNWRAP unwraps phases of SAR interferograms by using SNAPHU, an open-source, third-party software program. SNAPHU runs in either topographic or deformation mode. INSUNWRAP reads the InSARContent metadata tag that was set by INSTOPO to determine which mode should be run.
The unwrapped interferogram is written to the first channel in the output file. A second 32R channel is also created, with its content based on the InSARContent metadata tag. If InSARContent is set to TopographicInterferogram, the channel contains a real-valued digital elevation model (DEM), in meters. If the InSARContent is set to DeformationInterferogram, the output contains a real-valued deformation map, in meters.
In this example the tilesize parameter is set to 4000, 4000 so the interferogram is processed by Snaphu in a single tile. Unwrapping the interferogram as a single tile will help to avoid discontinuities between the tiles. However, for large subsets (15000 x 15000) and/or low RAM (ex. 8GB) unwrapping in a single tile may be unpractical as it will significantly increase the processing time.
** Note ** Make sure that you have extracted the SNAPHU files as outlined below:
The INSUNWRAP algorithm requires third-party software, SNAPHU, to be installed.
To install SNAPHU on Windows, go to https://github.com/PCIGeomatics/snaphu, download snaphu_win.zip, and then extract the files to the thirdparty\snaphu folder of your Geomatica installation.
To install SNAPHU on Linux, go to https://github.com/PCIGeomatics/snaphu, download snaphu_lx.tar.gz, and then extract the files to the thirdparty/snaphu directory of your Geomatica installation.
# INSUNWRAP - REQUIRES THE SNAPHU MODULE, CONVERTS THE INTERFEROGRAM INTO DISPLACEMENT
insunwrap_output = os.path.join(outdir, "%s_Insunwrap_Inscaldefo.pix" % step)
algo.insunwrap(fili=insmodgold_output, dbic=, filo=insunwrap_output, initalgo='MST', procmode='R',
Generate Calibration Points – NUWRIT/VREAD
User Manual – Pg. 44 (In DEMADJUST section)
Before running DEMADJUST or INSCALDEFO you require a vector layer which contains the geocoded position and elevation of the known adjustment points (DEMADJUST) or stationary points (INSCALDEFO).
NUMWRIT is run on the user provided DEM to extract the geographic coordinates and elevation information for points throughout the DEM file. These points are sampled from the DEM based on the sampling parameter value – [15, 15]. These points are the known points used as adjustment points or stationary points. A text file is created by NUMWRIT that contains these points.
DEMADJUST/INSCALDEFO requires a vector layer containing the adjustment/stationary points. VREAD will be run to convert the text file output from NUMWRIT to a vector segment in the INSUNWRAP output file. We need to make sure that the vecunit parameter is set correctly, based on the projection of the user provided DEM. This was set as the external_dem_projection variable at the beginning of the script.
# NUMWRIT - TAKES AN EXTERNAL DEM AND CREATES A POINT FILE WITH ELEVATION VALUES TO BE USED IN INSCALDEFO
numwrit_output = os.path.join(outdir, "Elevation_points." + "txt")
algo.numwrit(file=external_dem, dbic=, numform="GEOCOORD,Suppress", sampling=[15, 15], tfile=numwrit_output)
# VREAD - READS IN THE ELEVATION POINTS FILE AND CREATES A VECTOR CHANNEL, REQUIRED FOR THE ADJUSTMENT OF THE DEM
vector_file_transfer = numwrit_output
algo.vread(filv=vector_file_transfer, vecunit=external_dem_projection, file=insunwrap_output, dbsn="Stability_Points")
Add New Channel – PCIMOD
A new 32R channel needs to be added to the INSUNWRAP file. This channel will then be used as the output channel for DEMADJUST/INSCALDEFO. When pciop is set to “ADD”, the pcival option specifies the number of each channel to add. The channels are listed in the order: [8U, 16S, 16U, 32R]. We will set this parameter to [0, 0, 0, 1] as we want to add a single 32R channel.
# PCIMOD - ADD 32R CHANNEL, REQUIRED AS OUTPUT CHANNEL BY INSCALDEFO
algo.pcimod(file=insunwrap_output, pciop='add', pcival=[0, 0, 0, 1])
Topographic: Adjustment of Topographic Products – DEMADJUST
User Manual – Pg. 44
DEMADJUST applies a smooth correction to the unwrapped topographic output produced by INSUNWRAP. The unwrapping process introduces a random offset to the topographic product. This module adjusts the raster DEM to minimize the difference between the unwrapped output and the control point elevations stored in the vector layer. The adjmethod parameter specifies the adjustment algorithm that should be used. This controls which points should be used and if there should be any limits on magnitude of the adjustment that will be applied. In this tutorial we will use “UNDERSIGMA 2.0” - Only use calibration points with absolute values of differences under 2.0 sigma.
# DEMADJUST - ADJUSTS THE ELEVATIONS BASED ON REFERENCE ELEVATION DATA
step = step + 1
algo.demadjust(filv=insunwrap_output, dbvs=, fldnme="ZCOORD", filedem=insunwrap_output,
dbec=, dboc=, adjmethod="UNDERSIGMA 2.0")
Topographic: Convert DEM from Slant Range to Geocoded – GEOCODEDEM
User Manual – Pg. 47
GEOCODEDEM converts the digital surface model from slant range to geocoded output. The user has control over the sample size and projection of the output. The process is similar to orthorectification with the exception that the elevation values used as input are extracted from the available file.
# GEOCODEDEM - CONVERT THE OUTPUTS FROM SLANT TO GROUND RANGE
geocodedem_output = os.path.join(outdir, "%s_Geocoded_DEM.pix" % step)
algo.geocodedem(mfile=insunwrap_output, dbec=, failvalu=[-100], backelev=[-150], demedit="NO",
filo=geocodedem_output, mapunits=output_dem_projection, upleft=ul_coord, loright=lr_coord,
end = time.time()
print("Time to complete: ", end - start)
Deformation: Calibrate unwrapped deformation – INSCALDEFO
User Manual – Pg. 48
INSCALDEFO adjusts unwrapped displacement values to zero at all points stored in a vector layer. Typically, the points are stable ground locations (targets) that have no known movements. The input raster is created by unwrapping the interferogram in deformation mode. Deformation interferograms to calibrate can be either orthorectified or in the raw image geometry.
# INSCALDEFO - ADJUST UNWRAPPED DISPLACEMENT VALUES TO ZERO AT ALL CALIBRATION POINTS.
algo.inscaldefo(file=insunwrap_output, dbic=, dboc=, filv=insunwrap_output, dbvs=, calwin=[3, 3],
end = time.time()
print("Time to complete: ", end - start)
Deformation: Further Processing - ORTHO, INSSTACK
User Manual - Pg. 52
INSCALDEFO is the final step in the deformation workflow when processing a single pair of images. Once a number of calibrated deformation products have been generated, it is possible to stack them in chronological order for temporal analysis using the INSSTACK algorithm. Before stacking the deformation products, it is necessary to co-register all of the data sets of interest using ORTHO. The deformation maps may be derived from one or more reference images, but prior to stacking, they must be perfectly co-registered with identical image characteristics (number of pixel and lines), projection and sample sizes.