# Example 16 -- Grided bathymetry data, converted to image, then countoured
#
# Contributed 2001-May-14 by Peter S. Galbraith <GalbraithP@dfo-mpo.gc.ca>

`Set image colorscale gebco'
Set up colors for gebco charts of ocean/land altimetry.
{
    open "awk 'BEGIN { \
        for (i = 0; i < 256; i++) { \
            if (i <= 278) {R=0.820; G=0.604; B=0.361;}\
            if (i <= 232) {R=0.882; G=0.988; B=0.969;}\
            if (i <= 227) {R=0.749; G=0.949; B=0.925;}\
            if (i <= 223) {R=0.627; G=0.910; B=0.894;}\
            if (i <= 209) {R=0.514; G=0.871; B=0.871;}\
            if (i <= 185) {R=0.408; G=0.804; B=0.831;}\
            if (i <= 139) {R=0.310; G=0.733; B=0.788;}\
            if (i <=  93) {R=0.220; G=0.655; B=0.749;}\
            if (i <=  46) {R=0.133; G=0.573; B=0.710;}\
            if (i <=   0) {R=0.059; G=0.486; B=0.671;}\
            print(R, G, B); \
        } \
    }' |"
    read image colorscale rgb
    close
}

set trace on
set page landscape

# Set my own variable to emphasize the Gulf only?  Or the Ocean?
.gulf_emphasisGSL. = 1

# Map lat/lon corners
\xl = "73"
\xr = "50"
\yb = "43"
\yt = "53"

# Setup the map parameters
set x format "%.0lf$\circ$W"
set y format "%.0lf$\circ$N"
set x type map W
set y type map N
set x name ""
set y name ""
set y margin 3
set x margin 3
set tic size {rpn ..tic_size.. 2 /} # Make smaller tic sizes (by half)
set y size 14
if {rpn "/usr/local/bin/map-merc-scale-factor" file_exists}
    # Use my proj-based script to calc Mercator resizing if available
    \.factor. = system map-merc-scale-factor -a \yb -A \yt -o \xr -O \xl
    set x size {rpn ..ysize.. \.factor. *}
else
    # else use Gri (Less than 0.2% difference anyway)
    resize x for maps
end if
set x axis \xl \xr  -5 -1
# I want bigger tic marks at degree intervals, even if unlabelled
set font size 0
set y axis \yb \yt 1 {rpn 1 6 /} 
draw y axis 
# Reset to label every second one to draw the labelled axes.
set y axis \yb \yt 2 {rpn 1 6 /} 
set font to Helvetica
set font size 12
draw axes

# Get the x and y grid point locations from the pre-gridded `x y z' data file.
if {rpn "y.grid" file_exists !}
    system zcat topex-gsl.xyz.gz | awk '{print $2}' | uniq > y.grid
    system zcat topex-gsl.xyz.gz | awk '{print $1}' | sort | uniq > x.grid
end if
open x.grid
read grid x
close
open y.grid
read grid y
close

# Generate the grid matrix if not already done
if {rpn "topex-gsl.grid.gz" file_exists !}
    open topex-gsl.xyz.gz
    read columns x y z
    close
    
    set missing value 99999
    # `x y z' data file is regularly-spaced, pre-gridded, 
    # so use neighbour method
    convert columns to grid neighbor
    system rm -f topex-gsl.grid
    write grid to "topex-gsl.grid"
    system gzip --best topex-gsl.grid
else
    open topex-gsl.grid.gz
    read grid data
end if
show grid


if .gulf_emphasisGSL.
    Set image colorscale gebco
    set image range -500 50
    set x type linear
    set y type linear
    set x format "%.0f m"
    draw image palette left -600 right 100 increment 100 \
        box ..xmargin.. {rpn ..ymargin.. ..ysize.. + 2 + } \
        {rpn ..xmargin.. ..xsize.. +} {rpn ..ymargin.. ..ysize.. + 3 + }
else
    Set image colorscale gebco
    set image range -5000 500
    set x type linear
    set y type linear
    set x format "%.0f m"
    draw image palette left -6000 right 1000 increment 1000 \
        box ..xmargin.. {rpn ..ymargin.. ..ysize.. + 2 + } \
        {rpn ..xmargin.. ..xsize.. +} {rpn ..ymargin.. ..ysize.. + 3 + }
end if

# Convert our `x y z' data grid to an image following above specs
convert grid to image size 691 451
set clip postscript on          # Clip in case we go over the edge by a wee bit
draw image
# Draw selected depth contours on the data grid using a thin line
set line width .05
draw contour -50 unlabelled
draw contour -100 -500 -100 unlabelled
draw contour -1000 -3000 -1000 unlabelled

# Draw some fillable coastline data
set missing value -99
open gulfh_lewerIslands.lon_lat.gz
read columns x y
close
read colornames from RGB "/usr/lib/X11/rgb.txt"
set colour rgb 0.820 0.604 0.361
draw curve filled               # Draw filled-in continents in this colour.
set colour black
draw curve                      # Draw again in black, without filling-in.
draw axes frame

open "/data/atlas/gmt/gulff-rivers-all.lon_lat"
read columns x y
set colour rgb {rpn 204 255 /} {rpn 221 255 /} {rpn 225 255 /}
close
draw curve

open gulfh-level2.lon_lat.gz
read columns x y
close
draw curve filled               # Overlay lakes

open gulfh-level3.lon_lat
read columns x y
close
set colour rgb 0.820 0.604 0.361
draw curve filled               # Overlay islands in lakes.

set clip postscript off
set colour black
draw axes frame                 # Draw axes frame again to clean up

quit