#!/bin/bash

# edit clean.dat and remove the bad point in San Joaquin Valley by hand
#awk '$1>-119.4 && $1<-119.2 && $2>35.61 && $2<35.63 {print $0}' junk2.dat
#-119.349       35.623  0.8 0.34 -0.83 0.30     5.41244009103   -5.4890783

rm junk*
rm GPS_resid_east.dat GPS_resid_north.dat
awk '{print $1,$2,$3,$5,$4,$6 }' < gps.dat > junk0.dat
#awk '{print $1,$2,$3+15.7,$5,$4-17.3,$6 }' < clean.dat > junk0.dat
#
#   change the model to pixel registration
#
grdsample ve0.grd -Gve.grd -r
grdsample vn0.grd -Gvn.grd -r
#
#  extract the model from the grid
#
gmt grdtrack junk0.dat -Gve.grd > junk1.dat
gmt grdtrack junk1.dat -Gvn.grd > junk2.dat
#
# gridding the clean.dat
#
awk '{print ($1,$2,$3-$7,1./$4)}' < junk2.dat | gmt blockmedian -R-125/-113.333/32/40.336 -I.4 -r -Wi | awk '{print ($1, $2, $3)}'  >  GPS_resid_east.dat
awk '{print ($1,$2,$5-$8,1./$6)}' < junk2.dat | gmt blockmedian -R-125/-113.333/32/40.336 -I.4 -r -Wi | awk '{print ($1, $2, $3)}'  >  GPS_resid_north.dat

#
#  now grid the residuals and add back to the velocity
#
#gmt surface GPS_resid_east.dat -Gve_resid_surface.grd -T0.9 `gmt grdinfo -I ve.grd` `gmt grdinfo -I- ve.grd` -V
#gmt surface GPS_resid_north.dat -Gvn_resid_surface.grd -T0.9 `gmt grdinfo -I vn.grd` `gmt grdinfo -I- vn.grd` -V

#  perhaps a bit too smooth axi  mean, std, -6.01316 15.984
gmt greenspline GPS_resid_east.dat -Gve_resid_surface.grd -D2 -Cn100 `gmt grdinfo -I ve.grd` `gmt grdinfo -I- ve.grd` -r -V
gmt greenspline GPS_resid_north.dat -Gvn_resid_surface.grd -D2 -Cn100 `gmt grdinfo -I vn.grd` `gmt grdinfo -I- vn.grd` -r -V

#gmt greenspline GPS_resid_east.dat -Gve_resid_surface.grd -D2 -Cv99.99 `gmt grdinfo -I ve.grd` `gmt grdinfo -I- ve.grd` -V
#gmt greenspline GPS_resid_north.dat -Gvn_resid_surface.grd -D2 -Cv99.99 `gmt grdinfo -I vn.grd` `gmt grdinfo -I- vn.grd` -V
#
gmt grdmath ve.grd ve_resid_surface.grd ADD = ve_tot.grd
gmt grdmath vn.grd vn_resid_surface.grd ADD = vn_tot.grd

./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 vn_resid_surface 1
./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 ve_resid_surface 1

./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 vn_tot 5
./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 ve_tot 5

#
rm GPS_resid_east.dat GPS_resid_north.dat
rm junk*
gmt ps2raster ve_tot_ll.ps -A -Tj
gmt ps2raster vn_tot_ll.ps -A -Tj
gmt ps2raster ve_resid_surface_ll.ps -A -Tj
gmt ps2raster vn_resid_surface_ll.ps -A -Tj