```# 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); \
} \
}' |"
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
close
open y.grid
close

# Generate the grid matrix if not already done
if {rpn "topex-gsl.grid.gz" file_exists !}
open topex-gsl.xyz.gz
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
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
close
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"
set colour rgb {rpn 204 255 /} {rpn 221 255 /} {rpn 225 255 /}
close
draw curve

open gulfh-level2.lon_lat.gz
close
draw curve filled               # Overlay lakes

open gulfh-level3.lon_lat