New Mexico Tech
Earth and Environmental Science

ERTH 455 / GEOP 555 - Geodetic Methods

Lab 9: Processing w/ gmtSAR: SBAS

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on redoubt 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 $GEOP555/lab08
    $> mkdir SBAS
    $> cd SBAS
    $> ls ../raw | grep base
	baseline.ps
	baseline_table.dat

2. Setup SBAS

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

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 ../baseline_table.dat. My first line in that file looks like this: 

	S1A20170515_ALL_F3 2017134.0414552251 1229 0.000000 0.000000

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 /data/GEOP555/make_Bperp.sh $GEOP555/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: 

	$>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.

Deliverables: (submit via canvas as zip archive!)

rg <at> nmt <dot> edu | Last modified: October 25 2017 15:08.