N.B. The original version of the following script was written by David T. Sandwell for the Salton Sea region. This was adapted by Karen Watson for use in the Los Angeles region.
You will need to adapt the following lines (coordinates and DEM filenames) for your particular region:
dem2xyz santa_ana-w
dem2xyz san_bernadino-w
dem2xyz los_angeles-e
dem2xyz long_beach-e
cat long_beach-e.xyz los_angeles-e.xyz san_bernadino-w.xyz santa_ana-w.xyz > alldata.xyz
blockmedian alldata.xyz -I.001/.001 -R-118.8/-117.1/33.0/34.4 > block.xyz
surface block.xyt -I.001/.001 -R-118.8/-117.1/33.0/34.4 -T.35 -Gtopo.grd -V
# script to make topography for the salton sea area
# since much of the area is below sea level and the
# ermapper compilation does not contain negative elevations,
# one must go back to the original USGS files. In addition
# the USGS elevations are height above the geoid while the
# GPS points are height above the WGS84 ellipsoid.
# Thus we must add the geoid height to the USGS elevations.
#
# KMW: Edited to cover the LA region instead of Salton Sea (6/16/00)
#
# convert the relevant files to xyz
#
#KMW dem2xyz san_bernardino-e
#KMW dem2xyz needles-w
#KMW dem2xyz needles-e
#KMW dem2xyz santa_ana-e
#KMW dem2xyz salton_sea-w
#
dem2xyz santa_ana-w
dem2xyz san_bernadino-w
dem2xyz los_angeles-e
dem2xyz long_beach-e
#
# put all of the files together
#
#KMW cat san_bernardino-e.xyz needles-w.xyz needles-e.xyz santa_ana-e.xyz salton_sea-w.xyz > alldata.xyz
#
cat long_beach-e.xyz los_angeles-e.xyz san_bernadino-w.xyz santa_ana-w.xyz > alldata.xyz
#
# do a blockmedian to prepare for surface
#
#KMW blockmedian alldata.xyz -I.001/.001 -R-116.6/-114.9/33.1/34.5 > block.xyz
blockmedian alldata.xyz -I.001/.001 -R-118.8/-117.1/33.0/34.4 > block.xyz
#
# add a column of geoid height
#
#KMW interp_ship << END
/user2/sandwell/alt/bin/interp_ship << END
block.xyz
1 2 3
.01
block.xyzn
/topex/ftp/pub/global_geoid_2min/geoid.img.9.2
END
#
# add the geoid to the elevation
#
awk '{print $1, $2, $3+$4}' < block.xyzn > block.xyt
#
# now grid the heights and make a gips file
#
#KMW surface block.xyt -I.001/.001 -R-116.6/-114.9/33.1/34.5 -T.35 -Gtopo.grd -V
surface block.xyt -I.001/.001 -R-118.8/-117.1/33.0/34.4 -T.35 -Gtopo.grd -V
/usr/local/src/gips/bin/grd2gips topo.grd topo.gips r4
#
# now reproject the data from lon, lat, topo to rat
#KMW - note that this step takes the *LONGEST* time to run
#
latlon2radar master.head topo.gips > block.rat
#
# regrid the data using tsurface
#
blockmedian block.rat -R00/6144/00/28000 -I4/16 -V > temp.rat
surface temp.rat -R00/6144/00/28000 -I4/16 -T.30 -Gnode.grd -V
grdsample node.grd -T -Gpixel.grd
#
# convert to a gips file
#
/usr/local/src/gips/bin/grd2gips pixel.grd topo_rat.gips r4
#
# now apply a 9 x 9 gaussian filter with a sigma = 2
# to the topogrid. For an azimuth spacing of 16 x 4m = 64 m
# this will be a 0.5 cutoff of 64 x 2 x 5.3 = 680 m
#
ihconv 1 1 /opt/siosar/filters/gauss9x9 < topo_rat.gips > junk.gips
#
#
ihrot -ni junk.gips | ihm -mmaster.head > filtered_topo.gips
#
# filter it again
#
ihconv 1 1 /opt/siosar/filters/gauss9x9 < filtered_topo.gips > match_topo.gips
#
# convert topography to phase and compute x-phase and y-phase
#
topo2phase match_topo.gips phase.gips
#
ihconv 1 1 /opt/siosar/filters/xdir < phase.gips | \
iha "p = p1/2.;" | \
ihm -m"dmin=-.3 dmax=.3" > x_phase.gips
ihconv 1 1 /opt/siosar/filters/ydir < phase.gips | \
iha "p = p1/2.;" | \
ihm -m"dmin=-.3 dmax=.3" > y_phase.gips
#
#clean up
rm *.xyz
rm junk.gips temp.gips
rm *.grd filtered_topo.gips
Back to: InSAR class homepage
Last modified: 4 October 2000