Revisiting GMT, perl, and proj

The time has come to make a site map for our ship source level paper.  I always love this stage of a project — map making — because it still takes me back to some of the geekiest open source Unix stuff my PhD advisor, Russ McDuff, taught me: using public bathymetric data in Generic Mapping Tools (GMT) and Proj within Perl to make awesome maps of the ocean for free!  It’s always a chore to get all these tools working within your Li/Un/*ix operating system, but below is an example of what they can do from this working directory



In my latest attempt (within a Dreamhost.com account), I had to remember (again) the syntax and commands for setting path variables…

  1. edit .bashrc to include a path to your proj and GMT binaries
  2. export PATH=~/bin:~/usr/local:~/usr/local/bin:~/beamreach.org/maps/gmt/bin:$PATH
  3. execute: . .bashrc to reinitialize the PATH variable
  4. check PATH with env
  5. check all is well with which proj

Here’s the perl script for making the (above) rough draft map of Haro Strait and the San Juan Islands.

#!/usr/bin/perl

#`gmtset FRAME_PEN 1p BASEMAP_TYPE PLAIN DEGREE_FORMAT 2`;
#`gmtset X_AXIS_LENGTH 8.5i X_AXIS_LENGTH 11i`;
#`gmtset ANOT_FONT_SIZE 8 ANOT_OFFSET 0.037`;
#$gxscl=4; $gyscl=4;
#`gmtset GLOBAL_X_SCALE $gxscl  GLOBAL_Y_SCALE $gyscl`;

# Haro Strait region
#$prj=2;
#$prj=6;        # inches/degree
#$prj=13;
$prj=26;
#$prj=39;
#$prj=50;
$proj="-Jm$prj";
$lonmin=-123.5; $lonmax=-123; $latmin=48.25; $latmax=48.75;
$bound="-R$lonmin/$lonmax/$latmin/$latmax";
$tick="-Ba0.5f0.25/a0.5f0.25EwsN";
#$xoff=2.5;  $yoff=4;
$xoff=2;  $yoff=2;
#$xoff=0;  $yoff=0;

$out_nm="ship-sl-map";
$out_ps="$out_nm.ps";

#$basemap="/home/scottveirs/beamreach.org/maps/topography/etopo2-salish.grd";
$basemap="/home/scottveirs/beamreach.org/maps/topography/salishc_4938_data/salishc_4938/salishsea-noaa.grd";
#$cptfile="/home/scottveirs/beamreach.org/maps/scripts/cpts/hydronet-1000-3000-100m.cpt";
#$cptfile="/home/scottveirs/beamreach.org/maps/scripts/cpts/etopo2.cpt";
$cptfile="/home/scottveirs/beamreach.org/maps/scripts/cpts/terrain.cpt";
#$cptfile="/home/scottveirs/beamreach.org/maps/scripts/cpts/wiki-djerba.cpt";

##    PLOT BATHYMETRY
`psbasemap $proj $bound -K -P -V -X$xoff -Y$yoff $tick > $out_ps`;
`grdimage $basemap -C$cptfile $proj $bound -P -O -K -V >> $out_ps`;
`grdcontour $basemap -C50 $proj $bound -W.05p -Q10 -P -O -K -V >> $out_ps`;
#`grdcontour $basemap -Cfig1cnt $proj $bound -W1p/255/255/255 -Q10 -P -O -K -V >> $out_ps`;
#`grdcontour $basemap -C10 $proj $bound -W.05p -Q10 -P -O -K -V >> $out_ps`;
#`grdcontour $basemap -C20 -A100f3 $proj $bound -W1p -Q10 -P -O -K -V >> $out_ps`;
$cntpen="0.2p/100/100/100";
#`grdcontour $basemap -C500 $proj $bound -W$cntpen -Q10 -P -O -K -V >> $out_ps`;
#`grdcontour $basemap -C0 $proj $bound -W$cntpen -Q10 -P -O -K -V >> $out_ps`;
#`grdcontour $basemap -C2000 -Af3 $proj $bound -W0.5p/50/50/50 -Q10 -P -O -K -V >> $out_ps`;

#`pscoast $proj $bound -Dh -A9 -W -G255/255/255 -S240/240/240 -N1 -N2 -O -K -V >> $out_ps`;
#`pscoast $proj $bound -Dh -A9 -W -N1 -N2 -O -K -V >> $out_ps`;
# Land grey:
#`pscoast $proj $bound -Dh -A20 -W0.5p/10/10/10 -G200/200/200 -O -K -V >> $out_ps`;
`pscoast $proj $bound -Dh -A20 -W0.5p/10/10/10 -O -K -V >> $out_ps`;

# Label pertinent features...
open (WORDS, ">labels");
$wfont=22;
$nodefont=24;
print (WORDS "-123.75  48.445    $wfont 0 0 2 Victoria\n");
#print (WORDS "-124.690  48.245    $nodefont 0 0 TL Neah Bay\n");
#print (WORDS "-122.750  48.060    $nodefont 0 0 TR Marine Science Center\n");
#print (WORDS "-123.300  48.550    $nodefont 0 0 BR Orcasound lab\n");
print (WORDS "-123.065  48.525    $nodefont 0 0 TL Lime Kiln\n");
#print (WORDS "-122.700  48.425    $nodefont 0 0 TL State Park\n");
print (WORDS "-122.700  47.635    $nodefont 0 0 TC Seattle\n");
close(WORDS);
#$xtoff="0.1";
#$xyoff="0.1";
#`pstext -Jm -R -G0/0/0 -Dj$xtoff/$ytoff-O -K -V < labels >> $out_ps`;
# No background rectangles
#`pstext -Jm -R -G0/0/0 -O -K -V < labels >> $out_ps`;
# Background rectangles with rounded corners: [ -W[color,][o|O|c|C[pen]]
#`pstext -Jm -R -G0/0/0 -W255/255/255,Othick,255/0/0 -O -K -V < labels >> $out_ps`;
`rm labels`;

# Put symbols on nodes
open(SYM, ">nodes.lonlat") || die "Can't open nodes.lonlat!\n";
# print SYM "-124.675  48.390\n";  # Neah Bay
# print SYM "-122.750  48.125\n"; # Port Townsend
# print SYM "-123.200  48.550\n"; # Orcasound
 print SYM "-123.150  48.500\n"; # Lime Kiln
# print SYM "-122.350  47.600\n"; # Seattle
close(SYM);
`psxy nodes.lonlat -Jm -R -P -Sc0.5/0/0/0 -G255/0/0 -V -O -K >> $out_ps`;

#Put more symbols on...
`psxy catscradle.lonlat -Jm -R -P -Sc0.05/0/0/0 -G255/0/0 -V -O -K >> $out_ps`;



##    INSET LOCATION MAP WITH LABELS
`gmtset FRAME_PEN 4p BASEMAP_TYPE PLAIN`;
# Good offset if using WS lat/lon labels
#$locyoff=1.8; $locxoff=2.8;
$locyoff=0.75; $locxoff=0.75;
#`psbasemap -Jm0.2 -R-135/-105/35/55 -O -K -X$locxoff -Y$locyoff -P -Ba10f5/a5f5WeSn >> $out_ps`;
#`psbasemap -Jm0.2 -R-135/-105/35/55 -O -K -X$locxoff -Y$locyoff -P -Ba10f5/a5f5wesn >> $out_ps`;
`psbasemap -Jm2.5 -R-125/-122/47/48.5 -O -K -X$locxoff -Y$locyoff -P -Ba10f5/a5f5wesn >> $out_ps`;
# plot the coast line
#`pscoast -Jm -R -O -K -W -G210/180/140 -S135/206/250 -N1 -N2 >> $out_ps`;
#`pscoast -Jm -R -Di -A200 -O -K -W -G216/242/254 -S148/191/139 -N1 -N2 -V >> $out_ps`;
# Seattle Aquarium Kiosk colors
#`pscoast -Jm -R -Di -A200 -O -K -W -G172/208/165 -S186/232/254 -N1 -N2 -V >> $out_ps`;
`pscoast -Jm -R -Di -A200 -O -K -W -G200/200/200 -S255/255/255 -N1 -N2 -V >> $out_ps`;
`gmtset FRAME_PEN 1p BASEMAP_TYPE PLAIN`;

# Put symbols on key cities
open(SYM, ">cities.lonlat") || die "Can't open cities.lonlat!\n";
# print SYM "-123.45 48.45\n";
close(SYM);
`psxy cities.lonlat -Jm -R -P -Sc0.5/0/0/0 -G0/0/0 -V -O -K >> $out_ps`;

# Small box around the zoomed area
$lonboxmax=$lonmax; $lonboxmin=$lonmin;
$latboxmax=$latmax; $latboxmin=$latmin;
$marg="3";
open(BOX, ">magbox.lonlat") || die "Can't open magbox.lonlat!\n";
print BOX "$lonboxmax-$marg   $latboxmax+$marg\n";
print BOX "$lonboxmax-$marg      $latboxmin-$marg\n";
print BOX "$lonboxmin+$marg      $latboxmin-$marg\n";
print BOX "$lonboxmin+$marg      $latboxmax+$marg\n";
print BOX "$lonboxmax-$marg      $latboxmax+$marg\n";
close(BOX);
`psxy magbox.lonlat -Jm -R -P -W4p/0/0/0 -V -O >> $out_ps`;
`rm magbox.lonlat`;

`convert $out_ps $out_nm.png`;