Beyond the Mouse LAB 10: GMT 2
November 14, 16
Instructor: Jeff Freymueller
x7286 Elvey 413B jfreymueller@alaska.eduTA: Shanshan Li
Last Updated: November 7, 2017
Due: Tuesday Nov 22 before class
Reclaim your data! We follow up on GMT and introduce you to some of the tools that are beyond just mapping / plotting your data. GMT comes with a plethora of tools for data processing and analysis. This example focuses on one of those, for making a contour map.
Pre-lab reading
Read Smith and Wessel (1990) on the
procedure for fitting splines in tension. Original link is http://cosmos.hwr.arizona.edu/Technical/Smith-Wessel-1990.pdf.
Exercise 1
Here are script and data from the lecture: 10_gmt2_example.tar.gz, 10_gmt2_example.zip. Play with the example. Try some different filtering options, and other plotting options. Turn in one script and plot and explain what you changed relative to the default script.
Exercise 2
Now we will make a contour map. The map will contour uplift rates from a model of Glacial Isostatic Adjustment (GIA). Making a contour map is a fairly common thing, and this data set serves as a good example. There are three main steps to making a contour map, and this exercise will walk you through the steps:
- Organize the data into regular blocks (
blockmeanorblockmedian). - Fit a gridded surface to the data (
surface). - Plot a basemap and draw contours (
pscontour).
First, download the file Hu2013.uplift.txt,. This is your date file. It has four columns: longitude, latitude, uplift rate (mm per year), and a weight value. The 4th value is used as a weight for the data in the averaging process. Because these are all model computations, the weights are the same for all data.
The first step is to run blockmean or blockmedian
(you can choose, and maybe try both) to organize
the data into a more uniform grid. The data file contains the model values at
a regular grid of points plus the model values at irregularly sampled GPS site
locations. This step makes the data more uniform, which helps in the surface fitting.
In addition to the filename, you will need to specify some switches on the command:
- -R-180/-120/53/71 (this is the range of long-lat values we will consider)
- -Ix/y (the block spacing in long/lat. Try values in the range 0f 0.25-1 degree)
- -WI -V (this says to use the weights and be verbose)
Redirect the output to a file. You will use this as the input for the next step. You can also use a pipe.
Now run surface to make a grid file with the gridded fit to the data.
You will need to give it the file name and several switches on the command line.
- -R-180/-120/53/71 (this is the range of long-lat values we will consider)
- -Gfilename.grd (output file name, replace "filename" with something useful)
- -I5m/2m (the grid spacing in long/lat. Should be smaller than that used in blockmean)
- -T (tension value, between 0 and 1, try 0.25-0.75)
Finally, make your map, You will call pscoast to make a basemap, and then use
grdcontour to draw the contours. Here you have a lot of freedom. The lines of
code below are one way to draw them (you will need to change the input .grd file name
and the output file name to match your script). The first line plots all contours
for zero and positive values using red lines (-L selects which range of contours to
draw). It plots a contour line every 1 mm/yr (-C1), and every other contour line
is a major contour that is thicker and labeled (-Gd2i). The -W switches define the
line thickness and color for the minor (-Wc) and major (-Wa) contours. The second
and third lines draw contours from -1 to -2 and less than -2, respectively. I used three
calls instead of two so that the zero contour would end up being drawn only once.
gmt grdcontour ${base}.grd -S -R -J -L0/40 -Gd2i -C1 -A2s12 -Wc3,255/0/0 -Wa6,255/0/0 -V -O -K -P >> $PS_OUT
gmt grdcontour ${base}.grd -S -R -J -L-2/-1 -Gd2i -C1 -A2s12 -Wc3,0/0/255 -Wa3,0/0/255 -V -K -O -P >> $PS_OUT
gmt grdcontour ${base}.grd -S -R -J -L-40/-2 -Gd2i -C1 -A2s12 -Wc3,0/0/255 -Wa6,0/0/255 -V -O -P >> $PS_OUT
Once you have a working script, play around with some of the values and see how the contour map changes. Make a second version of the contour map in which you add symbols to show where the data values come from (use psxy to draw those symbols). Now assess your map. Are there areas where the contours seem unreasonable, or where countours are drawn that are not really constrained by the data? Turn in versions of the countour map with and without the data loctions plotted.
For extra credit, you can make another map that compares the input data and the gridded version of the data set. You can use the program grdtrack to compute values from the grid file that you have already created. After you have those values, you could either make a figure that compares the input and gridded data values plotted as vectors (using psvelo), or take the difference between the input and gridded values and plot a histogram of these with pshistogram. Doing either of these will require using awk to do some reformatting of files, or for the mathematical manipulations.
Dr. Jeffrey T. Freymueller
Professor of Geophysics
Geophysical Institute
University of Alaska, Fairbanks
Fairbanks, AK 99775-7320