GAMMA DIFF/GEO Software: Reference Manual

Terrain Geocoding (mk_geo, mk_geo_data)

Perl scripts: mk_geo, mk_geo_data

NAME:

mk_geo - Terrain geocoding of SAR images with lookup table refinement and transform DEM --> SAR Range-Doppler Coordinates (RDC)

mk_geo_data - Terrain geocoding of data sets in SAR slant-range coordinates

Description

This is a Quick-start guide on how to use the Perl scripts mk_geo and mk_geo_data to terrain geocode a SAR image, generate a DEM in SAR slant-range coordinates (RDC), and terrain geocode a SAR data that is in the same geometry as the SAR image that was terrain geocoded. These scripts are run from the command line. Think of theses scripts as as assistants to help you automate the geocoding task. mk_geo can perform ellipsoid coding of a SAR image, if no DEM is provided.

mk_geo	terrain geocode an MLI SAR image given a DEM, such as from SRTM, transform that DEM in SAR range-Doppler coordinates (RDC)
mk_geo_data	terrain geocode a SAR data that is in the same geometry as the MLI image that was terrain geocoded.

The processing of geocoding an image involves generating a lookup table that relates the SAR image coordinates to the DEM coordinates. A DEM in the Gamma processing paradigm consists of a text format DEM parameter file and the DEM data itself. The DEM parameter file contains information on the DEM projection. size, pixel spacing, and datum. A number projections are supported including Equiangular geographic (EQA), Universal transverse Mercator (UTM), and Polar stereographic. Often we use the SRTM DEM as input. These data are in EQA format with a height reference of the local geoid. The SRTM data are geoidal heights, relative to the local mean sea-level, rather than to a specific ellipsoid. The process of transforming the SRTM DEM into a UTM or other projection is not described here.

The lookup table is initially generated using just the image processing parameters and the DEM information. Invariably there is some residual shift that must be corrected to take into account inaccuracies in SAR processing or timing. This is done by a process of lookup table refinement. A simulated SAR image is generated from the the DEM that should match the actual SAR image. Offsets between the simulated and actual SAR image are measured and used to generate an offset polynomial model for the error. This polynomial model is used to adjust the lookup-table values to remove the residual geolocation error. In most cases only a single constant range and azimuth offset is requied for accurate geocoding The mk_geo script takes you through these different steps. Programs in the DIFF/GEO package are used to use the lookup table to resample the SAR image into the DEM geometry (geocode_back) and to transform the DEM into the SAR geometry (geocode).

1. Terrain geocoding with mk_geo

The arguments of mk_geo can be seen by typing mk_geo from the command line. You may need to copy mk_geo into your current directory.

*** ./mk_geo

  *** Copyright 2011, Gamma Remote Sensing, v3.3 28-Apr-2011 clw
  ***

  *** Terrain geocoding of SAR images with lookup table refinement
  and transform DEM --> SAR Range-Doppler Coordinates (RDC)
  ***

  

  usage: ./mk_geo <MLI> <MLI_par> <DEM>
  <DEM_par> <DEM_seg> <DEM_seg_par>
  <GEO_dir> <scene_id> <post> <mode>
  [ls_mode] [r_ovr] [n_ovr] [rad_max] [rlks] [azlks] [thres] [rpos]
  [azpos] [roff] [azoff] [r_patch] [az_patch] [nr] [naz]

  

     
  MLI          (input)
  MLI SAR image (including path)

      MLI_par      (input)
  ISP image parameter file of the MLI image (including path)

     
  DEM          (input)
  DEM in desired output projection (including path)

      DEM_par      (input)
  DEM parameter file (including path)

      DEM_seg      (output)
  DEM segment for output image products (including path)

      DEM_seg_par  (input/output) DEM parameter
  file for output image products (including path), regenerated each
  time

      GEO_dir     
  directory for output images, lookup tables and DEM products

      scene_id     scene name to
  identify output files

     
  post         output image
  sample spacing in meters or degrees for geographic (EQA)
  projection

     
  mode         mk_geo
  processing mode:   

                    
  0: generate initial lookup table, simulated SAR image, and DEM
  segment parameters

                    
  1: generate initial lookup table and simulated SAR image using
  existing DEM segment to determine image bounds

                    
  2: measure initial offset between simulated SAR image and actual
  SAR image

                    
  3: perform refinement of lookup table by offset measurement with
  respect to the simulated SAR image

                    
  4: update lookup table and produce terrain geocoded SAR image and
  DEM in SAR range-Doppler coordinates (RDC)

                    
  5: ellipsoid geocoding of the SAR image without a DEM: generate
  new DEM segment
  parameters         


                    
  6: ellipsoid geocoding of the the SAR image without a DEM: use
  existing DEM segment
  parameters         


      ls_mode     
  algorithm selection in gc_map for regions of layover, shadow, or
  DEM gaps:

                    
  0: set to (0.,0.)

                    
  1: linear interpolation across these regions

                    
  2: use actual value (default)

                    
  3: nearest neighbor thinned (nn-thinned)

     
  r_ovr        range
  over-sampling parameter for gc_map parameter ls_mode = 3
  (nn-thinned) (default: 2)

     
  n_ovr        interpolation
  oversampling factor in geocode (default = 2.0)

      rad_max      maximum
  interpolation search radius in geocode (enter - for default:
  4*n_ovr)

     

      NOTE: n_ovr and rad_max are parameters used by
  the program geocode for transformation of the simulated image

            and DEM
  into SAR geometry. The parameters rlks, azlks, thres, rpos,
  azpos, roff, azoff are used 

        for estimation of the initial offset of
  the SAR image with respected to the simulated SAR image.

              
       

     
  rlks         number of
  range looks for the initial offset estimate (default: 1)

     
  azlks        number of azimuth
  looks for the initial offset estimate (default: 1)

     
  thres        SNR threshold for
  offset estimates (default: 10)

     
  rpos         range
  position for initial offset (enter - for default)

     
  azpos        azimuth position
  for initial offset (enter - for default)

     
  roff         initial
  range offset estimate (enter - for current value in DIFF_par
  file)

     
  azoff        initial azimuth
  offset estimate (enter - for current value in DIFF_par file)

      r_patch      range
  patch size for offset estimation (default: 256 samples)

      az_patch     azimuth patch
  size for offset estimation (default: 256 lines)

     
  nr          
  number of range patches for offset estimation (default: 20)

     
  naz          number
  of azimuth patches for offset estimation (default: 20)

      -s scale  (option) set image display
  scale factor (default: 0.6)  

      -e exp    (option) set image
  display exponent (default: 0.4)

This is not nearly as daunting as it seems. mk_geo is run in a series of steps specified by the mode command line parameter. As input to running mk_geo you need to provide all command-line inputs with a <>. To understand the meaning of the command line parameters ls_mode, r_over, n_over and rad_max, consult the documentation on the gc_map program in the DIFF/GEO package. Nominally set ls_mode to 2 for differential SAR interferometry applications. These files are in the GEO_dir that you specify on the command line.

mk_geo generates a log file of its actions in the GEO_dir that be examined using any text editor. All script commands and the screen output from running the Gamma programs are captured in the log files. There are successive log files generated for each step in the mk_geo processing sequence with names like mk_geo_n.log. where n is the mode number.

Mode 0 or 1: generate initial geocoding lookup table and simulate SAR image from DEM

./mk_geo rmli_2_10/rmli_2_10.ave
  rmli_2_10/rmli_2_10.ave.par lan_utm.dem lan_utm.dem_par
  lan_utm_seg.dem lan_utm_seg.dem_par geo lan 40 0 2

The first step is to generate the initial lookup table. Modes 0, or 1 run the gc_map program and generate the initial lookup table, and a simulated SAR image in the SAR coordinates. To do this you can choose either mode 0, where you do not have a predefined DEM parameter file that describes the area for your geocoded image, or mode 1 where you have such a file available. Let us assume that you do NOT have such a file and that the terrain geocoded image will cover the entire are of your input SAR image. In the processing example that is part of this tutorial a small region in the California Mojave desert near the region of the Landers earthquake has been selected. A tandem pair of SLCs produced by ERS-1 and ERS-2 acquired 19960421 and 19960422 respectively were co registered, detected and averaged. The MLI image is located in the rmli_2_10 directory and is called rmli_2_10.ave. The associated MLI image parameter file is rmli_2_10.par and is the same as the image parameter file from the 19960421 ERS-1 image. The DEM is in the UTM projection and was generated from the SRTM global DEM.

The input DEM is in UTM and can be displayed using:


 disdem_par lan_utm.dem lan_utm.dem_par

The first step is generate the look up table to geocode the input image rmli_2_10.ave with parameter file rmli_2_10.ave.par using the DEM data in lan_utm.dem_par with parameter file lan_utm.dem_par:



  mk_geo rmli_2_10/rmli_2_10.ave rmli_2_10/rmli_2_10.ave.par
  lan_utm.dem lan_utm.dem_par lan_utm_seg.dem lan_utm_seg.dem_par
  geo lan 40 0 2

The portion of the the DEM covering the SAR image will be resampled to 40 meters and stored in the file lan_utm_seg.dem with parameter file lan_utm_seg.dem_par.
The output products will be in the geo directory .

The following commands are executed by the script in Mode 0 to accomplish this:

#generate lookup table and simulated SAR image lan_0.sim

gc_map rmli_2_10/rmli_2_10.ave.par - lan_utm.dem_par
  lan_utm.dem lan_utm_seg.dem_par lan_utm_seg.dem
  geo/lan_0.map_to_rdc 1.25 1.25 geo/lan_0.sim - - - - - - 8 2
  2

Note that the script calculates an oversampling factor of 1.25 to produce a resampled DEM at 40 meters from the 50 meter posting in the input DEM.



  #transform simulated SAR image lan_0.sim using the lookup table
  lan_0.map_to_rdc into RDC coordinates

  geocode geo/lan_0.map_to_rdc geo/lan_0.sim 1266 geo/lan_0.sim_rdc
  1148 574 2 0 1 1 2 8

  

  #generate raster image of simulated image in RDC and display on
  the terminal

  raspwr geo/lan_0.sim_rdc 1148 1 0 1 1

  dis2ras geo/lan_0.sim_rdc.ras rmli_2_10/rmli_2_10.ave.ras

You will see a display window pop up with the simulated image and the SAR image together. You can flicker between the images using the center mouse button.
If Mode 1 had been used, the lan_utm_seg.dem_par would have had to exist.

`Mode 2: Lookup table refinement, initial offset estimate using a simulated SAR image`

./mk_geo
  rmli_2_10/rmli_2_10.ave rmli_2_10/rmli_2_10.ave.par lan_utm.dem
  lan_utm.dem_par lan_utm_seg.dem lan_utm_seg.dem_par geo lan 40 2
  2

Residual errors in the geocoding are estimated and removed by measuring the offsets between the simulated SAR image produced from the DEM and the actual SAR image. The offset information is stored in a DIFF/GEO parameter file and created using the create_diff_par program in the DIFF/GEO package. The simulated SAR image is calculated in the map (DEM) geometry when running modes 0 or 1. by the call to gc_map. Running mode 2 creates the DIFF/GEO parameter file, and measures the initial offset using theinit_offsetm program. Try the defaults if they do not work, you can specify another region in the image by adding the command line parameters :

    rlks    number of
  range looks for the initial offset estimate (default: 1)

      azlks   number of azimuth looks
  for the initial offset estimate (default: 1)

      thres   SNR threshold for
  offset estimates (default: 10)

      rpos    range position for
  initial offset (enter - for default)

      azpos   azimuth position for initial
  offset (enter - for default)

      roff    initial range offset
  estimate (enter - for current value in DIFF_par file)

      azoff   initial azimuth offset
  estimate (enter - for current value in DIFF_par file)

These are arguments to the init_offsetm program. See the DIFF/GEO documentation on init_offsetm on how to set them to get a valid offset. The commands executed by mode 2 are:

#create DIFF/GEO parameter file using script generate lan.diff.in inputs

create_diff_par rmli_2_10/rmli_2_10.ave.par - geo/lan.diff
  1 < geo/lan.diff.in

  

  #determine initial offset between lan_0.sim_rdc and
  rmli_2_10.ave

  init_offsetm geo/lan_0.sim_rdc rmli_2_10/rmli_2_10.ave
  geo/lan.diff 1 1 574 287 0 0 10 - 1

Mode 3: Estimate offsets between the SAR and simulated image

./mk_geo
  rmli_2_10/rmli_2_10.ave rmli_2_10/rmli_2_10.ave.par lan_utm.dem
  lan_utm.dem_par lan_utm_seg.dem lan_utm_seg.dem_par geo lan 40 3
  2

Note that most cases it is sufficient for getting an accurate scene geocoding only to measure this one offset and go directly to mode 4. In some images with very few features there are only a few areas in the scene that can be used for measuring the offset. In that case, do not run mode 3 of mk_geo. Enter the center of this area on the command line in the MLI pixel coordinates for mode 2 and go directly to mode 4. Running mk_geo mode 3 will estimate the offsets over the entire scene using image patches (256x256) or as specified on the command line by the r_patch and az_patch parameters. The program offset_pwrm measures the offsets over the entire scene

#measure offsets of intensity images

  offset_pwrm geo/lan_0.sim_rdc rmli_2_10/rmli_2_10.ave
  geo/lan.diff geo/lan.offs geo/lan.snr

  

  #calculate range and azimuth offset polynomials with 3 terms
  each

  offset_fitm geo/lan.offs geo/lan.snr geo/lan.diff geo/lan.coffs -
  10 3

  

  In very few cases with sparse dem data or images with no
  contrast, it may not be possible to find enough points to
  estimate the offset polynomial parameters. In this case, select a
  single area in the image with some features and that are also
  visible in the simulated SAR image. Determine the coordinates of
  this point in the SAR image and then run mode 2 above specifying
  the coordinates with rpos and azpos command line parameters. The
  initial offset values are written as the constant values of
  offset polynomials. This should usually be sufficient to geocode
  the image and you can omit running Mode 3.

After this process the DIFF_GEO parameter file is as below, note the range and azimuth offset polynomials:

Gamma DIFF&GEO Processing Parameters

  title: lan

  initial_range_offset:                   
  3

  initial_azimuth_offset:                 
  0

  range_samp_1:                        
  1148

  az_samp_1:                            
  574

  first_nonzero_range_pixel_1:            
  0

  number_of_nonzero_range_pixels_1:    
  1148

  range_pixel_spacing_1:          
  15.809780

  az_pixel_spacing_1:             
  39.907180

  range_samp_2:                        
  1148

  az_samp_2:                            
  574

  first_nonzero_range_pixel_2:            
  0

  number_of_nonzero_range_pixels_2:    
  1148

  range_pixel_spacing_2:          
  15.809780

  az_pixel_spacing_2:             
  39.907180

  offset_estimation_starting_range:      
  48

  offset_estimation_ending_range:      
  1100

  offset_estimation_range_samples:       
  16

  offset_estimation_range_spacing:       
  70

  offset_estimation_starting_azimuth:    
  48

  offset_estimation_ending_azimuth:     
  526

  offset_estimation_azimuth_samples:     
  16

  offset_estimation_azimuth_spacing:     
  31

  offset_estimation_patch_width:        
  256

  offset_estimation_patch_height:       
  256

  offset_estimation_threshhold:       
  10.00

  range_offset_polynomial:         
  3.47071   8.7314e-05  -1.6293e-04  
  0.0000e+00

  azimuth_offset_polynomial:       
  0.33467   6.9837e-05  -6.3799e-04  
  0.0000e+00

  starting_azimuth_line:                  
  0

  map_azimuth_lines:                      
  0

  map_width:                              
  0

  first_map_range_pixel:                  
  0

  number_map_range_pixels:                
  0

  range_looks:                            
  0

  azimuth_looks:                          
  0

  diff_phase_fit:      
  0.00000   0.0000e+00   0.0000e+00  
  0.0000e+00   0.0000e+00  
  0.0000e+00

Mode 4: Update lookup table, generate terrain coded image and DEM in SAR RDC

./mk_geo rmli_2_10/rmli_2_10.ave
  rmli_2_10/rmli_2_10.ave.par lan_utm.dem lan_utm.dem_par
  lan_utm_seg.dem lan_utm_seg.dem_par geo lan 40 4 2

The final steps in terrain geocoding are updating the lookup table, generate the terrain geocoded image from the SAR image in radar coordinates, and the resampling the DEM into SAR coordinates. The resampled DEM is used for simulation of interferograms in the production of 2-pass differential interferograms. At each step, output image products are generated for comparison and display:

#refine lookup table
gc_map_fine geo/lan_0.map_to_rdc 1266 geo/lan.diff geo/lan_1.map_to_rdc 1
geocode geo/lan_1.map_to_rdc geo/lan_0.sim 1266 geo/lan_1.sim_rdc 1148 574 2 0 1 1 2 8

#project simulated SAR image into RDC and compare with original SAR image to confirm refinement worked correctly
geocode geo/lan_1.map_to_rdc geo/lan_0.sim 1266 geo/lan_1.sim_rdc 1148 574 2 0 1 1 2 8
raspwr geo/lan_1.sim_rdc 1148 1 0 1 1
dis2ras geo/lan_1.sim_rdc.ras rmli_2_10/rmli_2_10.ave.ras
dis2ras geo/lan_1.sim_rdc.ras rmli_2_10/rmli_2_10.ave.ras&

#project DEM into SAR range-Doppler coordinates and display
geocode geo/lan_1.map_to_rdc lan_utm_seg.dem 1266 geo/lan_dem.rdc 1148 574 2 0 1 1 2 8
rasshd geo/lan_dem.rdc 1148 40 40 1 0 1 1
disras geo/lan_dem.rdc.ras&

#terrain geocode the SAR image and display
geocode_back rmli_2_10/rmli_2_10.ave 1148 geo/lan_1.map_to_rdc geo/lan_map.mli 1266 818 1 0
raspwr geo/lan_map.mli 1266 1 0 1 1
disras_dem_par geo/lan_map.mli.ras lan_utm_seg.dem_par

The output SAR geocoded image is lan_map.mli, the DEM in SAR coordinates is called lan_dem.rdc, and the simulated image in RDC is lan_1.sim_rdc. There should be only very small residual offsets between lan_1.dm_rdc and the original SAR image. The output images have a geometry as defined in the lan_utm_seg.dem_par parameter file. shown here:

Gamma DIFF&GEO DEM/MAP parameter file

  title: Landers, California

  DEM_projection:     UTM

  data_format:        REAL*4

  DEM_hgt_offset:         
  0.00000

  DEM_scale:              
  1.00000

  width:               
  1266

  nlines:               
  818

  corner_north:  3897100.000   m

  corner_east:    489300.000   m

  post_north:    -40.0000000   m

  post_east:      40.0000000  
  m

  

  ellipsoid_name: WGS 84

  ellipsoid_ra:       
  6378137.000   m

  ellipsoid_reciprocal_flattening:  298.2572236

  

  datum_name: WGS 1984

  datum_shift_dx:             
  0.000   m

  datum_shift_dy:             
  0.000   m

  datum_shift_dz:             
  0.000   m

  datum_scale_m:        
  0.00000e+00

  datum_rotation_alpha:  0.00000e+00   arc-sec

  datum_rotation_beta:   0.00000e+00  
  arc-sec

  datum_rotation_gamma:  0.00000e+00   arc-sec

  datum_country_list Global Definition, WGS84, World

  

  projection_name: UTM

  projection_zone:                
  11

  false_easting:          
  500000.000   m

  false_northing:              
  0.000   m

  projection_k0:           
  0.9996000

  center_longitude:     
  -117.0000000   decimal degrees

  center_latitude:         
  0.0000000   decimal degrees

2. Terrain geocoding additional data files using the mk_geo_data script

An additional script mk_geo_data exists to generate terrain geocoded images of additional data products if mk_geo has been completed successfully. This script extracts values from the different parameter files and generates the call to geocode_back. the parameters are in the MLI_par and DEM_seg_par files.

*** mk_geo_data

  *** Copyright 2005, Gamma Remote Sensing, 1.1 13-Sep-2005 clw
  ***

  *** terrain geocoding of a data set in SAR RDC coordinates
  ***

  

  usage: ./mk_geo_data <MLI_par> <DEM_seg_par>
  <gc> <data> <data_geo> <interp_mode>
  <format>

   
  MLI_par      (input) MLI image parameter
  file with same dimensions as input data

      DEM_seg_par  (input) DEM parameter file
  for output geocoded product

     
  gc          
  (input) geocoding lookup table (map coordinates --> RDC), use
  refined lookup table if available

     
  data         (input)
  input data set with same dimensions as the MLI image

      data_geo     (output)
  output terrain geocoded data set

      interp_mode  interpolation mode:

                    
  0: nearest-neighbor

                    
  1: spline

                    
  2: spline-log

                    
  3: bilinear

                    
  4: bilinear_log

      format      
  input data format:

                    
  0: FLOAT

                    
  1: FCOMPLEX

                    
  2: SUN raster/BMP

                    
  3: UNSIGNED CHAR

                    
  4: SHORT

The MLI_par is used to define the dimensions of the SAR data set. This data set may be an interferometric product such as a deformation map. Also required is the lookup table that has been updated by gc_map_fine to account for geocoding location errors. The user can specify the interpolation mode and data format. Further documentation on these parameters is in the DIFF/GEO documentation of geocode_back. For example to geocode a float format height map in RDC into UTM

mk_geo_data rmli_2_10/rmli_2_10.ave.par
  lan_utm_seg.dem_par geo/lan_1.map_to_rdc
  diff1_2d/19960421_19960422.hgt diff1_2d/19960421_19960422_utm.hgt
  3 0

  

  This will generate a command to geocode_back:

  

  geocode_back diff1_2d/19960421_19960422.hgt 1148
  geo/lan_1.map_to_rdc diff1_2d/19960421_19960422_utm.hgt 1266 818
  3 0