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.