#!/bin/bash
set -x #Uncomment to echo these commands.
#######################################################
#Clarifying definitions. Do not change!################
start=" -K " #Should always redirect using > to write new PS.
middle=" -O -K " #Should always redirect using >> to append to PS.
end=" -O " #Should always redirect using >> to append to PS.
#######################################################
gmt5=0 #1/0 = GMT v5/4. GMTv5 support not finished.
if [ $gmt5 != 0 ]
then
  gmt_prefix="gmt "
else
  gmt_prefix=""
fi
color_scheme=1 #1/2=white/black background
montage=0 #1=left-justify titles, add (a),(b), run montage.sh.
prefixes=('(a)' '(b)' '(c)' '(d)' '(e)' '(f)' '(g)' '(h)' '(i)' '(j)' '(k)' '(l)' '(m)' '(n)' '(o)' '(p)' '(q)' '(r)' '(s)' '(t)' '(u)' '(v)' '(w)' '(x)' '(y)' '(z)')
index=-1 #Increments on each map, accesses prefixes above for montage.
png_options=" -P -Tg " #PDF default: -E720, else 300 dpi.
if [ $montage != 0 ]
then
  png_options=" -A"$png_options
fi
#Force off-white(dark gray) fore(back)ground color because
#flip_backgrounds.sh can change the maps' text from
#black to white, and their backgrounds from white to black.
$gmt_prefix gmtset COLOR_BACKGROUND=2/2/2 COLOR_FOREGROUND=253/253/253
digits=4
$gmt_prefix gmtset D_FORMAT=%.${digits}f
#######################################################
data_name="latlon_data_0001.nc"
plot_base="map_parameter_0001"
let index=$index+1
#######################################################
#Set resolution, coast_file, coast_thickness, and coastlines
#on first map only because they should be universal.
coast_file="/home/bryan/GAIA/data/ancillary/Rignot/InSAR_GL_Antarctica.txt"
resolution=" -E50 " #50/2000 is low/high quality.
pscoast_res=" -Di+ "
pscoast_res_orig=$pscoast_res #Don't want USA maps to repeatedly add -N2.
pscoast_thk="0.6"
coast_thk="0.009"
coastlines=1 #1:pscoast, 2:pscoast+InSAR.
#coast_color is gray82 for off-white, or gray10 for dark coastlines.
title_format="0 0 30 0 0 MC"
blurb_format="0 0 15 0 1 ML"
units_format="0 0 13 0 0 MC"
scale_units="meters"
misc_range=" -R0/1/0/1 -JX1c "
#grdcut requires actual limits, but if grdimage uses them: GMT Fatal Error: grdimage could not allocate memory [21.69 Gb, n_items = 5823567396]
minlon=0.000
maxlon=360.000
minlat=-90.000
maxlat=90.000
#Global projections:
#    1 - Robinson
#    2 - Winkel Tripel
#    3 - Mollweide
#    4 - Miller
#Polar projections:
#  101 - N. Azimuthal Equidistant
#  102 - S. Azimuthal Equidistant
#Specific regions:
# 1001 - North America
# 1002 - South America
# 1003 - Africa
# 1004 - Greenland
# 1005 - South Asia
# 1006 - Australia
# 1007 - Europe
# 1101 - Contiguous United States
# 1102 - California
projection_choice=2
standard_circle=0 #1=all specific regions use standard circular projection.
standard_rect=0 #1=all specific regions use standard rectangular projection.
. ./projections.sh
if [ $projection_choice == 101 ]
then
  minlat=0.000
  actual_range=" -R0.0/360.0/$minlat/90.0 "
  polar_radius=$(bc <<< "scale=5; 90-$minlat")
  projection=" -JE0/90.0/$polar_radius/${map_width}c " #N. Azimuthal Equidistant
elif [ $projection_choice == 102 ]
then
  maxlat=0.000
  actual_range=" -R0.0/360.0/-90.0/$maxlat "
  polar_radius=$(bc <<< "scale=5; 90+$maxlat")
  projection=" -JE0/-90.0/$polar_radius/${map_width}c " #S. Azimuthal Equidistant
fi
range=" -R$minlon/$maxlon/$minlat/$maxlat "
map_pos=" -Xa${map_x}c -Ya${map_y}c "
scale_pos=" -D${scale_x}c/${scale_y}c/${scale_length}c/${scale_width}c "
units_x=$(bc <<< "scale=5; $scale_x+$scale_width/2")
units_y=$(bc <<< "scale=5; $scale_y+$scale_length/2+0.8")
units_pos=" -Xa${units_x}c -Ya${units_y}c "
blurb_pos=" -Xa${blurb_x}c -Ya${blurbs_y}c "
blurb2_pos=" -Xa${blurb2_x}c -Ya${blurbs_y}c "
subset_name=$data_name"_subset" #Weird to append, but works with ../pl..
$gmt_prefix grdcut $data_name -G$subset_name $actual_range #Only use subset values for scale.
$gmt_prefix grdinfo -L0 $subset_name > subset_grdinfo #Extract data range from subset.
data_min_e=$(awk '/z_min: /{printf "%.'$digits'e\n", $3}' subset_grdinfo)
data_max_e=$(awk '/z_max: /{printf "%.'$digits'e\n", $5}' subset_grdinfo)
data_min_f=$(awk '/z_min: /{printf "%.'$digits'f\n", $3}' subset_grdinfo)
data_max_f=$(awk '/z_max: /{printf "%.'$digits'f\n", $5}' subset_grdinfo)
notation="f" #By default, numbers appear in floating point format.
data_min_print=$data_min_f #By default, numbers appear in floating point format.
data_max_print=$data_max_f #By default, numbers appear in floating point format.
palette=" -Cwysiwyg " #0-360phase=cyclic,amps=wysiwyg,trend=polar,topo=relief,pts=rainbow
symmetric_limit="-1.000000"
if [[ $(bc <<< "$symmetric_limit <= 0.0") == 1 || $(bc <<< "$data_min_f >= 0.0") == 1 || $(bc <<< "$data_max_f <= 0.0 ") == 1 ]]
then
  upper_limit=$data_max_f
  lower_limit=$data_min_f
else
  upper_limit=$symmetric_limit
  lower_limit=-$upper_limit
fi
#upper_limit=0.0
#lower_limit=0.0
limits=" -L$lower_limit/$upper_limit "
#limits="" #Default limits include all values (not just in subset), aren't always symmetric.
#limits=" -T= " #Symmetric limits that include all values (not just in subset).
. ./overflow.sh
if [[ $(bc <<< "$data_min_f >= 0.0") == 1 ]]
then
  #If all values are non-negative, they might cluster near bottom of
  #scale which is dark so use off-white coastlines.
  coast_color="gray82"
else
  #Otherwise, they might cluster near the scale center which is light so use dark coastlines.
  coast_color="gray10"
fi
numlevels="10"
$gmt_prefix grd2cpt $data_name $palette $limits -E$numlevels -Z > map.cpt
. ./notation.sh
blurb_contents="Data range: [$data_min_print, $data_max_print] $scale_units"
title="Elevation and bathymetry"
$gmt_prefix grdimage $data_name $boundary $resolution $range $projection $map_pos -Cmap.cpt $start > $plot_base.ps
$gmt_prefix pscoast -W$pscoast_thk/$coast_color $pscoast_res $range $projection $map_pos $middle >> $plot_base.ps
if [ $coastlines == 2 ]
then
  $gmt_prefix psxy -N $coast_file -: -Sc$coast_thk -W$coast_thk/$coast_color $range $projection $map_pos $middle >> $plot_base.ps
fi
#Uncomment to put a marker at echoed coords, given as lon lat:
#echo -85.19 -77.36 | $gmt_prefix psxy -N -S+0.5c -W5/244/164/96 -G244/164/96 $range $projection $map_pos $middle >> $plot_base.ps
$gmt_prefix psxy $data_name -bcmascon_lons/mascon_lats -Sc0.01c -G139/69/19 $range $projection $map_pos $middle >> $plot_base.ps
if [ $montage != 0 ]
then
  title=${prefixes[$index]}" "$title
  title_format="0 0 30 0 0 ML" #Left-justify so montage titles are uniform.
  title_x=$(bc <<< "scale=5; $blurb_x-0.1")
else
  title_x=$(bc <<< "scale=5; $map_x+$map_width/2")
fi
title_pos=" -Xa${title_x}c -Ya${title_y}c "
echo $title_format $title | $gmt_prefix pstext -N $title_pos $misc_range $middle >> $plot_base.ps
cpt_name="-Cmap.cpt "
$gmt_prefix psscale $cpt_name -L $scale_format $overflow $scale_pos -A $middle >> $plot_base.ps
#Print units manually, otherwise they're too close to numbers on scale.
echo $units_format $scale_units | $gmt_prefix pstext -N $units_pos $misc_range $middle >> $plot_base.ps
echo $blurb_format $blurb_contents | $gmt_prefix pstext -N $blurb_pos $misc_range $middle >> $plot_base.ps
blurb2_contents="" #Error bar: N/A $scale_units"
echo $blurb_format $blurb2_contents | $gmt_prefix pstext -N $blurb2_pos $misc_range $end >> $plot_base.ps
ps2raster $png_options $plot_base.ps #Convert PS to PNG format.
#ps2raster -P -Tf $plot_base.ps #Convert PS to PDF, if uncommented.
rm -f $plot_base.ps
mv map.cpt Zbackup_cpt_$plot_base.cpt
previous_data_name=$data_name
#######################################################
. ./trim.sh
#######################################################
if [ $color_scheme == 2 ]
then
  . ./flip_backgrounds.sh
fi
#######################################################
if [ $montage != 0 ]
then
  . ./montage.sh
fi