1 / 20

Applications of GMT David Sandwell August 16, 2016

Applications of GMT David Sandwell August 16, 2016. SRTM15_PLUS blockmedian, surface grd2kml GMTSAR lon,lat,topo to range azimuth API for read/write netcdf API for 2-D FFT post processing. Soundings used in SRTM15/30_PLUS multibeam, singlebeam, grids,. Data Sources. multi-beam

sarahall
Download Presentation

Applications of GMT David Sandwell August 16, 2016

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Applications of GMT David Sandwell August 16, 2016 • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

  2. Soundings used in SRTM15/30_PLUS multibeam, singlebeam, grids, . . .

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

  4. 12% in V7

  5. JJ Becker wrote most of this code if ($#argv < 3) then echo " " echo " example: `basename $0` median.xyz surfaceOpts unmasked.grd" echo " " exit endif set bxyz = $1; shift set out = $1; shift set opts = "$*" #/bin/rm -rf *out.grd # set B1 = -180./180./58./90.; set B2 = -180./180./28./62.; set B3 = -180./180./-2./32.; set B4 = -180./180./-32./02.; set B5 = -180./180./-62./-28.; set B6 = -180./180./-90./-58.; # set C1 = -180./180./60./90.; set C2 = -180./180./30./60.; set C3 = -180./180./00./30.; set C4 = -180./180./-30./00.; set C5 = -180./180./-60./-30.; set C6 = -180./180./-90./-60.; # # make the blend file # echo B1.grd -R$C1 1 > blend.txt echo B2.grd -R$C2 1 >> blend.txt echo B3.grd -R$C3 1 >> blend.txt echo B4.grd -R$C4 1 >> blend.txt echo B5.grd -R$C5 1 >> blend.txt echo B6.grd -R$C6 1 >> blend.txt # # do all the subgrids # # 1 # echo $bxyz -V -fg -bid -I15c -A0.5 -R$B1 $opts -GB1$out blockmedian $bxyz -I15c -bid -bod -R$B1 -V > B1.xyz surface B1.xyz -V -fg -bid -I15c -A.50 -Ll-800 -Lu800 -R$B1 $opts -GB1.grd # # # 2 # echo $bxyz -V -fg -bid -I15c -A.707 -R$B2 $opts -GB2$out blockmedian $bxyz -I15c -bid -bod -R$B2 -V > B2.xyz surface B2.xyz -V -fg -bid -I15c -A.707 -Ll-800 -Lu800 -R$B2 $opts -GB2.grd # # 3 # echo $bxyz -V -fg -bid -I15c -A.966 -R$B3 $opts -GB3.grd blockmedian $bxyz -I15c -bid -bod -R$B3 -V > B3.xyz surface B3.xyz -V -fg -bid -I15c -A.966 -Ll-800 -Lu800 -R$B3 $opts -GB3.grd # # 4 # echo $bxyz -V -fg -bid -I15c -A.966 -R$B4 $opts -GB4.grd blockmedian $bxyz -I15c -bid -bod -R$B4 -V > B4.xyz surface B4.xyz -V -fg -bid -I15c -A.966 -Ll-800 -Lu800 -R$B4 $opts -GB4.grd # # 5 # echo $bxyz -V -fg -bid -I15c -A.707 -R$B5 $opts -GB5.grd blockmedian $bxyz -I15c -bid -bod -R$B5 -V > B5.xyz surface B5.xyz -V -fg -bid -I15c -A.707 -Ll-800 -Lu800 -R$B5 $opts -GB5.grd # # 6 # echo $bxyz -V -fg -bid -I15c -A.50 -R$B6 $opts -GB6.grd blockmedian $bxyz -I15c -bid -bod -R$B6 -V > B6.xyz surface B6.xyz -V -fg -bid -I15c -A.50 -Ll-800 -Lu800 -R$B6 $opts -GB6.grd # # now blend all the files together # grdblend blend.txt -G$out -R-180/180/-90/90 -fg -I15c -V

  6. S&S 1 min V18

  7. SRTM30_PLUS V11

  8. SRTM15_PLUS V1

  9. Applications of GMT • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

  10. grd2kml 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: grd2kml.csh grd_file_stem cptfile [-R<west>/<east>/<south>/<north>] " echo " " echo "Example: grd2kml.csh phase phase.cpt " echo " " exit 1 endif # # set DX = `gmt grdinfo $1.grd -C | cut -f8` 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 -Nt0.7 -A60 if ($#argv == 3) then gmt grdimage $1.grd -Igrad.grd -C$2 $3 -Jx1id -P -Y2 -X2 -Q -V > $1.ps else if ($#argv == 2) then gmt grdimage $1.grd -Igrad.grd -C$2 -Jx1id -P -Y2 -X2 -Q -V > $1.ps endif # # now make the kml and png # gmt ps2raster $1.ps -W+k+t"$1" -E$DPI -TG -P -S -V -F$1.png rm $1.ps rm grad.grd rm ps2raster* rm psconvert* #

  11. Applications of GMT • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

  12. GMT code

  13. Applications of GMT • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

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

  15. Applications of GMT • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

  16. # include "gmtsar.h" # include "sbas.h" int main(int argc, char **argv) { /* Begin: Initializing new GMT5 session */ if ((API = GMT_Create_Session (argv[0], 0U, 0U, 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 open ",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]>=1e-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]<1e-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); sbas.c – Xiaopeng Tong

  17. Applications of GMT • SRTM15_PLUS • blockmedian, surface • grd2kml • GMTSAR • lon,lat,topo to range azimuth • API for read/write netcdf • API for 2-D FFT • post processing

  18. filter.csh – various # # 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 amp1.grd amp2.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 = tmp2.grd=bf conv 1 1 $filter3 tmp2.grd=bf corr.grd # # form the phase # echo "making phase..." gmt grdmath imagfilt.grd realfilt.grd ATAN2 mask.grd MUL FLIPUD = phase.grd # # make the Werner/Goldstein filtered phase # echo "filtering phase..." phasefilt -imag imagfilt.grd -real realfilt.grd -amp1 amp1.grd -amp2 amp2.grd -psize 16 gmt grdedit filtphase.grd `gmt grdinfo mask.grd -I- --FORMAT_FLOAT_OUT=%.12lg` 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 # # make the custom filter2 and set the decimation # make_gaussian_filter $1 $dec_rng $az_lks $3 > ijdec set filter2 = 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 $filter2 $idec $jdec # # filter the two amplitude images # echo "making amplitudes..." conv $az_lks $dec_rng $filter1 $1 amp1_tmp.grd=bf conv $idec $jdec $filter2 amp1_tmp.grd=bf amp1.grd rm amp1_tmp.grd conv $az_lks $dec_rng $filter1 $2 amp2_tmp.grd=bf conv $idec $jdec $filter2 amp2_tmp.grd=bf amp2.grd rm amp2_tmp.grd # # filter the real and imaginary parts of the interferogram # also compute gradients # echo "filtering interferogram..." conv $az_lks $dec_rng $filter1 real.grd=bf real_tmp.grd=bf conv $idec $jdec $filter2 real_tmp.grd=bf realfilt.grd conv $dec $dec $filter4 real_tmp.grd xreal.grd conv $dec $dec $filter5 real_tmp.grd yreal.grd rm real_tmp.grd rm real.grd conv $az_lks $dec_rng $filter1 imag.grd=bf imag_tmp.grd=bf conv $idec $jdec $filter2 imag_tmp.grd=bf imagfilt.grd conv $dec $dec $filter4 imag_tmp.grd ximag.grd conv $dec $dec $filter5 imag_tmp.grd yimag.grd rm imag_tmp.grd rm imag.grd

  19. Challenges of GMT Conversion of scripts from GMT4 to GMT5. Need simple examples. gmt grdvector GPS_u.nc GPS_v.nc -Ix15/20 -J $R -O -K -Q0.06i+e+n0.03i -Gblack -W.6,black -S1.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.

More Related