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

set trace on
set page landscape

# 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


# Clip the image range to a certain depth span of interest.
# And set the colour range:
#  e.g. "0 1 1"    is red at .Zmax. meters
#       ".666 1 1" is blue at .Zmin. meters depth
.Zmin. = -600
.Zmax. =  1000
set image range .Zmin. .Zmax.
set image colorscale hsb 0 1 1 .Zmax. hsb .666 1 1 .Zmin.
set x type linear
set y type linear
set x format "%.0f m"
draw image palette left .Zmin. right .Zmax. increment 200 \
    box ..xmargin.. {rpn ..ymargin.. ..ysize.. + 2 + } \
    {rpn ..xmargin.. ..xsize.. +} {rpn ..ymargin.. ..ysize.. + 3 + }

# 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 -300 -100 unlabelled
draw contour -500 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 OldLace
#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 DarkSeaGreen1
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