GMT scripts for Enderby Trip visualisation

antipodes.sh

#!/bin/bash
#		This map is a quick and dirty adaptation of GMT EXAMPLE 25
#   (using GMT 5.4.0_r17345)
# Purpose:	Display distribution of antipode types around NZ
# GMT modules:	gmtset, grdlandmask, grdmath, grd2xyz, gmtmath, grdimage, pscoast, pslegend
# Unix progs:	cat
#
# Create D minutes global grid with -1 over oceans and +1 over land
D=10
gmt grdlandmask -Rg -I${D}m -Dc -A500 -N-1/1/1/1/1 -r -Gwetdry.nc
# Manipulate so -1 means ocean/ocean antipode, +1 = land/land, and 0 elsewhere
gmt grdmath -fg wetdry.nc DUP 180 ROTX FLIPUD ADD 2 DIV = key.nc
# Calculate percentage area of each type of antipode match.
gmt grdmath -Rg -I${D}m -r Y COSD 60 $D DIV 360 MUL DUP MUL PI DIV DIV 100 MUL = scale.nc
gmt grdmath -fg key.nc -1 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
ocean=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
gmt grdmath -fg key.nc 1 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
land=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
gmt grdmath -fg key.nc 0 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
mixed=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
# Generate corresponding color table
gmt makecpt -Cblue,gray,red -T-1.5/1.5/1 -N > key.cpt
# Create the final plot and overlay coastlines
gmt set FONT_ANNOT_PRIMARY +10p FORMAT_GEO_MAP dddF
range="g"
proj="G-185/-45/4.5i"
psfile='enderby_antipodes.ps'
#plot antipodes raster
gmt grdimage key.nc -R${range} -J${proj} -Ckey.cpt -B10/10 -P -K > $psfile
#overplot landcover
gmt pscoast -R${range} -J${proj} -Dl -A0/0/1 -G35 -W0.4p -O -K >> $psfile
#add legend
gmt pslegend -R${range} -J${proj} -O -DJBC+w4.5i -Y-0.75i >> $psfile << END
N 2
S 0.15i s 0.2i gray  0.25p 0.3i Terrestrial Antipodes
S 0.15i s 0.2i blue 0.25p 0.3i Oceanic Antipodes 
END
#convert to PNG
gmt psconvert -Tg -A $psfile
#remove intermediate results
rm -f wetdry.nc scale.nc key.nc tmp.nc key.cpt key.b gmt.conf

enderby.sh

#map of Enderby Route
#
range="155/190/-55/-33"
proj="N12c"
psfile='enderby-overview.ps'
pngfile='enderby-overview.png'
GMT grdimage etopo5.grd -R${range} -J${proj} -Cibcao.cpt -B20/10 -P -K > $psfile
#plot route; need to replace invercargill coords with bluff coords and add waypoints
(echo "168.3475 -46.413056"; echo "166.533333 -48.016667"; echo "166.1 -50.7"; echo "169.145 -52.54"; echo "168.3475 -46.413056") | GMT psxy -R${range} -J${proj} -W2p,red+s -O -K >> $psfile
#GMT psbasemap -R${range} -J${proj}  -Gwhite -O -K >> $psfile
GMT pscoast -R${range} -J${proj} -Di -A0/0/1 -Glightgray -W0.4p -O  >> $psfile
#plot currents from talley
#GMT psxy 'data3.txt' -J${proj} -R${range} -Sc0.05 -Gred -O -K >>  $psfile
#GMT psxy 'atldata1.txt' -J${proj} -R${range} -Sc0.05 -Gyellow -O -K >>  $psfile
#GMT psxy 'atldata2.txt' -J${proj} -R${range} -Sc0.05 -Ggreen -O -K >>  $psfile
GMT psconvert -A -Tf $psfile
GMT psconvert -A -Qg1 -TG $psfile
#imagemagick overlay
#composite -gravity center $overlay $pngfile IOcurrents_globe.png
#composite -gravity center $overlay2 IOcurrents_globe.png IOcurrents_globe_shaded.png

#map of way to nz
#
range="g"
proj="G-135/-15/4.5i"
psfile='enderby_globe.ps'
pngfile='enderby_globe.png'
overlay='enderby_globe_overlay.png'
overlay2='Orthographic_gradient_1356x1356px.png'
GMT grdimage etopo5.grd -R${range} -J${proj} -Cibcao.cpt -B10/10 -P -K > $psfile
#GMT psbasemap -R${range} -J${proj}  -Gwhite -O -K >> $psfile
GMT pscoast -R${range} -J${proj} -Dl -A0/0/1 -Glightgray -W0.4p -O >> $psfile
#plot currents from talley
#GMT psxy 'data3.txt' -J${proj} -R${range} -Sc0.05 -Gred -O -K >>  $psfile
#GMT psxy 'atldata1.txt' -J${proj} -R${range} -Sc0.05 -Gyellow -O -K >>  $psfile
#GMT psxy 'atldata2.txt' -J${proj} -R${range} -Sc0.05 -Ggreen -O -K >>  $psfile
#plot chagos location
#echo "72 -6" | GMT psxy -R${range} -J${proj} -Ba0f0g0 -Sa0.25 -A5000 -Gblack -O -K >> $psfile
#plot seamounts
#GMT psxy 'airports.dat' -R${range} -J${proj} -W1p,red+s -O -K >> $psfile
#(echo "168.3475 -46.413056"; echo "166.533333 -48.016667"; echo "166.1 -50.7"; echo "169.145 -52.54"; echo "168.3475 -46.413056") | GMT psxy -R${range} -J${proj} -W1p,yellow+s -O -K >> $psfile
#(echo "168.3475 -46.413056"; echo "170.5 -45.866667"; echo "167.716667 -45.416667"; echo "168.3475 -46.413056") | GMT psxy -R${range} -J${proj} -W1p,green+s -O -K >> $psfile
#GMT psxy 'airports.dat' -J${proj} -R${range} -Sc0.08 -Gred -O >>  $psfile
GMT psconvert -A -Tf $psfile
GMT psconvert -A -Qg4 -TG $psfile
#imagemagick overlay
#composite -gravity center $overlay $pngfile IOcurrents_globe.png
composite -gravity center $overlay2 $pngfile enderby_globe_shaded.png
composite -gravity center $overlay enderby_globe_shaded.png enderby_globe_combo.png