#!/bin/bash
#
# make a plot of GPS velocity vectors 
#
rm -f map_uv.ps
#
gmt gmtset PS_MEDIA 620x500
makecpt -Cgray -T-2.8/1.4/.1 -Z > grad.cpt
gmt makecpt -Crainbow -T-4/44/1 -D -M -Z > v_tot.cpt
D=../02_models/gpsgridder/data.lluvenct
R=-R-123.0/-114/32.0/37.5
gmt grdcut topo.grd $R -Gtmp.grd
gmt grdgradient tmp.grd -A0 -Nt.4 -Ggrad.grd -V
gmt grdsample ve_mean.grd -Rgrad.grd -Gve.nc
gmt grdsample vn_mean.grd -Rgrad.grd -Gvn.nc
U=ve.nc
V=vn.nc
gmtset  FORMAT_GEO_MAP = dddF
gmt select $D $R -fg | awk '{ print($0," 0 ") }' > data.lluvenct
#
#   first make a mask
#
gmt grdlandmask -Gmask.grd -R$U -Df -NNaN/1
gmt grdmath $U mask.grd MUL = GPS_u.nc
gmt grdmath $V mask.grd MUL = GPS_v.nc
gmt grdmath GPS_u.nc GPS_v.nc HYPOT = GPS_tot.nc
#
#   make the map
#
gmt pscoast $R -JM7i -P -Glightgray -Ba1f30m -BWSne -K -Df -X0.7i -Wfaint > map_uv.ps
gmt grdimage grad.grd -Cgrad.cpt -R -J -O -K >> map_uv.ps
gmt grdimage GPS_tot.nc -Igrad.grd -Cv_tot.cpt -R -J -Q -O -K >> map_uv.ps
gmt psxy fault_data.gmt -J -R -W.5 -O -V -K >> map_uv.ps
#gmt psvelo data.lluvenct -J $R -Se.008i/0.95/8 -A9p -W.2,red -O -K >> map_uv.ps
gmt grdvector GPS_u.nc GPS_v.nc -Ix20/30 -J $R -O -K -Q0.06i+e+n0.03i -Gblack -W.6,black -S1.i --MAP_VECTOR_SHAPE=0.2 >> map_uv.ps
#gmt psvelo $D -J $R -Se.008i/0.95/8 -A9p -W.2,red -O -K >> map_uv.ps
psscale -Cv_tot.cpt -D5/2.0/8/0.4h -B5 -O -K >> map_uv.ps
pstext -R -J -O <<!>> map_uv.ps
-121.2 32.2 12 0 5 1 velocity (mm/yr)
!
rm *.nc