(15) Animating 2-D coupled gridding via SVD
Gridding via elastically coupled Green’s splines via gpsgridder is similar to greenspline gridding in that a SVD-based solution allows for the use of a subset of the eigenvalues. While a good first step is to use 25% of them, this animation shows the solutions for all choices of eigenvalues and tracks the reduction of misfit (both total and separately for the east and north components.) The movie shares the same data and setup as one of our test scripts (gpsgridder1.sh) but has been “weaponized” to do it via animation. Note as we include the contributions from the tiniest eigenvalues we rapidly “improve” the misfit while adding spurious variations to the solution. The moral is to not try to fit the data exactly.
#!/usr/bin/env bash
#
# Gridding via elastically coupled Green's splines is similar to greenspline gridding
# in that a SVD-based solution allows for the use of a subset of the eigenvalues.
# While a good first step is to use 25% of them, this animation shows the solutions
# for all choices of eigenvalues and tracks the reduction of misfit (both total and
# separately for the east and north components.) The movie shares the same data and
# setup as one of our test scripts (gpsgridder1.sh) but has been weaponized to do it
# via animation. Note as we include the contributions from the tiniest eigenvalues
# we rapidly "improve" the misfit while adding spurious variations to the solution.
# The moral is to not try to fit the data exactly.
#
# The finished movie is available in our YouTube channel as well:
# https://youtu.be/Pvvc4vb8G4Y
#
# The script took ~20 minutes to render on a 24-core MacPro 2013, with most of the
# time being used to solve for the cumulative and incremental grid components.
# 0. A few common variables to include
cat << 'EOF' > inc.sh
INC=5m # Desired grid spacing
DEC=2 # Let grdvector only plot every other node
R=-R122.5W/115W/32.5N/38N # Area of interest
RATE_CUM=50 # This is in mm/yr or km/Myr - change to use another scalebar
RATE_INC=5 # This is in mm/yr or km/Myr - change to use another scalebar
EOF
# 1. Background script doing the hard work plus making background map
cat << 'EOF' > pre.sh
gmt begin
# Prepare the GPS data set for use with a 1x1 arc minute grid
gmt set MAP_VECTOR_SHAPE 0.5
gmt select @wus_gps_final.txt -R122.5W/115W/32.5N/40N -fg -o0:5 > data.lluv
# Use blockmean to avoid aliasing
gmt blockmean $R -I${INC} data.lluv -fg -i0:2,4 -W > blk.llu
gmt blockmean data.lluv -fg -i0,1,3,5 -W > blk.llv
gmt convert -A blk.llu blk.llv -o0:2,6,3,7 > blk.lluv
gmt select blk.lluv $R -fg | awk '{ print($0," 0 ") }' > data.lluvenct
#
# Do the gridding. There are 2682 data and here we evaluate all possible solutions
gmt gpsgridder $R -I${INC} -Gtmp.grd blk.lluv -fg -Emisfit.txt -Fd8 -C+c+i -S0.5 -V
#
# Mask the grids
cat <<- END > corner.ll
-114.9 37.
-114.9 40.1
-118 40.1
-115.50 37
-115.25 32.5
-114.9 32.5
-114.9 37
END
gmt grdlandmask -Gtmp_mask1.grd -Df
gmt grdmask -Rtmp_mask1.grd corner.ll -Gtmp_mask2.grd -N1/0/0
gmt grdmath tmp_mask1.grd tmp_mask2.grd MUL 0 NAN = mask.grd
# Make a plot of GPS velocity vectors
gmt coast -Rmask.grd -JM4.2i -Glightgray -Ba1f30m -BWSne -Df -X0.4i -Y0.25i -Wfaint
gmt plot @CA_fault_data.txt -W0.5p
gmt velo data.lluvenct -Se.008i/0.95+f8p -A9p -W0.2p,red
# Place the scale using a geovector of length RATE_CUM
echo 121.5W 33N 90 ${RATE_CUM}k | gmt plot -S=0.06i+e+jc -Gblue -W1p,blue
echo 121.5W 33N ${RATE_CUM} mm/yr | gmt text -F+f8p+jCB -D0/0.07i
# Map for incremental vectors
gmt coast -Glightgray -Ba1f30m -BESnw -Df -X4.4i -Wfaint
gmt plot @CA_fault_data.txt -W0.5p
# Place the scale using a geovector of length RATE_INC (10x scaling)
echo 121.5W 33N 90 ${RATE_CUM}k | gmt plot -S=0.06i+e+jc -Gblue -W1p,blue
echo 121.5W 33N ${RATE_INC} mm/yr | gmt text -F+f8p+jCB -D0/0.07i
gmt end
EOF
# 2. Main script for the frames
cat << 'EOF' > main.sh
gmt begin
# Prepare the masked cumulative u,v grids
gmt grdmath tmp_u_cum_${MOVIE_ITEM}.grd mask.grd MUL = GPS_u.grd
gmt grdmath tmp_v_cum_${MOVIE_ITEM}.grd mask.grd MUL = GPS_v.grd
# Plot models
# Shrink down heads of vectors shorter than 10 km
gmt grdvector GPS_u.grd GPS_v.grd -Ix${DEC}/${DEC} -JM4.2i -Q0.06i+e+n10 -Gblue -W0.2p,blue -S100i -X0.4i -Y0.25i
# Prepare the masked incremental grids scaled by 10
gmt grdmath tmp_u_inc_${MOVIE_ITEM}.grd mask.grd MUL 10 MUL = GPS_du.grd
gmt grdmath tmp_v_inc_${MOVIE_ITEM}.grd mask.grd MUL 10 MUL = GPS_dv.grd
gmt grdvector GPS_du.grd GPS_dv.grd -Ix${DEC}/${DEC} -Q0.06i+e+n10 -Gblue -W0.2p,blue -S100i -X4.45i
gmt plot misfit.txt -i0,4 -qi0:${MOVIE_FRAME} -R0/${MOVIE_NFRAMES}/0/7 -JX8.6i/0.7i -X-4.45i -Y4.0i -B -BWSrt -W0.25p,green
# Show decrease in rms misfits
gmt plot misfit.txt -i0,5 -qi0:${MOVIE_FRAME} -W0.25p,blue
gmt plot misfit.txt -i0,3 -qi0:${MOVIE_FRAME} -W1p,red
gmt end
EOF
# 3. Run the movie
gmt movie main.sh -Iinc.sh -Sbpre.sh -CHD -Nanim15 -H8 -M210,png -Tmisfit.txt -Lf+jTC+t"Cumulative and Incremental Vector Contributions for Eigenvalue %4.4d" \
-Ls"rms@-u@-"+f12p,Helvetica,green+jTR+o0.5i/0.4i -Ls"rms@-v@-"+f12p,Helvetica,blue+jTR+o0.5i/0.55i -Ls"rms@-c@-"+f12p,Helvetica,red+jTR+o0.5i/0.7i -Fmp4 -V -Zs