Applications of GMT David Sandwell August 16 2016

Applications of GMT David Sandwell August 16, 2016 • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

Soundings used in SRTM 15/30_PLUS multibeam, singlebeam, grids, . . .

Data Sources multi-beam soundings predicted from satellite gravity single-beam soundings

source # points % flagged (500 m) NGDC_multi % seafloor @ 1 km 127901083 4. 75 5. 134 NOAA_geodas 40897565 11. 06 2. 506 US_multi 51187020 5. 32 2. 219 JAMSTEC_multi 79103040 0. 62 2. 04 SIO_various 40754645 13. 93 1. 325 IBCAO_various 18302390 0. 03 0. 773 GEBCO_various 8950614 0. 17 0. 523 AGSO_grid 12875795 2. 62 0. 503 DNC_points 5878651 1. 23 0. 49 CCOM_grid 10023471 0. 08 0. 195 GEOMAR_grid 18138868 0. 06 0. 181 NGA_single 4415125 19. 44 0. 179 NOAA_grids 6748376 0. 49 0. 162 IFREMER_single 7653537 10. 62 0. 151 3 DGBR_various 5523560 11. 36 0. 112 422089 0. 06 0. 009 NAVO_multi total 438, 775, 829 16. 5 % in V 11 12% in V 7

if ($#argv < 3) then echo " " echo " example: `basename $0` median. xyz surface. Opts unmasked. grd" echo " " exit endif set bxyz set out set opts = $1; shift = "$*" #/bin/rm -rf *out. grd # set B 1 = -180. /58. /90. ; set B 2 = -180. /28. /62. ; set B 3 = -180. /-2. /32. ; set B 4 = -180. /-32. /02. ; set B 5 = -180. /-62. /-28. ; set B 6 = -180. /-90. /-58. ; # set C 1 = -180. /60. /90. ; set C 2 = -180. /30. /60. ; set C 3 = -180. /00. /30. ; set C 4 = -180. /-30. /00. ; set C 5 = -180. /-60. /-30. ; set C 6 = -180. /-90. /-60. ; # # make the blend file # echo B 1. grd -R$C 1 1 > blend. txt echo B 2. grd -R$C 2 1 >> blend. txt echo B 3. grd -R$C 3 1 >> blend. txt echo B 4. grd -R$C 4 1 >> blend. txt echo B 5. grd -R$C 5 1 >> blend. txt echo B 6. grd -R$C 6 1 >> blend. txt # # do all the subgrids # #1 # echo $bxyz -V -fg -bid -I 15 c -A 0. 5 -R$B 1 $opts -GB 1$out blockmedian $bxyz -I 15 c -bid -bod -R$B 1 -V > B 1. xyz surface B 1. xyz -V -fg -bid -I 15 c -A. 50 -Ll-800 -Lu 800 -R$B 1 $opts -GB 1. grd # JJ Becker wrote most of this code # #2 # echo $bxyz -V -fg -bid -I 15 c -A. 707 -R$B 2 $opts -GB 2$out blockmedian $bxyz -I 15 c -bid -bod -R$B 2 -V > B 2. xyz surface B 2. xyz -V -fg -bid -I 15 c -A. 707 -Ll-800 -Lu 800 -R$B 2 $opts -GB 2. grd # #3 # echo $bxyz -V -fg -bid -I 15 c -A. 966 -R$B 3 $opts -GB 3. grd blockmedian $bxyz -I 15 c -bid -bod -R$B 3 -V > B 3. xyz surface B 3. xyz -V -fg -bid -I 15 c -A. 966 -Ll-800 -Lu 800 -R$B 3 $opts -GB 3. grd # #4 # echo $bxyz -V -fg -bid -I 15 c -A. 966 -R$B 4 $opts -GB 4. grd blockmedian $bxyz -I 15 c -bid -bod -R$B 4 -V > B 4. xyz surface B 4. xyz -V -fg -bid -I 15 c -A. 966 -Ll-800 -Lu 800 -R$B 4 $opts -GB 4. grd # #5 # echo $bxyz -V -fg -bid -I 15 c -A. 707 -R$B 5 $opts -GB 5. grd blockmedian $bxyz -I 15 c -bid -bod -R$B 5 -V > B 5. xyz surface B 5. xyz -V -fg -bid -I 15 c -A. 707 -Ll-800 -Lu 800 -R$B 5 $opts -GB 5. grd # #6 # echo $bxyz -V -fg -bid -I 15 c -A. 50 -R$B 6 $opts -GB 6. grd blockmedian $bxyz -I 15 c -bid -bod -R$B 6 -V > B 6. xyz surface B 6. xyz -V -fg -bid -I 15 c -A. 50 -Ll-800 -Lu 800 -R$B 6 $opts -GB 6. grd # # now blend all the files together # grdblend. txt -G$out -R-180/-90/90 -fg -I 15 c -V

S&S 1 min V 18

SRTM 30_PLUS V 11

SRTM 15_PLUS V 1

Applications of GMT • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

grd 2 kml with help from Joaquim #!/bin/csh –f # unset noclobber # # script to convert a grd file to a kml file for Google Earth # if ($#argv < 2 || $#argv > 3) then echo " " echo "Usage: grd 2 kml. csh grd_file_stem cptfile [-R<west>/<east>/<south>/<north>] " echo "Example: grd 2 kml. csh phase. cpt " echo " " exit 1 endif # # set DX = `gmt grdinfo $1. grd -C | cut -f 8` set DPI = `gmt gmtmath -Q $DX INV RINT = ` #echo $DPI gmt set COLOR_MODEL = hsv gmt set PS_MEDIA = letter # gmt grdgradient $1. grd -Ggrad. grd -V -Nt 0. 7 -A 60 if ($#argv == 3) then gmt grdimage $1. grd -Igrad. grd -C$2 $3 -Jx 1 id -P -Y 2 -X 2 -Q -V > $1. ps else if ($#argv == 2) then gmt grdimage $1. grd -Igrad. grd -C$2 -Jx 1 id -P -Y 2 -X 2 -Q -V > $1. ps endif # # now make the kml and png # gmt ps 2 raster $1. ps -W+k+t"$1" -E$DPI -TG -P -S -V -F$1. png rm $1. ps rm grad. grd rm ps 2 raster* rm psconvert* #

Applications of GMT • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

GMT code

Applications of GMT • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

dem 2 topo_ra. csh – Matt Wei if ($#argv < 2) then echo " " echo "Usage: dem 2 topo_ra. csh master. PRM dem. grd" echo " ” endif # set PRF = `grep PRF *. PRM | awk 'NR == 1 {printf("%d", $3)}'` # # use SAT_llt 2 rat # gmt grd 2 xyz --FORMAT_FLOAT_OUT=%lf $2 -s | SAT_llt 2 rat $1 1 -bod > trans. dat # # use an azimuth spacing of 2 for low PRF data such as S 1 A TOPS # if ($PRF < 1000) then gmtconvert trans. dat -o 0, 1, 2 -bi 5 d -bo 3 d | gmt blockmedian -R 0/$XMAX/0/$YMAX -I$rng/2 -bi 3 d -bo 3 d -r -V > temp. rat gmt surface temp. rat -R 0/$XMAX/0/$YMAX -I$rng/2 -bi 3 d -T. 50 -N 1000 -Gpixel. grd -r -V else gmtconvert trans. dat -o 0, 1, 2 -bi 5 d -bo 3 d | gmt blockmedian -R 0/$XMAX/0/$YMAX -I$rng/4 -bi 3 d -bo 3 d -r -V > temp. rat gmt surface temp. rat -R 0/$XMAX/0/$YMAX -I$rng/4 -bi 3 d -T. 50 -N 1000 -Gpixel. grd -r -V endif # # flip to bottom for both ascending and descending passes # gmt grdmath pixel. grd FLIPUD = topo_ra. grd # # plotting # gmt grd 2 cpt topo_ra. grd -Cgray -V -Z > topo_ra. cpt gmt grdimage topo_ra. grd $scale -X 0. 2 i -P -Ctopo_ra. cpt -V > topo_ra. ps #

Applications of GMT • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

# include "gmtsar. h" # include "sbas. h" sbas. c – Xiaopeng Tong int main(int argc, char **argv) { /* Begin: Initializing new GMT 5 session */ if ((API = GMT_Create_Session (argv[0], 0 U, NULL)) == NULL) return EXIT_FAILURE; /* reading in the table files */ read_table_data_ts(API, infile, datefile, gfile, cfile, H, bperp, flag, var, phi, S, N, xdim, ydim, &Out, L, time); for (i=0; i<N; i++) { if ((CC = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, cfile[i], NULL)) == NULL) die("Can't open ", cfile[i]); xin = CC->header->nx; yin = CC->header->ny; if (xin != xdim || yin != ydim) die("dimension don't match!", cfile[i]); if ((GG = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, gfile[i], NULL)) == NULL) die("Can't ope ", gfile[i]); xin = GG->header->nx; yin = GG->header->ny; if (xin != xdim || yin != ydim) die("dimension don't match!", gfile[i]); if ((CC = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, NULL, cfile[i], CC)) == NULL) die("Can't read ", cfile[i]); if ((GG = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, NULL, gfile[i], GG)) == NULL) die("Can't read ", gfile[i]); corin = CC->data; grdin = GG->data; for (k=0; k<ydim; k++) { for (j=0; j<xdim; j++) { phi[i*xdim*ydim+ydim*j+k] = grdin[j+k*xdim]; if (isnan(grdin[j+k*xdim]) != 0) { flag[j+xdim] = 1; } if (corin[j+k*xdim]>=1 e-2 && corin[j+k*xdim]<=0. 99) { /* Rosen et al. , 2000 IEEE */ var[i*xdim*ydim+ydim*j+k] = sqrt((1. 0 -corin[j+k*xdim]*corin[j+k*xdim])/(corin[j+k*xdim]*corin[j+k*xdim])); } else if (corin[j+k*xdim]<1 e-2) { var[i*xdim*ydim+j*ydim+k] = 99. 99; } else { var[i*xdim*ydim+j*ydim+k] = 0. 1; } } } /* write output */ write_output_ts(API, Out, argc, argv, xdim, ydim, S, flag_rms, flag_dem, disp, vel, res, dem);

Applications of GMT • SRTM 15_PLUS blockmedian, surface grd 2 kml • GMTSAR lon, lat, topo to range azimuth API for read/write netcdf API for 2 -D FFT post processing

filter. csh – various # # make the custom filter 2 and set the decimation # make_gaussian_filter $1 $dec_rng $az_lks $3 > ijdec set filter 2 = gauss_$3 set idec = `cat ijdec | awk -v dc="$dec" '{ print dc*$1 }'` set jdec = `cat ijdec | awk -v dc="$dec" '{ print dc*$2 }'` echo $filter 2 $idec $jdec # # filter the two amplitude images # echo "making amplitudes. . . " conv $az_lks $dec_rng $filter 1 $1 amp 1_tmp. grd=bf conv $idec $jdec $filter 2 amp 1_tmp. grd=bf amp 1. grd rm amp 1_tmp. grd conv $az_lks $dec_rng $filter 1 $2 amp 2_tmp. grd=bf conv $idec $jdec $filter 2 amp 2_tmp. grd=bf amp 2. grd rm amp 2_tmp. grd # # filter the real and imaginary parts of the interferogram # also compute gradients # echo "filtering interferogram. . . " conv $az_lks $dec_rng $filter 1 real. grd=bf real_tmp. grd=bf conv $idec $jdec $filter 2 real_tmp. grd=bf realfilt. grd conv $dec $filter 4 real_tmp. grd xreal. grd conv $dec $filter 5 real_tmp. grd yreal. grd rm real_tmp. grd rm real. grd conv $az_lks $dec_rng $filter 1 imag. grd=bf imag_tmp. grd=bf conv $idec $jdec $filter 2 imag_tmp. grd=bf imagfilt. grd conv $dec $filter 4 imag_tmp. grd ximag. grd conv $dec $filter 5 imag_tmp. grd yimag. grd rm imag_tmp. grd rm imag. grd # # form amplitude image # echo "making amplitude. . . " gmt grdmath realfilt. grd imagfilt. grd HYPOT = amp. grd gmt grdmath amp. grd 0. 5 POW FLIPUD = display_amp. grd # # form the correlation # echo "making correlation. . . " gmt grdmath amp 1. grd amp 2. grd MUL = tmp. grd gmt grdmath tmp. grd $thresh GE 0 NAN = mask. grd gmt grdmath amp. grd tmp. grd SQRT DIV mask. grd MUL FLIPUD = tmp 2. grd=bf conv 1 1 $filter 3 tmp 2. grd=bf corr. grd # # form the phase # echo "making phase. . . " gmt grdmath imagfilt. grd realfilt. grd ATAN 2 mask. grd MUL FLIPUD = phase. grd # # make the Werner/Goldstein filtered phase # echo "filtering phase. . . " phasefilt -imagfilt. grd -realfilt. grd -amp 1. grd -amp 2. grd -psize 16 gmt grdedit filtphase. grd `gmt grdinfo mask. grd -I- --FORMAT_FLOAT_OUT=%. 12 lg` gmt grdmath filtphase. grd mask. grd MUL FLIPUD = phasefilt. grd rm filtphase. grd # # form the phase gradients # echo "making phase gradient. . . " gmt grdmath amp. grd 2. POW = amp_pow. grd gmt grdmath realfilt. grd ximag. grd MUL imagfilt. grd xreal. grd MUL SUB amp_pow. grd DIV mask. grd MUL FLIPUD = xphase. grd gmt grdmath realfilt. grd yimag. grd MUL imagfilt. grd yreal. grd MUL SUB amp_pow. grd DIV mask. grd MUL FLIPUD = yphase. grd

Challenges of GMT Conversion of scripts from GMT 4 to GMT 5. Need simple examples. gmt grdvector GPS_u. nc GPS_v. nc -Ix 15/20 -J $R -O -K -Q 0. 06 i+e+n 0. 03 i -Gblack -W. 6, black -S 1. i --MAP_VECTOR_SHAPE=0. 2 >> map_uv. ps When I kill a job sometimes a gmt keeps running forever. Some scripts are in inches and some in cm. Yes my problem!!!! I don’t understand the gmt. conf file. Need examples.