Note that you DO have to work on redoubt today!
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!
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.
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:
psbasemap to set up ... basemap
#!/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_map.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_map.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!)
In principle, there are 2 steps to plotting 2d data:
$> 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 normaization, 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?
Use the last command above to create a DEM illuminated map of your phase_mask_ll.grd.
Then 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 phase_mask_ll.grd
with a rainbow color palette colorbar, label things, etc.
Realize that this will require some experimentation and man-page reading.
us_map.gmt scriptrg <at> nmt <dot> edu | Last modified: November 11 2015 21:48.