#!/bin/csh -f

# Amplitude dispersion from a stack of aligned SLCs

# Creator: Xiaopeng Tong, David Sandwell

# Date: July 4th 2013

# input: PRM.list

# output: amplitude dispersion

# The principle: first calibrate each images by a weighting

# factor then compute the mean amplitude M_A

# and the standard deviation of the amplitude sig_A

# then take the ratio of the M_A and sig_A

# Reference: Lyons and Sandwell, Fault creep along the southern San Andreas from InSAR, permanent scatterers, and stacking, J. Geophys. Res., 108 (B1), 2047, doi:10.1029/ 2002JB001831 2003

alias rm 'rm -f'

if ($#argv != 3) then

echo ""

echo "Usage: dispersion.csh PRM.list scatter.grd <rng0>/<rngf>/<azi0>/<azif>"

echo ""

echo " Amplitude dispersion from a stack of aligned SLCs"

echo ""

echo " PRM.list -- a list of PRM files of the aligned SLC"

echo " scatter.grd -- output scattering coeffient grid"

echo " <rng0>/<rngf>/<azi0>/<azif> — region of interests in radar coordinates"

echo ""

exit 1

endif

if (! -e $1) then

echo ""

echo "no input file found: $1"

echo ""

exit 1

endif

set list = $1

set outgrd = $2

set region = $3

#set filter = $GMT5SAR/filters/gauss1x1

set filter = $GMT5SAR/filters/gauss5x3

set namearray =

set weightarray =

echo ""

echo "Start -- compute amplitude dispersion"

echo ""

#

# conv to get amplitude of the SLCs and compute their sum and average

#

echo "compute amplitude from SLCs ..."

@ num = 1

foreach prm (`cat $list`)

set name = `grep input_file $prm | awk '{print $3}' | sed 's/\.raw//'`

set namearray = ($namearray $name)

if (-e $prm) then

if (! -e $name".grd") then

echo "running conv on file $name"

conv 1 1 $filter $prm $name".grd"

gmt grdmath $name".grd" FLIPUD = $name".grd"

endif

# gmt grd2cpt $name".grd" -Cgray -Z -D > cpt

gmt makecpt -Cgray -T1e-8/1e-7/1e-8 -Z > a.cpt

gmt grdimage $name".grd" -Ca.cpt -JX6i -P -Xc -Yc -V -R0/6144/0/25200 -Bf2500a5000WSen > $name".ps"

else

echo ""

echo "no PRM file found: $prm"

echo ""

exit 1

endif

gmt grdcut $name".grd" -R$region -Gtmp.grd

mv tmp.grd $name".grd"

if ($num == 1) then

gmt grdmath $name".grd" = sum.grd

else

gmt grdmath sum.grd $name".grd" ADD = sumtmp.grd

mv sumtmp.grd sum.grd

endif

@ num ++

end

@ num --

gmt grdmath sum.grd $num DIV = ave.grd

set avemean = `gmt grdinfo ave.grd -L2 -C | cut -f 12`

echo $avemean

#

#

# compute the M_A term

#

echo ""

echo "compute the M_A term ..."

@ num = 1

foreach name ($namearray)

set ampmean = `gmt grdinfo $name".grd" -L2 -C | cut -f 12`

set weight = `echo $ampmean $avemean | awk '{print $1/$2}'`

echo $weight

set weightarray = ($weightarray $weight)

gmt grdmath $name".grd" $weight DIV = tmp.grd

if ($num == 1) then

gmt grdmath $name".grd" $weight DIV = sum.grd

else

gmt grdmath $name".grd" $weight DIV sum.grd ADD = sumtmp.grd

mv sumtmp.grd sum.grd

endif

@ num ++

end

@ num --

gmt grdmath sum.grd $num DIV = M_A.grd

#rm tmp.grd

#

#

# compute the sig_A term

#

echo "compute the sig_A term ..."

@ num = 1

foreach name ($namearray)

if ($num == 1) then

gmt grdmath $name".grd" $weightarray[$num] DIV M_A.grd SUB SQR = sum2.grd

else

gmt grdmath $name".grd" $weightarray[$num] DIV M_A.grd SUB SQR sum2.grd ADD = sum2tmp.grd

mv sum2tmp.grd sum2.grd

endif

@ num ++

end

@ num --

gmt grdmath sum2.grd $num DIV SQRT = sig_A.grd

#

#

# compute the amplitude dispersion

#

echo ""

echo "compute the scattering amplitude ... "

gmt grdmath sig_A.grd M_A.grd DIV = $outgrd

gmt makecpt -Cgray -T0.1/0.6/0.1 -Z -D > scatter.cpt

gmt grdimage $outgrd -Cscatter.cpt -JX6i -P -X1i -Y3i -V -R$region -Bf5000a10000WSen -K > scatter.ps

gmt psscale -D3/-0.8/5/0.5h -Cscatter.cpt -Ba0.2f0.1:"amplitude dispersion": -O >> scatter.ps

echo ""

echo "Finish -- compute amplitude dispersion"

echo ""

#

# clean up

#

foreach name ($namearray)

# rm $name".grd"

end

rm sum.grd sumtmp.grd sum2.grd sum2tmp.grd ave.grd

rm scatter.cpt M_A.grd sig_A.grd