How can I plot IMS-24km ASCII data using General Mapping Tools (GMT)?

There are two "tricks" to using this dataset (g02156) with GMT: reading the data into a GMT grdfile, and knowing the projection to use when plotting data
over the gridded image.

First, the ASCII files have a 30-line header, then single-digit values all strung together. GMT's xyz2grd command needs the data reformatted a bit
to work. This shell script will do it:

#!/bin/sh
infile=$1
tail +31 $infile | perl -pe 'chomp; s/(.)/$1\n/g' |
xyz2grd new_data.txt -Gt.grd -ZBLa -V \
-R-12288000/12288000/-12288000/12288000 -I24000 -F
# end of script

The 'tail' command skips the header. The perl command puts each single digit on its own line.

Visualizing this data set with coastlines overplotted is tricky. First, the grid is defined beyond the horizon of the Earth's disk, and GMT doesn't handle that well. So, the following shell script makes a subset of the grid using the GMT program 'grdcut'. See the comments in the script for details on that issue. If you want
to plot the whole disk with coastlines, try trimming the grid so that the disk is tangent with the edge of the grid, then plot the coastlines with -R-180/180/$minlat/90, where $minlat is the latitude where you trimmed the disk to.

#!/bin/sh
if [ "$1" = "" ] ; then
echo "Usage: $0 input_ascii_file"
exit 1
fi
infile=$1
tail +31 $infile | perl -pe 'chomp; s/(.)/$1\n/g' |
xyz2grd -R-12288000/12288000/-12288000/12288000 \
-Gt.grd -ZBLa -V -I24000 -F
# Cut out middle of grid. Output grid is 272 x 272.

grdcut t.grd -R-3264000/3264000/-3264000/3264000 -Gt3.grd -V


# Plot the subsetted grid and overlay coastlines. Plot projection is linear
# for the grid, and polar stereographic for the coastlines. Plots are done at
# the same width to ensure the same scale.

out='ims2.ps'
width='16c'

# Plot grid
grd2cpt t3.grd -Cpolar -V > my_polar.cpt
grdimage t3.grd -Cmy_polar.cpt -B1000000WS -JX$width -V -K > $out

# In order to compute the lower-left and upper-right corner coordinates in
# lon/lat, the Proj.4 utility 'proj' can be used:

# $ proj +proj=stere +lon_0=-80 +lat_0=90 +lat_ts=60 +ellps=WGS84 -f %12.6f -I
# -3264000 -3264000
# -125.000000 47.667586
# 3264000 3264000
# 55.000000 47.667586

pscoast -R-125.000000/47.667586/55.000000/47.667586r \
-JS-80/90/$width -B30EN -W0.5p/green -Di -V -O >> $out

# End of script

In the script above, I showed how to use the Proj.4 utility "proj" to compute the lon/lat coordinates from the subsetted grid corners.
Here's how to do it with the GMT program "mapproject":

$ gmtset ELLIPSOID WGS-84 D_FORMAT %lf
$ mapproject -R-180/180/30/90 -Js-90/90/60/1:1 -C -F -I
-3264000 -3264000
225.000000 47.667585
3264000 3264000
45.000000 47.667585

If you have further questions that are general GMT questions, you might consider sending a message to the GMT-Help mailing list.

Provided by Bruce Raup

Last Updated May 2017