# Draw map of Maritime provinces of Canada, to demonstrate
# 'set map projection' command. (The data file used here
# is available upon request.)

# Use lambert-conformal-conic (lcc) projection, with the Clark 1966
# formula describing the ellipsoidal shape of the earth. The projection
# is centered at a longitude of 63.5 west (near Halifax, Nova Scotia).

# This projection requires intersection latitudes (circles at which
# the cone touches the earth surface). Here I use 43N and 49N to
# match Canadian Hydographic Service charts of the region.
#
# NB. The -m argument sets a scale factor of 0.001, so that the unit
# will be kilometers rather than the default meters (see axis commands
# below).

set map projection "+proj=lcc +ellps=clrk66 +lon_0=63.5w +lat_0=44.6n +lat_1=43 +lat_2=49 -m 1/1000";

# Set up axes to view rectangle 500 kilometers on a side.
#
# NB. It is crucial that the x and y ranges be identical or the
# map will be distored.

set x axis -300 400 100;
set y axis -250 350 100;
draw axes none; # Don't draw axes

# Set clip on, so that coastline doesn't plot outside the
# box containing the map.
set clip on;

# The data file is in a weird format, with header lines that
# have the letter "P" in them. Compensate by passing
# data though a gawk script first.
set missing value -999;
open(IN, "gawk '{if($1 == \"P\"){print -999,-999}else{print}}' /users/dek/kelley/data/Coastline/NS_1.dat |") or die;

read columns IN x=1 y=2;
close(IN);

# Draw a black coastline surrounding an 80-percent
# gray land-mass.
set graylevel 0.8;
draw curve filled;
set color black;
draw curve;

# Draw frame around box
set clip off;
draw axes frame;

# Draw graticule crosses, with labels
for ($lat = 43; $lat < 48; $lat += 1) {
        for ($lon = -67; $lon < -58; $lon += 2) {
                draw line from ($lon - 1/10) $lat to ($lon + 1/10) $lat;
                draw line from $lon ($lat - 1/10) to $lon ($lat + 1/10);
        }
}
set font size 9;
for ($lon = -67; $lon < -58; $lon += 2) {
        draw label (-$lon . "W") whiteunder at ($lon - 0.2) 42.7;
}
for ($lat = 43; $lat < 48; $lat += 2) {
        draw label ($lat . "N") whiteunder at -58.8 ($lat - 0.04);
}