.. _gmt2: .. highlight:: tcsh ********************************* Generic Mapping Tools 2 ********************************* Now that we have a basic idea of how GMT works, we can start actually making maps. As with :math:`{x}`-:math:`{y}` plots, there are many options to specify, particularly the map projections. This lab gives several examples of how to plot different map projections, which we discuss before getting to the details of how to plot these projections in GMT. ================== Map Projections ================== Because the earth is spherical, it is not possible to draw a perfect map projection onto a flat surface. Different map projections are designed to preserve a particular quantity. Some of the most common quantities include: * Area * Shape * Distance * Bearing * Direction These things can all be measured independently, and we can choose a map projection that preserves at least one of these quantities. Because we may need different maps for different purposes, not all maps will want to preserve the same quantities. For instance, we may want to preserve areas in order to show some quantity per unit area. Equal area projections tend to skew shapes far from certain regions. If we want to maintain shapes (aka have a conformal map), then areas are not usually well-represented. Some maps are better for navigational purposes in that they show bearings correctly, or we may need to get accurate distance estimates. Thus, when considering which map projection to use, you need to consider the data that will be represented on the map. For instance, most of you have probably seen the Mercator projection. The Mercator projection projects the earth's surface onto a cylinder, and is conformal (it preserves small shapes and angles). It also has the property that straight lines are lines of constant course, useful for nautical navigation: .. image:: gmt2_0.* It is well known, partially because of its longtime use for navigation, partially because it is a nice rectangular projection (i.e. nice for hanging on walls or putting on rectangular book pages), and most recently, because it is the projection used in Google Maps and other online map applications. However, it greatly exaggerates areas near the poles, and makes Greenland appear to be the same size as Africa (when I was in Elementary School, I had a placemat showing a Mercator projection, and was amazed to learn that the area "rankings" I had deduced from the placemat were wildly inaccurate when I compared it to a globe). There are 5 types of projections available in GMT. We will consider an example of each of them (four in this lab, with the fifth type being the non-geographic :math:`{x}`-:math:`{y}` "projections" from the last lab). You can get more information on the various map projections on the web (the GMT manual has basic information on each projection, and Wikipedia also has good information on various projections). They break down into the following categories: * **Cylindrical Projections** (project the earth's surface onto a cylinder): Cassini, Cylindrical Sterographic, Miller Cylindrical Projectioin, Mercator, Oblique Mercator, Cylindrical Equidistant, Transverse Mercator, Universal Transverse Mercator, Cylindrical Equal-Area * **Conic Projections** (project the earth's surface onto a cone): Albers, Conic Equidistant, Lambert, (American) Polyconic * **Azimuthal Projections** (directions from central point are preserved): Lambert, Azimuthal Equidistant, Gnomic, Orthographic, General Perspective, General Stereographic * **Miscellaneous** (other qualities, or compromise projections that do not preserve any particular quantity): Hammer, Sinusoidal, Eckert IV, Eckert VI, Robinson, Winkel Tripel, Van der Grinten, Mollweide * **Non-Geographic Projections**: Polar Coordinates, Cartesian (:math:`{x}`-:math:`{y}`) coordinates ====================== A Simple Example ====================== In GMT, we use ``pscoast`` to draw landforms on maps. ``pscoast`` takes many of the same options as ``psbasemap``, and in many cases it is not necessary to call ``psbasemap`` because all of the necessary details are provided when calling ``pscoast``. However, I tend to call both on most of my maps, and specify ``-B`` when I call ``psbasemap`` rather than on ``pscoast`` (as I think of the grid lines belonging more to the axis frame than to the map within). Here I will give examples of both, and you should feel free to do whichever you prefer. The provided shell script ``gmt2.csh`` should produce four maps. The first one, ``gmt2_1.pdf`` is a simple example of a map of a Mercator projection of Oceania: .. image:: gmt2_1.* This map requires only one GMT call, to ``pscoast``. In GMT 4: :: pscoast -R90/180/-50/20 -JM5i -Ba10g10:."Mercator": -W -P > gmt2_1.ps and in GMT 5: :: gmt pscoast -R90/180/-50/20 -JM5i -Bxa10g10 -Bya10g10 -B+t"Mercator" -W -P > gmt2_1.ps As with the previous set of notes, I use shell variables extensively to make these commands more tractable. The actual command in the shell script is :: set code = "pscoast -R${xmin}/${xmax}/${ymin}/${ymax} -JM$size $bvals $continents $portrait > $outfile" eval "$code" Here I use the trick of setting the entire command to a variable in case we need to include a title that includes spaces. The command options should all be familiar to you from the first GMT lab, though their exact usage is a bit different here. Let's look at them one at a time: * ``-R90/180/-50/20`` are the ``-R`` options for this map: we give minimum and maximum longitude and latitude values to determine the map region. Here I am going from :math:`{90^\circ}` to :math:`{180^\circ}` in longitude, and :math:`{-50^\circ}` to :math:`{-20^\circ}` in latitude. Unlike when doing an :math:`{x}`-:math:`{y}` plot, the region can only take valid sets of geographical coordinates. Longitude must be in the range :math:`{(-180,180)}` or :math:`{(0,360)}`, and latitude must be in the range :math:`{(-90,90)}`. Certain map projections like the Mercator cannot include the full globe, as they are not defined at the poles. If you try to extend the region of such a projection to a pole, you will get an error message. * ``-JM5i`` are the ``-J`` options for this map, which sets the map projection. ``M`` stands for the Mercator projection, and the fact that it is a capital letter means that we set the map width of 5 inches (lower case means we would specify a length per degree longitude). Note that the map projection determines the map height, so we only need to specify the width. Some projections require additional information to define the projection that is independent of the region, though the Mercator projection is not among these (the min/max lat/lon from the ``-R`` option is all that is needed). * ``-Ba10g10:."Mercator":`` (GMT 4) or ``-Bxa10g10 -Bya10g10 -B+t"Mercator"`` (GMT 5) sets the boundary annotation to every 10 degrees, and the grid lines to every 10 degrees. The title of the plot is "Mercator" * ``-W`` tells GMT to draw coastlines using the specified line properties (none specified sets it to the defaults). You need to specify at least one additional option in ``pscoast`` regarding how to draw or fill the continents; here I set how to draw them (we will see later how to fill land and water areas). * ``-P`` sets output to portrait mode, and ``> gmt2_1.ps`` redirects output to file. I then convert to PDF using ``psconvert`` and delete the postscript file. This example is quite simple, as the Mercator projection is fairly easy to understand. ================================ A Non-Rectangular Projection ================================ The second example is of a Lambert Conic projection, which would produce a non-rectangular map if we specify the region using minimum and maximum longitudes and latitudes. However, what if you want your map to be rectangular? We can still make the final map output rectangular for this projection by changing the syntax for specifying ``-R``. This produces a map ``gmt2_2.pdf`` that looks like the following: .. image:: gmt2_2.* Note that while the map is rectangular, the boundaries are not lines of constant latitude and longitude. This effect can be produced in GMT as follows (GMT 4): :: pscoast -R-90/20/-60/35r -JL-75/27.5/20/35/5i \ -Gpalegreen -Sblue -Cskyblue -Na -Ir -W -Di \ -K -P > gmt2_2.ps psbasemap -R -J -Ba5g5:."Lambert Conic": -O -P >> gmt2_2.ps and for GMT 5: :: gmt pscoast -R-90/20/-60/35r -JL-75/27.5/20/35/5i \ -Gpalegreen -Sblue -Cskyblue -Na -Ir -W -Di \ -K -P > gmt2_2.ps gmt psbasemap -R -J -Bxa5g5 -Bya5g5 -B+t"Lambert Conic" -O -P >> gmt2_2.ps The commands using variables is the following: :: pscoast -R${xmin}/${ymin}/${xmax}/${ymax}r -JL${clon}/${clat}/${minlat}/${maxlat}/$size \ -G$land -S$ocean -C$lake $natbounds $rivers $continents -D$resolution \ $firstplot $portrait > $outfile set code = "psbasemap -R -J $bvals $finalplot $portrait >> $outfile" eval "$code" I have split this into two commands, one ``pscoast`` specifying how to draw the map, and ``psbasemap`` specifying how to draw the frame since the ``pscoast`` command is already fairly complicated. This example introduces several new syntax items: * ``-R-90/20/-60/35r`` In order to make a non-rectangular projection into a rectangular map, you need to specify upper left and upper right coordinates, rather than minimum and maximum values for each coordinate. The final lower case ``r`` tells GMT to interpret the ``-R`` command as lower left/upper right coordinates. Thus, I am telling GMT to plot the rectangular portion defined with the lower left coordinate at :math:`{(-90^\circ,20^\circ)}` and upper right coordinate at :math:`{(-60^\circ,35^\circ)}`. * ``-JL-75/27.5/20/35/5i`` tells GMT to make a Lambert Conic projection centered at :math:`{(-75^\circ,27.5^\circ)}`, with minimum latitude of :math:`{20^\circ}` and a maximum latitude of :math:`{35^\circ}`, and width of 5 inches. Note that this defines how the projection is constructed, but does not determine the final area that will be shown on the map (that is what ``-R`` specifies). Note that the two are independent of one another, though it is likely that the values will be related, as my central point is at the center of the map and the latitude range matches the range in the region that I have selected. * The next three options tell GMT how to color the map. ``-Gpalegreen`` says to color all land areas pale green, ``-Sblue`` says to color all water areas blue, and ``-Cskyblue`` says to color all lakes and rivers sky blue. If ``-C`` is not specified, lakes and rivers have the same color as water areas. You do not need to specify all three of these -- you can make a map that only colors land, or only colors water -- but you must specify one of ``-G``, ``-S``, or ``-W``. * The three options following the coloring options tell GMT how to draw map boundaries, including political borders (``-Na``), rivers (``-Ir``), and coastlines (``-W``). There are several options for political borders (``-N1`` draws international borders, ``-N2`` draws national borders within the Americas, ``-N3`` draws marine borders, and ``-Na`` draws all three). There are many options for drawing rivers. ``-Ir`` draws all permanent rivers, and for more options, see the manual. As we saw above, the ``-W`` option tells GMT how to draw coastlines, and I call ``-W`` with no options to use the default line style to draw the coastlines. * ``-Di`` specifies the map resolution. The following letter can be ``f`` (full), ``h`` (high), ``i`` (intermediate), ``l`` (low), or ``c`` (crude). The default is low, which is reasonable for maps at the scale of the first example. This example depicts a smaller region, so I have increased the resolution to intermediate. The other options should be familiar to you. The output is sent to the file ``gmt2_2.ps``. We then add the axis frame and grid lines by calling ``psbasemap``. These options should be familiar to you, here we call ``-Ba5g5:."Lambert Conic":`` or ``-Bxa5g5 -Bya5g5 -B+t"Lambert Conic"`` to annotate and draw grid lines every 5 degrees. We also give a title. The other options should be familiar, and note that we just give ``-R -J`` to use the same region and projection as the ``pscoast`` command. One important thing to recognize here was that we needed to put ``psbasemap`` after ``pscoast``. This is because successive GMT commands are plotted on top of any previous commands. Since we colored the water and continents, if we drew the grid lines before coloring the land and water, the grid lines would not show up. This is important to keep in mind when you construct complicated maps, as you need to build up the plot layer by layer with different commands. ========================== Plotting Data on a Map ========================== While making pretty maps is nice, in doing scientific research we would like to plot data on maps. Here is a basic example using ``psxy`` to draw plate boundaries (included in the file ``nuvel_plates.txt``) on a global map ``gmt2_3.pdf`` in the included shell script: .. image:: gmt2_3.* This plot uses the Robinson projection, a compromise projection designed to look "right" when plotting the entire globe and avoid the area distortion in the Mercator projection. The commands to make this map are as follows (GMT 4): :: pscoast -Rd -JN5i -N1 -A10000 -Dc -Gcoral -Sblue \ -X1.5i -Y1.5i -K -P > gmt2_3.ps psbasemap -R -J -Ba60g30/a30g15:."Robinson with Plate Boundaries": \ -K -O -P >> gmt2_3.ps psxy nuvel_plates.txt -R -J -W1p,red -m">" -O -P >> gmt2_3.ps For GMT 5: :: gmt pscoast -Rd -JN5i -N1 -A10000 -Dc -Gcoral -Sblue \ -X1.5i -Y1.5i -K -P > gmt2_3.ps gmt psbasemap -R -J -Bxa60g30 -Bya30g15 -B+t"Robinson with Plate Boundaries" \ -K -O -P >> gmt2_3.ps gmt psxy nuvel_plates.txt -R -J -W1p,red -O -P >> gmt2_3.ps where the actual commands in the shell script are :: pscoast -R$region -JN$size $natbounds $area -D$resolution -G$land -S$ocean \ -X$xshift -Y$yshift $firstplot $portrait > $outfile set code = "psbasemap -R -J $bvals $midplot $portrait >> $outfile" eval "$code" psxy $infile -R -J -W${pen},$color $marker $finalplot $portrait >> $outfile First, we draw the coastlines using ``pscoast``. Here we give a couple of options that you are not familiar with: * ``-Rd`` is shorthand for a global domain: ``-Rd`` is equivalent to ``-R-180/180/-90/90``. Also available is ``-Rg``, which is the same as ``-R0/360/-90/90``. * ``-JN5i`` means a 5 inch wide Robinson projection. * ``-A10000`` means to not draw any land areas with an area of less than 10000 km :math:`{^2}`. This option is useful on large global projections to avoid trying to draw tiny islands. * The other options should be familiar. We color the continents and water, draw international boundaries, translate the map (to prevent cutting off the title), choose portrait mode, and intend to add more output in later commands. Next, we draw the basemap frame to draw grids and label the axes. These commands should all be familiar to you. Remember that we do this second so that the grid lines are not obscured by the continents and oceans. The third command uses ``psxy`` to draw lines on the map. The command reads data from file (``nuvel_plates.txt``) and plots red lines of width 1 point on the map (``-W1p,red``). ``-R -J`` should be familiar to you, as well as ``-O -P``. The additional option for this command is ``-m">"``, which tells ``psxy`` how to separate different line segments in the input file. If you look at the input file, you will see it is a long list of longitude and latitude coordinates, with each segment separated by a line that begins with ``>``. The option ``-m">"`` tells GMT to use this to separate the different segments. The ``-m`` option has been deprecated in GMT 5, and the ``>`` character is the default for separating line segments in a data file, so that argument is not necessary in GMT 5. To handle this, I set the ``marker`` shell variable to be an empty string if GMT 5 is detected. This is a relatively simple example of how to plot data on a map. We will explore this topic in much greater detail in the next lab, showing how to plot additional types of data on map projections with much more complexity. ============================ A More Complicated Map ============================ The final example shows how to construct a map that is more difficult to make in GMT. This is mostly because the default version of the projection does not make the exact map that I desire, but there are sufficient tools that I can piece it together using several different commands (this is an example of how "Unix makes difficult tasks possible"). This is pretty tricky, so do not worry if you do not understand every detail. If you need to make this kind of map, you can use my script as a "cookbook" example and change the necessary pieces. One nice map projection for determining earthquake distances and back azimuths to a particular seismic station is the azimuthal equidistant projection. This projection is centered on a particular point, and all lines to that point are great circle paths (this is the "azimuthal" part) and all distances from the central point are accurately represented (this is the "equidistant" part). For instance, here is azimuthal equidistant projection centered on Memphis: .. image:: gmt2_4.* If we plot an earthquake epicenter on this projection, we can easily determine the distance and back azimuth. How do we make this map? :: psbasemap -R0/360/-90/0 -JS0/-90/90/5i \ -Ba15:."Azimuthal Equidistant": -K -P > gmt2_4.ps pscoast -Rg -JE-89.97/35.12/175/5i -Dc \ -Gdarkgray -K -O -P >> gmt2_4.ps psbasemap -Rg -JE0/-90/175/5i -Bg15/g30 -O -P >> gmt2_4.ps For GMT 5: :: gmt psbasemap -R0/360/-90/0 -JS0/-90/90/5i \ -Bxa15 -B+t"Azimuthal Equidistant" -K -P > gmt2_4.ps gmt pscoast -Rg -JE-89.97/35.12/175/5i -Dc \ -Gdarkgray -K -O -P >> gmt2_4.ps gmt psbasemap -Rg -JE0/-90/175/5i -Bxg15 -Byg30 -O -P >> gmt2_4.ps The code containing variable names is :: set code = "psbasemap -R0/360/-90/0 -JS0/-90/90/$size $bvals $firstplot $portrait > $outfile" eval "$code" pscoast -Rg -JE${clon}/${clat}/${dist}/$size -D$resolution \ -G$land $midplot $portrait >> $outfile psbasemap -Rg -JE0/-90/${dist}/$size $bvals2 $finalplot $portrait >> $outfile This is a bit tricky, because GMT's azimuthal equidistant projection does not draw the frame or grids as I want them to appear. However, we can piece together additional commands to draw the frame and grid lines. The first command draws the exterior basemap frame, and nothing else. Because this frame is not available in the azimuthal equidistant projection, I am instead using a polar stereographic projection centered at the South Pole (``-JS0/-90/90/5i``), but only plotting the frame. Here, the first two numbers are the longitude and latitude of the center point :math:`{(0,-90^\circ)}`, the third number is the horizon (chosen here to be :math:`{90^\circ}` so that the map is for the entire southern hemisphere). The final number is the map size (5 inches). ``-R0/360/-90/0`` says to plot the southern hemisphere centered at the south pole, which gives the frame with the correct azimuth annotations that I desire (every 15 degrees, specified by ``-Ba15:."Azimuthal Equidistant":`` or ``-Bxa15 -B+t"Azimuthal Equidistant"``), but without any grid lines (which I will add later). I choose portrait mode (``-P``) and will add more plotting commands later so I add ``-K``. Output is sent to ``gmt2_4.ps``. Once I have created the frame with ``psbasemap``, I now plot the actual map. I use ``pscoast`` to plot an azimuthal equidistant projection centered on Memphis. The option ``-JE-89.97/35.12/175/5i`` says that I want my projection centered at :math:`{(-89.97^\circ,35.12^\circ)}` (the geographic coordinates of Memphis), with a horizon of :math:`{175^\circ}` (meaning that I only plot points up to :math:`{175^\circ}` away, with a width of 5 inches. The depicted region is global (``-Rg``), and resolution is set to be coarse (``-Dc``). Continents are colored dark gray (``-Gdarkgray``), I use portrait mode (``-P``), and I append to an existing file, with the intention of continuing with further commands (``-K -O``). A few comments about this projection: in GMT 4 it is known to have some bugs that can lead to problems in depicting the continents. These have supposedly been solved in GMT 5. It also has trouble coloring the map if the antipode (the point :math:`{180^\circ}` away) is land. I had trouble getting this map to plot in GMT 4 when I colored both the land and the water, though some of the problems seem to be resolved in GMT 5, and I was able to successfully plot this map in both GMT versions. The final ``psbasemap`` command is a trick to draw the desired gridlines on the plot. GMT is only able to draw gridlines that are meridians or parallels on map projections. However, what I wanted for this plot was a set of great circles radiating from the center of the map, and lines of constant distance from the map center. To do this, I used a trick where I plot the grid for the azimuthal equidistant projection centered on the south pole. By choosing this particular map projection, the "meridians" become the great circles, and the "parallels" become the distance contours. (Bob Smalley deserves credit for this, I got this from his GMT notes.) To draw these gridlines, I call ``psbasemap`` for a global domain (``-Rg``), with the projection defined by ``-JE0/-90/175/5i``. This is the same size and horizon as the ``pscoast`` command, but because it is centered at the south pole :math:`{(0,-90^\circ)}` it will draw the grid lines in the places that I want. Note that ``-Bg15/g30`` or ``-Bxg15 -Byg30`` only tells it to draw the grid lines, and the :math:`{15^\circ}` increment in longitude matches up with the annotation interval specified in the first ``psbasemap`` command. The grid commands are appended onto the file. As I said earlier, I do not expect you to totally understand what I did here. This should illustrate how you can stick together different pieces of code to do something complex, and if you need to make this map projection in the future you can use this code as a cookbook example starting point. ============ Exercises ============ As with the last lab, you should try changing different options as you go through the scripts. Here are some ideas of things to try to change on some of the maps: * Plot the same projection for a different region. * Plot a different projection for the same region. Either use one of the other ones we covered today, or learn the syntax for a new projection (see the manual for details). * Change the grid and annotation spacing, or remove the grid entirely. * Add or take away political boundaries and rivers. * Change the fill colors or the coastline outline colors. * Change the resolution (it is best to do this on a small map, as global maps at high resolution create huge files). * Plot symbols on the map using ``psxy`` and standard input. * Add a text annotation to a map using ``pstext``.