Generic Mapping Tools Hints and Tips

If you are unaware of what GMT is then stop reading and visit the GMT homepage. In short GMT is what I always wanted when I was making maps as a graduate student. Although slightly esoteric (it is a map making tool after all, and cartography is a complex science, ask a cartographer or school teacher about Mercator projections!) it is probably the most useful coding I have ever learnt (excluding Vi and Vim of course!).

Bascially the most efficient way of coding in GMT is to write a series of GMT commands in a shell script. I use the tcsh shell, but bash and sh are just fine. Shell scripting is waaaaaayyyy beyond the scope of these tips, so please don't ask for help in writing shell scripts.

What follows are just some general tips (primarily to remind me!) in how to write efficient GMT scripts. I owe much of my initial GMT learning to my co-workers at UCSD, especially Jen Bowen and Debi Kilb. My acolyte period has been helped by Kris Walker and Kent Lindquist.

Define as many variables as you can

  1. set TEMPDIR = .
  2. set PS = $TEMPDIR/anza_jacinto.ps
  3. set BASE = /opt/gmt/bin
  4. set GRIDS = /path/to/source/grdfiles/
  5. set PS_GLOBE = $TEMPDIR/global.ps
  6. set REGION = -117/33.2/-116.2/33.8r
  7. set CENTER = -116.6/33.5/28

Psbasemap scale

These have been notoriously tricky for me to master. The manual page says:

-L[f][x]lon0/lat0[/slon]/slat/length[m|n|k][:label:just][+ppen][+ffill]
-L[f][x]///[m|n|k], append m, n, or k for miles, nautical miles, or km [Default]

So I typically code for a fancy scale...

  1. $BASE/psbasemap -X0 -Y0 -P -R$REGION -JE$CENTER -V -Bg0.5a0.5:"":/g0.5a0.5:""::".":WeSN -K -Lf-116.8/33.5/34/20k >> $PS

Here f means a 'fancy' scale, -116.8 is the longitude of the center of the scale, 33.5 is the latitude of the center of the scale, 34 is the latitude at which the scale is calculated. This will differ if you locate it in different parts of the map, depending on the type of projection you use. I typically use an equidistant (-JE) projection, therefore this latitude is extremely important. The final integer is the length of the scale and in this example I defined 20 kilometers.

Clip regions

Clip regions are for places on the planet where you need to overlay, or mask, what is plotted underneath. I most commonly use it for recreating 'land-only' topography over an area that is below sea level (and therefore plots as a 'wet' area). These locations, although uncommon, are found in several areas around the world, such as the Salton Sea and the Caspian Sea. You need an xy file that defines the polygon that will be clipped. This can as simple as 4 points (a rectangle) or as complex as a multi-sided polygon.

  1. # -- START: PLOT TOPO,BATHYMETRY CORRECTLY (including Salton Sea) --
  2. # Salton Sea co-ords -R-116.8/-115/32/34
  3. # Define a clip region
  4. $BASE/psclip saltonsea.xy -P -R$REGION -JE$CENTER -V -K >> $PS
  5. # Make the Salton Sea be "land-only" and put into the clipping region
  6. $BASE/grdimage saltonsea.grd -V -P -R$REGION -JE$CENTER -Clandonly.cpt -Isaltonsea.grad -O -K >> $PS
  7. # Color the Salton Sea Blue
  8. $BASE/pscoast -V -P -R$REGION -JE$CENTER -C201/244/255 -Df -O -K >> $PS
  9. # Close psclip
  10. $BASE/psclip -C -K -O >> $PS
  11. # -- END: PLOT TOPO,BATHYMETRY CORRECTLY (including Salton Sea) --

GRD problems

Check a grd file is valid with grdinfo. You might need to add a '=' to correctly parse the grd file:

  1. $BASE/grdimage anza_3.grd=2 -P -R$REGION -JE$CENTER -Clandonly.cpt -Ianza_3.grad -V -K -E100 -X2 -Y2 >> $PS

Using psmeca

I plot earthquake focal mechanisms on maps occasionally. These are derived and defined in a variety of formats: the Harvard CMT; the SCEC; and Aki & Richard (1980). To plot the 'beach balls' (as they are represented on a map) I use GMT's built in command psmeca (manual page here). I always opt for the easiest method, which turned out to be using psmeca's partial data focal mechanism option (-Sp). Using the SCEC catalog is better for me as well because the Harvard catalog only provides focal mechanisms for events with magnitude > 5.5. This is not so useful for local events in our Anza network area, where events < 4 are most common.

Using an example from SCEC, you might have:

** SCSN Moment Tensor Solution Message **

REAL-TIME SOLUTION: OPERATOR REVIEWED
Reviewed On: 05/10/2008 3:9:57

Inversion Method: Complete Waveform
Number of Stations used: 4
Stations: CI.GLA CI.GMR CI.DVT CI.SWS 
                   
Real-Time Solution:
-------------------
Event ID    : 14367532
Magnitude   : 4.07
Depth (km)  : 2.5
Origin Time :  05/09/2008 22:38:07:880
Latitude    : 33.45
Longitude   : -116.46
Further Information at: http://pasadena.wr.usgs.gov/recenteqs/Quakes/ci14367532.htm
                            
SCSN Moment Tensor Solution:
----------------------------
Moment Magnitude     : 4.39
Depth (km)           : 14
Variance Reduction(%): 17.81
Quality Factor       : C
   (A : Mw, MT good enough for distribution)
   (B : Mw only good enough for distribution)
   (C : Solution needs review before distribution)
                                                  
                                             
Best Fitting Double Couple and CLVD Solution:
---------------------------------------------------
                                                                
Moment Tensor: Scale = 10**21 Dyne-cm
    Component   Value
       Mxx      -36.3
       Mxy      -11.8
       Mxz      -18.1
       Myy      29.5
       Myz      -27.9
       Mzz      6.82

Best Fitting Double Couple Solution:              
--------------------------------------------------
 Moment Tensor: Scale = 10**22 Dyne-cm
    Component   Value
       Mxx     -3.626
       Mxy     -1.177
       Mxz     -1.739
       Myy      2.927
       Myz     -2.817
       Mzz      0.699
 
 Principle Axes:
    Axis    Value   Plunge   Azimuth
      T     4.842      34      270
      N     0.000      46      136
      P    -4.842      24       18
 
 Best Fitting Double-Couple:
    Mo = 4.84E+22 Dyne-cm
    Plane   Strike   Rake   Dip
     NP1      322     136    84
     NP2       58       8    46
 
 Moment Magnitude = 4.39
                                               
                                               
                    -------                    
              -------------------              
           ----------------   ------           
         ###--------------- P --------         
       #######-------------   ----------       
      ##########-------------------------      
     #############-----------------------#     
    ###############----------------------##    
    #################-------------------###    
   ###################-----------------#####   
   ######   ############---------------#####   
   ###### T #############------------#######   
   ######   ###############---------########   
   #########################-------#########   
    #########################----##########    
    ##########################-############    
     #######################----##########     
      ##################---------########      
       ---------------------------######       
         --------------------------###         
           -------------------------           
              -------------------              
                    -------                    
                                               
     Lower Hemisphere Equiangle Projection
 
=============  Station Information ==============
                                                 
Name       Distance     Azimuth     VR      ZCore
-------------------------------------------------
CI.GLA      158.446    105.747   12.014      6.00
CI.GMR      165.396     26.282   23.070     62.00
CI.DVT       93.894    158.958   19.487     41.00
CI.SWS       83.806    132.175   18.571     13.00

The part I use for psmeca is:

 Best Fitting Double-Couple:
    Mo = 4.84E+22 Dyne-cm
    Plane   Strike   Rake   Dip
     NP1      322     136    84
     NP2       58       8    46

Onto the example. You need to put your focal mechanisms in an external data file:

  1. # So use last resort, Sp partial data
  2. # lon lat depth strike1 dip1 strike2 -1/+1 magnitude bblon bblat text
  3. -116.59 33.63 14 36 77 298 1 5.2 -116.59 33.63 20050612Mw5.2
  4. # -116.51 33.51 15.2 48 90 318 1 5.1 -116.51 33.51 20011031Mw5.1
  5. # Guess you need to have a dip for one of the planes
  6. -116.51 33.51 15.2 318 33 48 1 5.1 -116.51 33.51 20011031Mw5.1
  7. -116.46 33.44 5.8 49 88 318 1 3.78 -116.46 33.44 20080501Mw4.2
  8. -116.46 33.45 2.5 322 84 58 1 4.39 -116.46 33.45 20080509Mw4.1

Then in your GMT script, you call psmeca with the correct option to the 'S' argument:

  1. $BASE/psmeca moment_tenors.txt -R$REGION -JE$CENTER -V -Sp0.5 -P -O >> $PS

Et voilà, you get focal mechanisms added to your map at the location specified in the bblon and bblat fields. The 'beach balls' are scaled by magnitude. If you want them all the same size, use the same magnitude in the data file.

Further reading