Global distribution of antipodes
using Printf
Gwetdry = grdlandmask(region=:global360, inc="30m", res=:crude,
area=500, N="-1/1/1/1/1", reg=true);
# Manipulate so -1 means ocean/ocean antipode, +1 = land/land, and 0 elsewhere
Gkey = gmt("grdmath -fg ? DUP 180 ROTX FLIPUD ADD 2 DIV =", Gwetdry);
# Calculate percentage area of each type of antipode match.
Gscale = gmt("grdmath -Rg -I30m -r Y COSD 60 30 DIV 360 MUL DUP MUL PI DIV DIV 100 MUL =");
Gtmp = gmt("grdmath -fg ? -1 EQ 0 NAN =", Gkey);
Gtmp = Gtmp * Gscale;
key = grd2xyz(Gtmp, skip_NaN=true, flags=:TLf);
ocean = gmt("gmtmath -bi1f -Ca -S ? SUM UPPER RINT =", key);
Gtmp = gmt("grdmath -fg ? 1 EQ 0 NAN =", Gkey)
Gtmp.z = Gtmp.z .* Gscale.z
key = grd2xyz(Gtmp, skip_NaN=true, flags=:TLf)
land = gmt("gmtmath -bi1f -Ca -S ? SUM UPPER RINT =", key)
Gtmp = gmt("grdmath -fg ? 0 EQ 0 NAN", Gkey)
Gtmp = Gtmp * Gscale;
key = grd2xyz(Gtmp, skip_NaN=true, flags=:TLf)
mixed = gmt("gmtmath -bi1f -Ca -S ? SUM UPPER RINT =", key)
# Generate corresponding color table
C = makecpt(color="blue,gray,red", range=(-1.5,1.5,1))
# But unfortunately this palette is not correct, so lets patch it
C.colormap[1,1:2] .= 0; C.colormap[1,3] = 1;
C.colormap[3,2:3] .= 0; C.colormap[3,1] = 1;
# Create the final plot and overlay coastlines
grdimage(Gkey, proj=(name=:EckertVI, center=180), figsize=20,
frame=(axes=:WsNE, title="Antipodal comparisons"), xaxis=(annot=60,),
yaxis=(annot=30,), yshift=3, interpol=:n, portrait=false)
coast!(shore=:thinnest, area=500)
# Place an explanatory legend below
legend!(pos=(anchor=:BC, width=15), box=(pen=:thick,), yshift=-0.5,
par=(:FONT_ANNOT_PRIMARY,10),
text_record(["N 3"
@sprintf("S 0.15i s 0.2i red 0.25p 0.3i Terrestrial Antipodes [%d %%]", land[1].data[1])
@sprintf("S 0.15i s 0.2i blue 0.25p 0.3i Oceanic Antipodes [%d %%]", ocean[1].data[1])
@sprintf("S 0.15i s 0.2i gray 0.25p 0.3i Mixed Antipodes [%d %%]", mixed[1].data[1])]), show=1)
See also GMT ex25