New Mexico Tech
Earth and Environmental Science

ERTH 455 / GEOP 555 - Geodetic Methods

Lab 12: GMT - Generic Mapping Tools

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on redoubt today!

Introduction

We've been dealing a lot with data that, given the nature of this course, has some geospatial meaning. Therefore, it is important to be able to create maps that present the data in space; maybe spiced up with some auxiliary information such as topography, coastlines, national and international boundaries, rivers and so forth. Ideally, these maps would also be aesthetically pleasing and of high enough quality to be used in publications, on posters, slides, etc.

GMT or the Generic Mapping Tools fit the bill. In fact, it does all of the above and provides a range of tools for data processing (you've already used grdmath!). From their website:

"GMT is an open source collection of about 80 command-line tools for manipulating geographic and Cartesian data sets (including filtering, trend fitting, gridding, projecting, etc.) and producing PostScript illustrations ranging from simple x–y plots via contour maps to artificially illuminated surfaces and 3D perspective views; the GMT supplements add another 40 more specialized and discipline-specific tools. GMT supports over 30 map projections and transformations and requires support data such as GSHHG coastlines, rivers, and political boundaries and optionally DCW country polygons."

Now here's the deal: As it says above, GMT is a collection of about 70-80 programs / modules that all certain jobs. You won't get around spending a lot of time working with the man-pages (e.g, $> man gmt) and/or the online documentation. This lab is an attempt to introduce you to the concept of working with GMT; we won't even get past scratching the surface of GMT's abilities. Consider achieving mastery of GMT a lifelong journey - there's always some new stuff to learn! To convince you, here's a list of the current GMT modules (type $> man gmt to get this list); each of these has their own man page!

       blockmean      L2 (x,y,z) data filter/decimator
       blockmedian    L1 (x,y,z) data filter/decimator
       blockmode      Mode (x,y,z) data filter/decimator
       filter1d       Filter 1-D data sets (time series)
       fitcircle      Finds the best-fitting great circle to a set of points
       gmt2rgb        Convert Sun rasterfile or grid to r, g, b grids
       gmtconvert     Convert between ASCII and binary 1-D tables
       gmtdefaults    List the current default settings
       gmtmath        Mathematical operations on data tables
       gmtset         Set individual default parameters
       gmtselect      Extract data subsets based on spatial criteria
       grdfilter      Filter 2-D data sets in the space domain
       grd2cpt        Make a color palette table from grid files
       grd2xyz        Conversion from 2-D grid file to table data
       grdblend       Blend several partially over-lapping grid files onto one grid
       grdclip        Limit the z-range in gridded data
       grdcontour     Contouring of 2-D gridded data
       grdcut         Cut a sub-region from a grid file
       grdedit        Modify header information in a 2-D grid file
       grdfft         Operate on grid files in the wavenumber (or frequency) domain
       grdgradient    Compute directional gradient from grid files
       grdhisteq      Histogram equalization for grid files
       grdimage       Produce images from 2-D gridded data
       grdinfo        Get information about grid files
       grdlandmask    Create mask grid file from shoreline data base
       grdmask        Reset nodes outside a clip path to a constant
       grdmath        Mathematical operations on grid files
       grdpaste       Paste together grid files along a common edge
       grdproject     Project gridded data onto a new coordinate system
       grdreformat    Converting between different grid file formats
       grdsample      Resample a 2-D gridded data set onto a new grid
       grdtrend       Fits polynomial trends to grid files
       grdtrack       Sampling of 2-D data set along 1-D track
       grdvector      Plot vector fields from grid files
       grdview        3-D perspective imaging of 2-D gridded data
       grdvolume      Volume calculations from 2-D gridded data
       greenspline    Interpolation using Green's functions for splines in 1-3 dimensions
       makecpt        Make color palette tables
       mapproject     Forward or inverse map projections of table data
       minmax         Find extreme values in data tables
       nearneighbor   Nearest-neighbor gridding scheme
       project        Project data onto lines/great circles
       ps2raster      Crop and convert PostScript files to raster images, EPS, and PDF
       psbasemap      Create a basemap plot
       psclip         Use polygon files to define clipping paths
       pscoast        Plot coastlines and filled continents on maps
       pscontour      Contour xyz-data by triangulation
       pshistogram    Plot a histogram
       psimage        Plot images (EPS or Sun raster files) on maps
       pslegend       Plot legend on maps
       psmask         Create overlay to mask out regions on maps
       psrose         Plot sector or rose diagrams
       psscale        Plot gray scale or color scale on maps
       pstext         Plot text strings on maps
       pswiggle       Draw time-series along track on maps
       psxy           Plot symbols, polygons, and lines on maps
       psxyz          Plot symbols, polygons, and lines in 3-D
       sample1d       Resampling of 1-D table data sets
       spectrum1d     Compute various spectral estimates from time-series
       splitxyz       Split xyz-files into several segments
       surface        A continuous curvature gridding algorithm
       trend1d        Fits polynomial or Fourier trends to y = f(x) data
       trend2d        Fits polynomial trends to z = f(x,y) data
       triangulate    Perform optimal Delaunay triangulation and gridding
       xyz2grd        Convert equidistant xyz data to a 2-D grid file

The way this will work is that I'll build an example that you will follow along. Then you'll take some of your own products (description below) and make your own map!

1. A simple map of the US

Let's get our feet wet! We'll make a map of the US in Cassini projection.

    $> gmt pscoast -R-125/-70/24/50 -Jc-97.5/37/0.3 -Ggrey -Sblue -P -Dc > us.ps
    $> evince us.ps &
	

You will see a simple map with a crude coastline of the US. Here is what the command line parameters mean:

    -R-125/-70/24/50   The map region be between -125 West, -70 deg East, 24 deg south and 50 degrees North. 
                       Order is obviously important!
    -Jc-97.5/37/0.3    Select a Cassini projection, centered on longitude -97.5 and latitude 37 with 
                       a scaling factor of 0.3 (not the best projection, we'll change that later)
    -Ggrey             Make the land area grey
    -Sblue             Make the wet areas blue
    -P                 Plot in Portrait mode
    -Dc                Make a crude coastline
 	

A crusde coastline is not that great! -D comes with a few parameters (check man page!) let's try full resolution - tip: keep evince open, it should reload your new map automatically!

    $> gmt pscoast -R-125/-70/24/50 -Jc-97.5/37/0.3 -Ggrey -Sblue -P -Df > us.ps
	

There's a lot more water! And the coastlines look better now! Let's add International and State boundaries

    $> gmt pscoast -R-125/-70/24/50 -Jc-97.5/37/0.3 -Ggrey -Sblue -P -Df -N1/1,black -N2/0.5,100/100/100 > us.ps
	

Here -N has the first parameter 1,2,3,a for International, State, Marine, or all boundaries. What follows is a separator / and a pen definition which has the comma separated arguments size,color,texture. Here we are only changing the size and the color. GMT has a list of predefined color names, or you can fall back on RED/GREEN/BLUE triplets; where RED,GREEN,BLUE are values between 0-255.

Let's finish this up with the -B argument to add some labels and the map axes. This command line argument has been changed for the new version 5 of GMT, but it remains probably the most complex argument. It is fully described in the psbasemap man pages $> man psbasemap. Please read this for the full info.

    $> gmt pscoast -R-125/-70/24/50 -Jc-97.5/37/0.3 -Ggrey -Sblue -P -Df -N1/1,black -N2/0.5,100/100/100 -Ba10f5:."Map of the contiguous US":WenS > us.ps
	

Here we set the axes labels for every 10 degrees and the fine tickmarks at 5 degrees (a10f5) and we add a Map title (:."Map of the contiguous US":); not that the colons are important!! And lastly we turn axes labels on for the West and South boundaries, and off for east and north (WenS). Let's change the spacing of the labels for the longitude and make axis labels appear at the East border too.

    $> gmt pscoast -R-125/-70/24/50 -Jc-97.5/37/0.3 -Ggrey -Sblue -P -Df -N1/1,black -N2/0.5,100/100/100 -Ba15f1/a10f5:."Map of the contiguous US":WEnS > us.ps
	

It is important to note that, by default, GMT always expects longitude arguments first and then latitude arguments, which conforms with the math X,Y order of coordinates. So, we added a15f1/ to the -B string and changed the small "e" to "E" to turn axes labels on in the east. As one last step, change the projection to mercator:

    $> gmt pscoast -R-125/-70/24/50 -Jm-97.5/37/0.3 -Ggrey -Sblue -P -Df -N1/1,black -N2/0.5,100/100/100 -Ba15f1/a10f5:."Map of the contiguous US":WEnS > us.ps
	

What changed you may ask? The -Jc became -Jm. The arguments behind this depend on the projection and I refer you, again, to the psbasemap documentation to learn about offered projections.

2. Add Stuff to the map

So far, we've made use of GMT's internal databases and not really added much of our own information to the map. As soon as we start doing this, we should move to a shell script. Here's a basic one that's based on the above command, but changes a few things:

#!/bin/tcsh

#some setup - check gmtdefaults for meaning.
gmt gmtset PS_MEDIA             Custom_500x600
gmt gmtset FONT_ANNOT_PRIMARY   14
gmt gmtset MAP_FRAME_TYPE       PLAIN

#that's the filename where the results go.
set OUTPUT = us.ps

# Make a basemap to define map region (-R) and projection (-J)
# in following commands you can use only -R -J without the arguments
# of this first call. GMT will then use the previous values.
# this step is not always necessary
gmt psbasemap -Y2i -R-125/-90/24/50 -Jm-107.5/37/0.5 -P \
            -Ba15f5/a10f5:."Map of the western US":WenS -K               >  $OUTPUT

# add a topography. I downloaded the e.grd at some point
# from the web. I created e_i5.grd with some aximuth to 
# get illumination
gmt grdimage /usr/local/GMT5/share/dbase/e.grd      \
            -I/usr/local/GMT5/share/dbase/e_i5.grd  \
            -C/usr/local/GMT5/share/cpt/usa.cpt \
            -R -J -O -K                                                 >> $OUTPUT

#add state borders (-N flag)
gmt pscoast -R -J -B -Slightblue -Df -N1/1,black -N2/0.25,black -O -K   >> $OUTPUT

#add white circles (-Sc) with a certain radius at the 
#locations given between the <<END and END labels 
gmt psxy -R -J -B -Sc0.15 -Gwhite -O -K <<END                           >> $OUTPUT 
-106.609991 35.110703
-106.899424 34.061759 
END

# add some text labels (where are the cities?) 
gmt pstext -R -J -O -F+f10p,Helvetica,white+jLM+a0 -N -K <<END         >> $OUTPUT
-106.3 35.110703 Albuquerque
-106.5 34.061759 Socorro
END

# add scale stick - how long etc - NOTE THAT THERE IS NO -K here!!!
gmt psbasemap -R -J -B -L-120/26/37/500 -O -K                           >> $OUTPUT

#OLD -D convention!!
gmt psscale -C/usr/local/GMT5/share/cpt/usa.cpt -D2.75i/-0.5i/5.5i/0.25ih -O -B1000:topography:/:m:  >> $OUTPUT

# convert to PDF 
ps2pdf $OUTPUT
	

You can download the script above from here. Save it to a convenient spot and execute it:

    $> chmod +x us_map.gmt
    $> ./us_map.gmt
	

The chmod command makes sure the script is executable for you, the second line executes it. If your evince window with us.ps is still open, the new file should be displayed automatically if you are working in the same directory, else you need to open us.ps to see the result.

NOTE: For some reason psscale is still executed with the old -D syntax. I don't know exactly when that changed (check GMT 4.5 documentation for that!)

3. Plotting Grid files

The topography we're plotting in the US map is a gridded file. In principle, there are 2 steps to plotting such 2d data:

Let's do this for the DEM you used in homework 2 - start out by copying this into your current directory and then make a color palette:

    $> gmt grd2cpt dem.grd -Crelief -Z > t.cpt
	

You now have your own color palette t.cpt that is based on the color scheme of the relief standard palette. An overview of standard palettes is given, for example, here. Note that on redoubt you can find the files that define the standard palettes in /usr/local/GMT5/share/cpt. Rainbow is probably best for interferograms, but feel free to experiment!

Next, let's plot this!

    $> gmt grdimage dem.grd -Ct -JM10 `gmt grdinfo -I.0001 dem.grd` -Ba0.5f0.01:."DEM":WSen -P > dem.ps
	

Note that the call in the backticks extracts the map region from the DEM and returns it in the "-RW/E/S/N" format, here to the nearest 0.0001 degree. (execute this on the command line by itself to see the result)

This looks alright. Let's try adding some illumination!

    $> gmt grdgradient dem.grd -A0 -Nt -Gdem_int.grd
	

grdgradient computes the directional derivative for a 2D grid file; the result can be used to illuminate the DEM. The -A flage here is the aximuth (in degrees clockwise from north). -N is the normalization, here we use the topography values, i.e. we normalize by the elevation of the terrain. The input file is dem.grd, the result is saved to dem_int.grd. Let's use that for our DEM:

    $> gmt grdimage dem.grd -Ct -JM10 `gmt grdinfo -I.0001 dem.grd` -Ba0.5f0.01:."DEM":WSen -P -Idem_int.grd > dem.ps
	

Looks pretty fancy, right?

Your task

Use the last command above to create a DEM illuminated map of your phasefilt_ll.grd. Here are the steps:

Beyond that, get as fancy as you want, for example: use the us_map.gmt script as a blue print for your own map in which you will plot the DEM from homework 2 and then add the phasefilt_ll.grd with a rainbow color palette colorbar, label things, etc. Realize that this will require some experimentation on your end, and man-page reading.

Deliverables: (submit via canvas!)

rg <at> nmt <dot> edu | Last modified: November 21 2017 21:17.