University of Alaska Fairbanks
Department of Geosciences

GEOS F493 / F693 - Geodetic Methods

Lab 9: Processing w/ gmtSAR: SBAS

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on LuNGS for some things today!

Introduction

Last week you processed some Sentinel 1 data and prepared them for time series analysis. Today, we'll run SBAS on this (short baseline analysis). I instructed you to select interferogram pairs that have short temporal and spatial baselines in the baseline plot. I didn't give very precise instructions, so we all can end up with fairly different results (somewhat by design). Let's try this anyway and see what we all get. This lab is shorter as I am on travel and you should use the remaining time to work on your project.

1. Preparation

Change to your working directory to that of the previous week, make a working directory that you can call whatever you want, but I'll refer to it as SBAS and make sure all the files you need are there. 

    $> cd ~/geodesy
    $> mkdir lab09
    $> cd lab09
    $> ln -s ../lab08/intf_all intf_all
    $> ln -s ../lab08/topo topo
    $> ln -s ../lab08/raw raw
    $> mkdir SBAS

A directory listing should look similar to this: 

    $> ls -l
total 0
lrwxrwxrwx. 1 roon users 17 Oct 30 10:35 intf_all -> ../lab08/intf_all
lrwxrwxrwx. 1 roon users 12 Oct 30 10:35 raw -> ../lab08/raw
drwxr-xr-x. 2 roon users 10 Oct 30 10:38 SBAS
lrwxrwxrwx. 1 roon users 13 Oct 30 10:35 topo -> ../lab08/topo

2. Setup SBAS

Go into the SBAS directory and make sure you can see the baseline_table file in the raw directory from there: 

    $> cd SBAS
    $> ls ../raw | grep base
	baseline.ps
	baseline_table.dat

Now, you're in the SBAS directory. Run sbas, which is the analysis code that we'll use today: 

	$> sbas

 USAGE: sbas intf.tab scene.tab N S xdim ydim [-atm ni] [-smooth sf] [-wavelength wl] [-incidence inc] [-range -rng] [-rms] [-dem]

 input: 
  intf.tab             --  list of unwrapped (filtered) interferograms:
   format:   unwrap.grd  corr.grd  ref_id  rep_id  B_perp 
  scene.tab            --  list of the SAR scenes in chronological order
   format:   scene_id   number_of_days 
   note:     the number_of_days is relative to a reference date 
  N                    --  number of the interferograms
  S                    --  number of the SAR scenes 
  xdim and ydim        --  dimension of the interferograms
  -smooth sf           --  smoothing factors, default=0 
  -atm ni              --  number of iterations for atmospheric correction, default=0(skip atm correction) 
  -wavelength wl       --  wavelength of the radar wave (m) default=0.236 
  -incidence theta     --  incidence angle of the radar wave (degree) default=37 
  -range rng           --  range distance from the radar to the center of the interferogram (m) default=866000 
  -rms                 --  output RMS of the data misfit grids (mm): rms.grd
  -dem                 --  output DEM error (m): dem.grd 

 output: 
  disp_##.grd          --  cumulative displacement time series (mm) grids
  vel.grd              --  mean velocity (mm/yr) grids 

 example:
  sbas intf.tab scene.tab 88 28 700 1000

This tells you that you will need to create 2 files: intf.tab and scene.tab. The latter contains all your interferograms and the number of days with respect to some reference date. That information is in ../raw/baseline_table.dat. My first line in that file looks like this: 

s1b-iw3-slc-vv-20170515t005941-20170515t010006-005602-009d02-006 2017134.0414552251 1229 0.000000000000 0.000000000000

The scene id is 2017134 (2nd column, before the decimal point); the number of days is 1229 making the reference date 2014-01-02 (year of Sentinel 1A launch). This turns into the first line in scene.tab: 

	2017134 1229

Do this for all the other scenes identified in ../raw/baseline_table.dat and save everything to scene.tab.

The second file intf.tab is a tad more involved, but not too bad. For each scene it wants:

The latter is a bit of a tough nut, but not too bad. The information is in the last column of baseline_table.dat, but you would have to do some math subtracting the right baselines from one another as everything in there is given with respect to the first line. I provide you with make_Bperp.sh, which you should copy into your SBAS directory: 

	$>cp /InSAR/make_Bperp.sh $HOME/geodesy/lab08/SBAS/
	$>./make_Bperp.sh

This creates a file bperp.txt in your current directory, which has in it's columns:ref_id rep_id Bperp. The first line for my case looks like this: 

	2017134 2017158 33.993

You should realize that this is already 3/5 of the intf.tab file! All you need to do now is add the paths to the two grid files for each interferogram to each line such that intf.tab (which could be named anything) contains all the necessary information. My first line looks like this: 

	../intf_all/2017134_2017158/unwrap.grd  ../intf_all/2017134_2017158/corr.grd    2017134 2017158 33.993

Note that the path to the .grd files should be with respect to the location from where you run sbas (or you could try using absolute paths.).

Now you'll need to determine 4 other parameters: N, S, xdim, ydim to run sbas:

Note that you may have to swap out the interferogram id if you didn't interfere 2017138 with 2017158. With these values you can now run sbas: 

	$>sbas intf.tab scene.tab $N $S  $xdim $ydim -smooth 1.0 -wavelength 0.0554658 -incidence 30 -range 800184.946186 -rms -dem

Note that we're setting the wavelength, radar wave incidence angle and range to values for S1A, the defaults are for ALOS-1. We also want to look at the quality of fit and DEM error, hence the last two flags. The result you're interested in for now is vel.grd; note that this is in radar coordinates, but you should be able to transform into lat-lon using (all characters are important): 

	$>ln -s ../topo/trans.dat .
	$>ln -s ../intf_all/2017134_2017158/gauss_*
	$>proj_ra2ll.csh trans.dat vel.grd vel_ll.grd
	$>gmt grd2cpt vel_ll.grd -T= -Z -Cjet > vel_ll.cpt
	$>grd2kml.csh vel_ll vel_ll.cpt

Note that proj_ra2ll.csh requires the size of the pixel footprint, which it gets from the gauss_* file ( whichever kind of filtering you used in the interferogram generation), that's why we link it here.

To get an idea of the magnitude of the velocities, find the z_min and z_max fields. The units are given in the SBAS help. 

	$> grdinfo vel_ll.grd

Deliverables: (submit via canvas as zip archive!)

rgrapenthin <at> alaska <dot> edu | Last modified: November 06 2019 02:08.