The Generic Mapping Tools
Version 4.5.11
Technical Reference and Cookbook
by
Pål (Paul) Wessel
School of Ocean and Earth Science and Technology
University of Hawai’i at Mānoa
and
Walter H. F. Smith
Laboratory for Satellite Altimetry
NOAA/NESDIS
November 2013
The Generic Mapping Tools (GMT) could not have been designed without the generous support of several people. We gratefully acknowledge A. B. Watts and the late W. F. Haxby for supporting our efforts on the original version 1.0 while we were their graduate students at Lamont-Doherty Earth Observatory. Doug Shearer and Roger Davis patiently answered many questions over e-mail. The subroutine gauss was written and supplied by Bill Menke. Further development of versions 2.0–2.1 at SOEST would not have been possible without the support from the HIGP/SOEST Post-Doctoral Fellowship program to Paul Wessel. Walter H. F. Smith gratefully acknowledges the generous support of the C. H. and I. M. Green Foundation for Earth Sciences at the Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California at San Diego. GMT series 3.x, 4.x, and 5.x owe their existence to grants EAR-93-02272, OCE-95-29431, OCE-00-82552, OCE-04-52126, and OCE-1029874 from the National Science Foundation, which we gratefully acknowledge.
We would also like to acknowledge the feedback we have received from many of the users of earlier versions. Many of these suggestions have been implemented, and the bug reports have been useful in providing more robust programs. Specifically, we would like to thank Michael Barck, Manfred Brands, Stephan Eickschen, Ben Horner-Johnson, John Kuhn, Angel Li, John Lillibridge, Andrew Macrae, Alex Madon, Greg Neumann, Lloyd Parkes, Ameet Raval, Georg Schwarz, Richard Signell, Peter Schimidt, Dirk Stoecker, Eduardo Suárez, Mikhail Tchernychev, Malte Thoma, David Townsend, Garry Vaughan, William Weibel, Florian Wobbe, and many others, including their advice on how to make GMT portable to a wide range of platforms. John Lillibridge provided the original example 11; Hanno von Lom helped resolve early problems with DLL libraries for Win32; Lloyd Parkes enabled indexed color images in PostScript; Kurt Schwehr maintains the Fink packages; Wayne Wilson implemented the full general perspective projection; and William Yip helped translate GMT to POSIX ANSI C and incorporate netCDF 3. The SOEST RCF staff (Ross Ishida, Pat Townsend, and Sharon Stahl) provided valuable help on Linux, web, and CGI script issues.
Honolulu, HI, Silver Spring, MD, Cornish, NH, and Faro, Portugal, November 2013
Starting with GMT version 3.2, all GMT documentation was converted from Microsoft Word to LATEX files. This step was taken for a number of reasons:
Please send email to the GMT help list if you find errors or inconsistencies in the documentation.
If you feel it is appropriate, you may consider paying us back by citing our EOS articles on GMT (and perhaps
also our Geophysics article on the GMT program surface) when you publish papers containing results or
illustrations obtained using GMT. The EOS articles on GMT are
The article in Geophysics on surface is
GMT includes some code supplied by others, in particular the Triangle code used for Delaunay triangulation. Its author, Jonathan Shewchuk, says
“If you use Triangle, and especially if you use it to accomplish real work, I would like very much to hear from you. A short letter or email (to jrs@cs.cmu.edu) describing how you use Triangle will mean a lot to me. The more people I know are using this program, the more easily I can justify spending time on improvements and on the three-dimensional successor to Triangle, which in turn will benefit you.”
A few GMT users take the time to write us letters, telling us of the difference GMT is making in their work. We appreciate receiving these letters. On days when we wonder why we ever released GMT we pull these letters out and read them. Seriously, as financial support for GMT depends on how well we can “sell” the idea to funding agencies and our superiors, letter-writing is one area where GMT users can affect such decisions by supporting the GMT project.
Copyright ©1991 – 2013 by Paul Wessel and Walter H. F. Smith
The Generic Mapping Tools (GMT) is free software; you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software Foundation.
The GMT package is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See
the file LICENSE.TXT in the GMT directory or the GNU General Public License for more details.
Permission is granted to make and distribute verbatim copies of this manual provided that the copyright notice
and these paragraphs are preserved on all copies. The GMT package may be included in a bundled distribution of
software for which a reasonable fee may be charged.
The Generic Mapping Tools (GMT) does not come with any warranties, nor is it guaranteed to work on your computer. The user assumes full responsibility for the use of this system. In particular, the School of Ocean and Earth Science and Technology, the National Oceanic and Atmospheric Administration, the National Science Foundation, Paul Wessel, Walter H. F. Smith, or any other individuals involved in the design and maintenance of GMT are NOT responsible for any damage that may follow from correct or incorrect use of these programs.
In reading this documentation, the following provides a summary of the typographic conventions used in this document.
While GMT has served the map-making and data processing needs of scientists since 19881, the current global use was heralded by the first official release in EOS Trans. AGU in the fall of 1991. Since then, GMT has grown to become a standard tool for many users, particularly in the Earth and Ocean Sciences. Development has at times been rapid, and numerous releases have seen the light of day since the early versions. For a detailed history of the changes from release to release, see file ChangeLog in the main GMT directory. For a nightly snapshot of ongoing activity, see the online ChangeLog page.
The success of GMT is to a large degree due to the input of the user community. In fact, most of the capabilities and options in GMT programs originated as user requests. We would like to hear from you should you have any suggestions for future enhancements and modification. Please send your comments to the GMT help list.
GMT 4.x will continue to see corrections of legacy bugs and problems, including Windows DLL silliness. The GMT 4.x releases are mostly bug-fixes as all development is now focussed on GMT 5; this new series, released Nov-5, 2013, is distinguished by being completely restructured so as to allow developers to call high-level GMT processes from a variety of programming environments. Below is a brief history of the development milestones in the 4.x series. We expect to maintain the GMT 4 series for some time until GMT 5 has reached a stable state. As of July 27, 2011, GMT 4 is under subversion control.
]
Note: Due to a few issues we had an aborted update to 4.5.10 that was briefly released, then retracted; hence the 4.5.11-numbered version. Another service release with bug-fixes only. The only non-bug change was adding the latest dimensions for recent Sandwell/Smith img files that go up to 85°, and adding definition file dat.def for mgd77 ASCII DAT format to the x2sys supplement. We also had to modify the –S option in pscontour.c to address a bug. This GMT release also coincides with the latest GSHHG release version 2.2.4 which adds a few missing lakes to California and fixes an error in the Baffin Island coastline and removes skinny spikes from numerous features. Below is the list of bug corrections for individual library files or programs:
]
Predominantly a bug-fix release, we also have made some changes to GSHHS. First, GSHHS is now called GSHHG, the “Global Self-consistent Hierarchical High-resolution Geography”, since GSHHG contains both shorelines as well as political boundaries and rivers. GSHHG, being required by both GMT 4 and 5, is now released separately from GMT. Second, we have made some minor changes to a few islands that have shown to be offset with respect to modern data (Tahiti, Moorea, and Mehetia in South Pacific and Agalega Islands in the Indian Ocean). Third, we have removed the item known as “Sandy Island” in the Corel Sea since available satellite data show no evidence of land in this area. Finally, we have purged 48 duplicates of very small islands (mostly in the Red Sea, Persian Gulf, and the Cook-Austral archipelago) where inaccurate WDBII versions and accurate WVS versions of the same features had survived our initial processing. This GMT release thus coincides with the latest GSHHG release version 2.2.2. We also fixed a typo in the MJD date in the gmtdefaults.txt documentation and updated the CM4 coefficient files for the mgd77 supplement.
Below is the list of bug corrections for individual programs:
]
Another bug-fix release, except for the mgd77 supplement where we now have added support for the new MGD77T tab-delimited format introduced by NGDC, and for ps2raster.c under Windows where the Ghostscript executable path is now fetched from the registry. In case that fails, we fall back to the old “get it from the path” mechanism. One common bug shared by several programs was the failure to consult FIELD_DELIMITER and/or D_FORMAT for ASCII output formatting. Below is the list of bug corrections for individual programs:
]
This is another bug-fix release, including an update to GSHHS (now at 2.2.0) which fixes a truncation error for the polygon areas (which only affected users of the gshhs supplement and not GMT itself). The supplement tool gshhs now has a few more options to allow better feature extraction for GSHHS users. Below is the list of bug corrections:
]
This is another bug-fix release, including an update to GSHHS which fixes error is the Germany-Poland border and a few “spiky” islands. Therefore, this version requires the new GSHHS 2.1.1 release. We also patched some errors in the “jet” color table. Below is the list of bug corrections:
]
This is again mostly a bug-fix release, and coincides with the availability of GMT 5.0.0. Due to a few issues we had an aborted update to 4.5.4 that was never announced; hence the 4.5.5-numbered version. A few minor improvements have been added:
Here is the list of bug corrections:
Also, Appendix F had missing shading for two items in the Standard+ table, and example 23 placed the city names at an angle of 1 degree rather than horizontally.
]
A few minor technical issues in the distribution led us to make a few changes and increment the version to 4.5.5.
]
This is mostly a bug-fix release, including more corrections to the political boundaries distributed via the GSHHS netCDF files (these affect the Syria-Israel, Israel-Jordan, Moldova-Ukraine, and the Eritrea-Ethiopia borders) as well as missing river-lake metadata in the GSHHS distribution. Therefore, this version requires the new GSHHS 2.1.0 release.
Here is the list of bug corrections:
Here is a list of the recent enhancement to various programs; these were introduced to correct mistakes or overcome limitations:
]
This is mostly another bug-fix release, including one that required us to add more meta-data to the GSHHS coastline netCDF files. Therefore, this version requires GSHHS 2.0.2 or higher. As was the case for 4.5.1, note that the GSHHS polygons themselves have not changed (still at version 2.0). We also added in the relatively recent Nunavut province boundary in Canada. However, some enhancements were added as well, most notably a new graph frame mode for linear projections (to add arrow heads to math axes) and a new symbol in psxy.c (to draw a circular arrow used to indicate angles); these capabilities are demonstrated in a new (and final) example 30. Finally, we fixed the long-standing problem of psxy -SE requesting major and minor axes but actually treating them as if they are semi-axes. We now consistently expect and use major and minor axes; you may thus notice a scaling of two if you continue to give semi-major/minor axes. Here is the list of bug corrections:
Here is a list of the recent enhancement to various programs:
]
This is almost entirely a bug-fix release where we address several 64-bit incompatibilities and rebuild the netcdf GSHHS library to include some attributes from GSHHS that were needed by new options in pscoast and other programs. Note that the GSHHS polygons have not changed (still at version 2.0), but we had to update the derived netcdf repackaging used by GMT to 2.0.1. However, some enhancements were added as well, most significantly support for the polyconic projection (-JPoly), experimental support for grid and image imports via GDAL (requires –enable-gdal during configure and properly installed GDAL libraries and include files), and allowing -JXwidth/height to recompute a height given as zero based on the width (or vice versa) and the aspect ratio of the region.
Here is the list of bug corrections:
Here is a list of the recent enhancement to various programs:
]
This is another significant update of the official distribution and hence it has a mix of bug fixes and program enhancements. We have added a new supplement (sph) which offers interpolation, triangulation (Delaunay and Voronoi), and distance calculations on a spherical surface. The hard work is done by the original effort of Robert Renka who developed the Fortran-77 SSRFPACK and STRIPACK libraries; these are here supplied via a f2c-assisted translation. The imgsrc supplement has a new Bourne script img2google, which simplifies making Google Earth tiles from Sandwell and Smith bathymetry. The mgd77 supplement has a new program mgd77magref, which is used to evaluate either the CM4 comprehensive geomagnetic model, a more sophisticated alternative to IGRF, or the IGRF. The misc supplement has received two new tools (gmt2kml and kml2gmt) that simplify the presentation of GMT data in Google Earth, and one (dimfilter) that offers directional spatial filtering of grids. The x2sys supplement has a new tool (x2sys_merge) to merge updated COEs table into a main COE table database. Finally, ps2raster.c has evolved further and can now be used to create simple KML files for Google Earth.
A major new enhancement is the global option -g, which is used to determine if excessive spacing between data points (“gaps”, to be defined in a variety of ways) should be used to segment an otherwise continuous line. We expect to enable -g in several programs during the next revision; at the moment it is available in gmtconvert, mapproject, psxy and psxyz. Given that all the lower-case GMT options deal with low-level data i/o settings we have decided to rename the -M option (which controls the presence of multiple segment headers) to -m; this allows us to promote this ubiquitous option to global status (i.e., has the same meaning in all GMT programs). Use of -M will remain valid for the rest of GMT 4.x but results in a warning about the new usage. Related to this is the introduction of a new parameter (NAN_RECORDS) that determines if NaNs in key columns (such as longitude, latitude) should constitute a line break or bad data to be skipped.
We have revised how ellipsoids are specified. When importing an ellipsoid file, we allow a,b,f as ellipsoid parameters, where b or f could both be zero. If file does not exist, attempt to read name as a[/[b=f=]f], meaning semi-major axis, b=semi-minor axis, f = flattening, or inverse flattening. We have also added parameters for the TOPEX ellipsoid and for the Moon and planets (IAU2000).
This release of GMT coincides with the release 2.0 of GSHHS, the coastline data used by GMT. In addition to general improvements to the data, we have expanded the -A option that controls the limits on what features to extract. New modifiers allow users to exclude “river-lakes” and any feature whose area is less than a fraction of the original full resolution feature.
Finally, our configure script continues to evolve and now better supports installation on 64-bit systems and can automatically detect if and where netCDF exists on your system.
Here is the list of bug corrections:
Here is a list of the recent enhancement to various programs:
]
This is a significant update of the official distribution and hence has a mix of bug fixes and program enhancements. We have added a new program (greenspline.c) which offers interpolation and gridding in 1–3 dimensions using Green’s functions of various splines. Also, the misc supplement has a new tool (gmtdp.c) which offers line-reduction using the Douglas-Peucker algorithm we used for the various shoreline resolutions. The mex supplement has a new Matlab/Octave function (imgread.m) to directly read Sandwell/Smith *.img files. The x2sys supplement has three new programs: x2sys_list.c can extract a subset of crossovers from the list produced by x2sys_cross.c, x2sys_report.c reports statistics of crossovers, whereas x2sys_solve.c will determine systematic trends from a set of crossover errors. These programs are intended to replace the old x_system tools x_list.c, x_report.c and x_solve_dc_drift.c. We have also temporarily added GMT_qsort which is a 64-bit compliant version of qsort. The latter is broken under OS X 64-bit and is thus substituted on that platform only for 64-bit compilations until Apple fixes the problem. Finally ps2raster.c can now be used to create geotiff images if gdal is installed on your system. Here is the list of bug corrections:
Here is a list of the recent enhancement to various programs:
Finally, we have added three new examples to demonstrate plotting of *.img grids, mixing UTM grids and geographic projections, and using greenspline.c for gridding on a spherical surface.
]
This quick update only 2 weeks after the release of version 4.3.0 was prompted by the discovery of three serious bugs; two of which were quite old but had caused no harm until tested under Fedora 9. The third critical bug prevented the wholesale reading and writing of GRD98 format grids. In addition a few minor bugs were discovered; this is the list of all corrections:
In addition, many of the supplements did not work properly under Windows due to internal problems with the DLL. Finally, one enhancement snuck in before the decision to issue this update was made:
]
Changes are once again a mix of structural improvements, bug fixes, and a few enhancements. The coastline files (now GSHHS 1.10) have seen minor modifications, the mex supplement now offers support for Octave, all source code is now fully 64-bit compliant, we have added an isolation mode option (if GMT_TMPDIR is defined, write temporary and hidden files to that directory), and the configure/make setup has been further improved (such as honoring CFLAGS and LDFLAGS set by user). Colors may now be specified as hexadecimal codes (e.g., #ff0000 for red), and projections can be specified by name (similar to Proj4). Finally, binary table data can now be COARDS-compliant netCDF files. As for documentation, we have now switched from C shell to Bourne shell (although the csh examples are still distributed).
The following lists specific enhancements or new program options:
A long list of bugs has been squashed since the last release, the most important are listed below:
]
Changes in GMT 4.2.1 once again address many structural issues as well as numerous bug fixes. System-wide changes include a revamping of the entire configure/make setup for both regular installations and CVS users, an improvement to how the BCR 2-D interpolations for images and grids are done by adding B-spline and nearest neighbor as optional interpolants, introduction of a new PostScript Level 2 pattern machinery in pslib.c, an updated GSHHS coastline version (which also includes Australia internal state boundaries, fixes to the Yemeni and Lebanese borders, and more river lines), and general improvements and corrections to the documentation, such as placing all man pages in section 1 (except pslib which goes in section 3). Starting with GMT 4.2.1 we will also begin naming the coastline-related archives by the GSHHS prefix and use the actual GSHHS version number (now 1.9).
Individual programs have also seen some new options or enhancements:
Below is a list of previous problems that we have identified and corrected in the current release:
A few bug-fixes applies to the supplements as well:
Finally, as far as CVS users are concerned, the old "gurumake" system has gone. To compile from CVS, users need to use a GNU compatible make program. A combination of GNUmakefile and makefile files make sure that those components not in the tarballs are created from scratch. Type make in the GMT directory to get a list of targets.
]
Changes in GMT 4.2.0 address many structural issues as well as many bug fixes. We have consolidated user initialization files in the /.gmt directory, continued to replace tiling with bitmaps, and have performed a myriad of under-the-hood changes. One imporant and more visible new feature is the fact that grdimage and pscoast now can use the general perspective projection with arbitrary elevation (-JG has been enhanced to handle the extra arguments required – see the new example 26 for details). Also, the coastline files have been updated to use GSHHS version 1.5 which fixes minor inconsistencies in the coastline database. We have also corrected issues that made the Windows DLL explode in 4.1.4. Finally, a few enhancements have been made to these programs:
Below is a list of previous problems that we have identified and corrected in the current release:
A few bug-fixes applies to the supplements as well:
]
Changes in GMT 4.1.4 are again relatively minor and predominantly bug fixes. One imporant new feature is that GMT can now automatically recognize the format of the grid file given to a program. The use of the “=id” mechanism is now only needed when writing an output file in a grid format other than the netCDF default or when reading using custom scaline and translation is required. We have also added a new user directory pointed to by GMT_USERDIR (default directory is /.gmt) where items such as .gmtdefaults4 will be looked for. Additionally, a few enhancements have been made to overcome limitations in the previous versions:
Below is a list of previous problems (a few accidently introduced in GMT 4.1.3) that we have identified and corrected in the current release:
A few bug-fixes applies to the supplements as well:
]
Changes in GMT 4.1.3 are relatively minor and predominantly bug fixes. However, a few enhancements have been made to overcome limitations in the previous versions:
Below is a list of previous problems (some accidentily introduced in GMT 4.1.2) that we have identified and corrected in the current release:
]
On the surface, changes in GMT 4.1.2 are relatively minor. Most of the work has involved realignment of code and parsing of arguments to simplify the upcoming port to GMT 5; a brief listing of more visible changes include
A few programs or options have received minor updates and new features, such as
Below is a list of previous problems that we have identified and corrected in the current release:
A few bug-fixes applies to the supplements as well:
]
Changes in GMT 4.1.1 are mostly minor; a brief listing include
A few programs or options have received minor updates and new features, such as
Inevitably, when new features are added, new bugs come along with them. Below is a list of problems that we have identified and corrected in the current release:.
A few bug-fixes applies to the supplements as well:
]
Most changes in GMT 4.1 are improvements “under the hood”. The most significant of these are
Many programs or options have received minor updates and new features, such as
Inevitably, when new features are added, new bugs come along with them. Below is a list of problems that we have identified and corrected in the current release:
]
GMT 4 represents a major overhaul of the package, hence the major version number increment. There are four categories of changes that have been implemented:
Most scientists are familiar with the sequence: raw data processing final illustration. In order to finalize papers for submission to scientific journals, prepare proposals, and create overheads and slides for various presentations, many scientists spend large amounts of time and money to create camera-ready figures. This process can be tedious and is often done manually, since available commercial or in-house software usually can do only part of the job. To expedite this process we introduce the Generic Mapping Tools (GMT for short), which is a free4, software package that can be used to manipulate columns of tabular data, time-series, and gridded data sets, and display these data in a variety of forms ranging from simple x-y plots to maps and color, perspective, and shaded-relief illustrations. GMT uses the PostScript page description language [Adobe Systems Inc., 1990]. With PostScript, multiple plot files can easily be superimposed to create arbitrarily complex images in gray tones or 24-bit true color. Line drawings, bitmapped images, and text can be easily combined in one illustration. PostScript plot files are device-independent: The same file can be printed at 300 dots per inch (dpi) on an ordinary laserwriter or at 2470 dpi on a phototypesetter when ultimate quality is needed. GMT software is written as a set of UNIX tools5 and is totally self-contained and fully documented. The system is offered free of charge and is distributed over the computer network (Internet) [Wessel and Smith, 1991; 1995a,b; 1998].
The original version 1.0 of GMT was released in the summer of 1988 when the authors were graduate students at Lamont-Doherty Earth Observatory of Columbia University. During our tenure as graduate students, L-DEO changed its computing environment to a distributed network of UNIX workstations, and we wrote GMT to run in this environment. It became a success at L-DEO, and soon spread to numerous other institutions in the US, Canada, Europe, and Japan. The current version benefits from the many suggestions contributed by users of the earlier versions, and now includes more than 50 tools, more than 30 projections, and many other new, more flexible features. GMT provides scientists with a variety of tools for data manipulation and display, including routines to sample, filter, compute spectral estimates, and determine trends in time series, grid or triangulate arbitrarily spaced data, perform mathematical operations (including filtering) on 2-D data sets both in the space and frequency domain, sample surfaces along arbitrary tracks or onto a new grid, calculate volumes, and find trend surfaces. The plotting programs will let the user make linear, log, and x–y diagrams, polar and rectangular histograms, maps with filled continents and coastlines choosing from many common map projections, contour plots, mesh plots, monochrome or color images, and artificially illuminated shaded-relief and 3-D perspective illustrations.
GMT is written in the highly portable ANSI C programming language [Kernighan and Ritchie, 1988], is fully POSIX compliant [Lewine, 1991], has no Year 2000 problems, and may be used with any hardware running some flavor of UNIX, possibly with minor modifications. In writing GMT, we have followed the modular design philosophy of UNIX: The raw data processing final illustration flow is broken down to a series of elementary steps; each step is accomplished by a separate GMT or UNIX tool. This modular approach brings several benefits: (1) only a few programs are needed, (2) each program is small and easy to update and maintain, (3) each step is independent of the previous step and the data type and can therefore be used in a variety of applications, and (4) the programs can be chained together in shell scripts or with pipes, thereby creating a process tailored to do a user-specific task. The decoupling of the data retrieval step from the subsequent massage and plotting is particularly important, since each institution will typically have its own data base formats. To use GMT with custom data bases, one has only to write a data extraction tool which will put out data in a form readable by GMT (discussed below). After writing the extractor, all other GMT modules will work as they are.
GMT makes full use of the PostScript page description language, and can produce color illustrations if a color PostScript device is available. One does not necessarily have to have access to a top-of-the-line color printer to take advantage of the color capabilities offered by GMT: Several companies offer imaging services where the customer provides a PostScript plot file and gets color slides or hardcopies in return. Furthermore, general-purpose PostScript raster image processors (RIPs) are now becoming available, letting the user create raster images from PostScript and plot these bitmaps on raster devices like computer screens, dot-matrix printers, large format raster plotters, and film writers6. Because the publication costs of color illustrations are high, GMT offers 90 common bit and hachure patterns, including many geologic map symbol types, as well as complete graytone shading operations. Additional bit and hachure patterns may also be designed by the user. With these tools, it is possible to generate publication-ready monochrome originals on a common laserwriter.
GMT is thoroughly documented and comes with a technical reference and cookbook which explains the purpose of the package and its many features, and provides numerous examples to help new users quickly become familiar with the operation and philosophy of the system. The cookbook contains the shell scripts that were used for each example; PostScript files of each illustration are also provided. All programs have individual manual pages which can be installed as part of the on-line documentation under the UNIX man utility or as web pages. In addition, the programs offer friendly help messages which make them essentially self-teaching – if a user enters invalid or ambiguous command arguments, the program will print a warning to the screen with a synopsis of the valid arguments. All the documentation is available for web browsing and may be installed at the user’s site.
The processing and display routines within GMT are completely general and will handle any (x,y) or (x,y,z) data as input. For many purposes the (x,y) coordinates will be (longitude, latitude) but in most cases they could equally well be any other variables (e.g., wavelength, power spectral density). Since the GMT plot tools will map these (x,y) coordinates to positions on a plot or map using a variety of transformations (linear, log-log, and several map projections), they can be used with any data that are given by two or three coordinates. In order to simplify and standardize input and output, GMT uses two file formats only. Arbitrary sequences of (x,y) or (x,y,z) data are read from multi-column ASCII tables, i.e., each file consists of several records, in which each coordinate is confined to a separate column7. This format is straightforward and allows the user to perform almost any simple (or complicated) reformatting or processing task using standard UNIX utilities such as cut, paste, grep, sed and awk. Two-dimensional data that have been sampled on an equidistant grid are read and written by GMT in a binary grid file using the functions provided with the netCDF library (a free, public-domain software library available separately from UCAR, the University Corporation of Atmospheric Research [Treinish and Gough, 1987]). This XDR (External Data Representation) based format is architecture independent, which allows the user to transfer the binary data files from one computer system to another8. GMT contains programs that will read ASCII (x,y,z) files and produce grid files. One such program, surface, includes new modifications to the gridding algorithm developed by Smith and Wessel [1990] using continuous splines in tension.
Most of the programs will produce some form of output, which falls into four categories. Several of the programs may produce more than one of these types of output:
GMT is available over the Internet at no charge. To obtain a copy, visit the GMT home page gmt.soest.hawaii.edu, and select DOWNLOAD. This page contains all the information you need to obtain GMT for your platform. We also maintain two electronic mailing lists you may subscribe to in order to stay informed about bug fixes and upgrades (See Chapter 7).
For those without net-access that need to obtain GMT: Geoware (http://www.geoware-online.com) makes and distributes CD-R and DVD-R media with the GMT package, compatible supplements, and several Gb of useful Earth and ocean science data sets. For more information send e-mail to geoware@geoware-online.com.
GMT has served a multitude of scientists very well, and their responses have prompted us to develop these programs even further. It is our hope that the new version will satisfy these users and attract new users as well. We present this system to the community in order to promote sharing of research software among investigators in the US and abroad.
The following is a summary of all the programs supplied with GMT and a very short description of their purpose. For more details, see the individual UNIX manual pages or the online web documentation. For a listing sorted by program purpose, see Section 3.2.
blockmean | L (x,y,z) table data filter/decimator |
blockmedian | L (x,y,z) table data filter/decimator |
blockmode | Mode estimate (x,y,z) table data filter/decimator |
filter1d | Filter 1-D table data sets (time series) |
fitcircle | Finds the best-fitting great or small circle for a set of points |
gmt2rgb | Convert Sun raster or grid file to red, green, blue component grids |
gmtconvert | Convert data tables from one format to another |
gmtdefaults | List the current default settings |
gmtmath | Mathematical operations on table data |
gmtselect | Select subsets of table data based on multiple spatial criteria |
gmtset | Change selected parameters in current .gmtdefaults4 file |
grd2cpt | Make color palette table from a 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 sets |
grdcontour | Contouring of 2-D gridded data sets |
grdcut | Cut a sub-region from a grid file |
grdedit | Modify header information in a 2-D grid file |
grdfft | Perform operations on grid files in the frequency domain |
grdfilter | Filter 2-D gridded data sets in the space domain |
grdgradient | Compute directional gradient from grid files |
grdhisteq | Histogram equalization for grid files |
grdimage | Produce images from 2-D gridded data sets |
grdinfo | Get information about grid files |
grdlandmask | Create masking grid files from shoreline data base |
grdmask | Reset grid nodes in/outside a clip path to constants |
grdmath | Mathematical operations on grid files |
grdpaste | Paste together grid files along a common edge |
grdproject | Project gridded data sets onto a new coordinate system |
grdreformat | Converts grid files into other grid formats |
grdsample | Resample a 2-D gridded data set onto a new grid |
grdtrack | Sampling of 2-D gridded data set along 1-D track |
grdtrend | Fits polynomial trends to grid files |
grdvector | Plotting of 2-D gridded vector fields |
grdview | 3-D perspective imaging of 2-D gridded data sets |
grdvolume | Calculate volumes under a surface within specified contour |
greenspline | Interpolation using Green’s functions for splines in 1–3 dimensions |
makecpt | Make color palette tables |
mapproject | Transformation of coordinate systems for table data |
minmax | Report extreme values in table data files |
nearneighbor | Nearest-neighbor gridding scheme |
project | Project table data onto lines or 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 (and fill) coastlines, borders, and rivers on maps |
pscontour | Contour or image raw table data by triangulation |
pshistogram | Plot a histogram |
psimage | Plot Sun raster files on a map |
pslegend | Plot a legend on a map |
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 table data 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 series |
trend2d | Fits polynomial trends to series |
triangulate | Perform optimal Delauney triangulation and gridding |
xyz2grd | Convert an equidistant table xyz file to a 2-D grid file |
Instead of an alphabetical listing, this section contains a summary sorted by program purpose. Also included is a
quick summary of the standard command line options and a breakdown of the -J option for each of the over 30
projections available in GMT.
FILTERING OF 1-D AND 2-D DATA | |
blockmean | L estimate () data filters/decimators |
blockmedian | L estimate () data filters/decimators |
blockmode | Mode estimate () data filters/decimators |
filter1d | Filter 1-D data (time series) |
grdfilter | Filter 2-D data in space domain |
PLOTTING OF 1-D and 2-D DATA | |
grdcontour | Contouring of 2-D gridded data |
grdimage | Produce images from 2-D gridded data |
grdvector | Plot vector fields from 2-D gridded data |
grdview | 3-D perspective imaging of 2-D gridded data |
psbasemap | Create a basemap frame |
psclip | Use polygon files as clipping paths |
pscoast | Plot coastlines, filled continents, rivers, and political borders |
pscontour | Direct contouring or imaging of xyz data by triangulation |
pshistogram | Plot a histogram |
psimage | Plot Sun raster files on a map |
pslegend | Plot a legend on a map |
psmask | Create overlay to mask specified regions of a map |
psrose | Plot sector or rose diagrams |
psscale | Plot gray scale or color scale |
pstext | Plot text strings |
pswiggle | Draw anomalies along track |
psxy | Plot symbols, polygons, and lines in 2-D |
psxyz | Plot symbols, polygons, and lines in 3-D |
GRIDDING OF (X,Y,Z) TABLE DATA
| |
greenspline | Interpolation using Green’s functions for splines in 1–3 dimensions |
nearneighbor | Nearest-neighbor gridding scheme |
surface | Continuous curvature gridding algorithm |
triangulate | Perform optimal Delauney triangulation on xyz data |
SAMPLING OF 1-D AND 2-D DATA | |
grdsample | Resample a 2-D gridded data onto new grid |
grdtrack | Sampling of 2-D data along 1-D track |
sample1d | Resampling of 1-D data |
PROJECTION AND MAP-TRANSFORMATION
| |
grdproject | Transform gridded data to a new coordinate system |
mapproject | Transform table data to a new coordinate system |
project | Project data onto lines or great circles |
INFORMATION | |
gmtdefaults | List the current default settings |
gmtset | Command-line editing of parameters in the .gmtdefaults4 file |
grdinfo | Get information about the content of grid files |
minmax | Report extreme values in table data files |
MISCELLANEOUS | |
gmtmath | Reverse Polish Notation (RPN) calculator for table data |
makecpt | Create GMT color palette tables |
spectrum1d | Compute spectral estimates from time-series |
triangulate | Perform optimal Delauney triangulation on xyz data |
CONVERT OR EXTRACT SUBSETS OF DATA | |
gmt2rgb | Convert Sun raster or grid file to red, green, blue component grids |
gmtconvert | Convert table data from one format to another |
gmtselect | Select table data subsets based on multiple spatial criteria |
grd2xyz | Convert 2-D gridded data to table data |
grdcut | Cut a sub-region from a grid file |
grdblend | Blend several partially over-lapping grid files onto one grid |
grdpaste | Paste together grid files along common edge |
grdreformat | Convert from one grid format to another |
splitxyz | Split () table data into several segments |
xyz2grd | Convert table data to 2-D grid file |
DETERMINE TRENDS IN 1-D AND 2-D DATA
| |
fitcircle | Finds best-fitting great or small circles |
grdtrend | Fits polynomial trends to grid files () |
trend1d | Fits polynomial or Fourier trends to series |
trend2d | Fits polynomial trends to series |
OTHER OPERATIONS ON 2-D GRIDS
| |
grd2cpt | Make color palette table from grid file |
grdclip | Limit the –range in gridded data sets |
grdedit | Modify grid header information |
grdfft | Operate on grid files in frequency domain |
grdgradient | Compute directional gradients from grid files |
grdhisteq | Histogram equalization for grid files |
grdlandmask | Creates mask grid file from coastline database |
grdmask | Set grid nodes in/outside a clip path to constants |
grdmath | Reverse Polish Notation (RPN) calculator for grid files |
grdvolume | Calculate volume under a surface within a contour |
ps2raster | Crop and convert PostScript files to raster images, EPS and PDF |
STANDARDIZED COMMAND LINE OPTIONS WITH OLD PROJECTION CODES
| |
-B[ps]xinfo[/yinfo[/zinfo]][WESNZwesnz+][:.title:] Tickmarks. Each info is
| |
[t]stride[±phase][u][lp][:"label":][:="prefix":][:,"unit":], where l and p apply to axes only,
| |
and type t = {a, A, f, g, i, I}; unit u = {c, C, d, D, h, H, K, k, m, M, o, O, r, R, u, U, y, Y}
| |
The leading ps selects primary [Default] or secondary axis items
| |
-H[i][n_headers] | ASCII [input] tables have header record[s] |
-J (upper case for width, lower case for scale) | Map projection |
-JA/[/horizon]/width | Lambert azimuthal equal area |
-JB////width | Albers conic equal area |
-JC//width | Cassini cylindrical |
-JCyl_stere/[/[/]]width | Cylindrical stereographic |
-JD////width | Equidistant conic |
-JE/[/horizon]/width | Azimuthal equidistant |
-JF/[/horizon]/width | Azimuthal gnomonic |
-JG/[/horizon]/width | Azimuthal orthographic |
-JG//alt/azim/tilt/twist/W/H/width | General perspective |
-JH[/]width | Hammer equal area |
-JI[/]width | Sinusoidal equal area |
-JJ[/]width | Miller cylindrical |
-JKf[/]width | Eckert IV equal area |
-JKs[/]width | Eckert VI equal area |
-JL////width | Lambert conic conformal |
-JM[/[/]]width | Mercator cylindrical |
-JN[/]width | Robinson |
-JOa//azim/width | Oblique Mercator, 1: origin and azimuth |
-JOb////width | Oblique Mercator, 2: two points |
-JOc////width | Oblique Mercator, 3: origin and pole |
-JP[a]width[/origin] | Polar [azimuthal] () (or cylindrical) |
-JPoly[/[/]]width | (American) polyconic |
-JQ[/[/]]width | Equidistant cylindrical |
-JR[/]width | Winkel Tripel |
-JS//[/horizon]/width | General stereographic |
-JT/[/]width | Transverse Mercator |
-JUzone/width | Universal Transverse Mercator (UTM) |
-JV[/]width | Van der Grinten |
-JW[/]width | Mollweide |
-JXwidth[lpexpTt][/height[lpexpTt]][d] | Linear, log, –, and time |
-JY//width | Cylindrical equal area |
-K | Append more PS later |
-O | This is an overlay plot |
-P | Select Portrait orientation |
-Rwest/east/south/north[/zmin/zmax][r] | Specify Region of interest |
-U[[just]/dx/dy/][label] | Plot time-stamp on plot |
-V | Run in verbose mode |
-X[acr]off [u] | Shift plot origin in -direction |
-Y[acr]off [u] | Shift plot origin in -direction |
-b[io][csSdD][ncol] | Select binary input or output |
-ccopies | Set number of plot copies [1] |
-f[io]colinfo | Set formatting of ASCII input or output |
-g[+]xXyYdDgap[u] | Segment data by detecting gaps |
-m[io]flag | Set multi-segment data mode |
-:[io] | Expect y/x input rather than x/y |
-B[ps]xinfo[/yinfo[/zinfo]][WESNZwesnz+][:.title:] Tickmarks. Each info is
| |
[t]stride[±phase][u][lp][:"label":][:="prefix":][:,"unit":], where l and p apply to axes only,
| |
and type t = {a, A, f, g, i, I}; unit u = {c, C, d, D, h, H, K, k, m, M, o, O, r, R, u, U, y, Y}
| |
The leading ps selects primary [Default] or secondary axis items
| |
-H[i][n_headers] | ASCII [input] tables have header record[s] |
-J (upper case for width, lower case for scale) | Map projection |
-Jaea/////scale | Albers conic equal area |
-Jaeqd//[/horizon]/scale | Azimuthal equidistant |
-Jcass///scale | Cassini cylindrical |
-Jcea///scale | Cylindrical equal area |
-Jcyl_stere/[/[/]]scale | Cylindrical stereographic |
-Jeqc/[/[/]]scale | Equidistant cylindrical |
-Jeqdc/////scale | Equidistant conic |
-Jgnom//[/horizon]/scale | Azimuthal gnomonic |
-Jhammer/[/]scale | Hammer equal area |
-Jeck4/[/]scale | Eckert IV equal area |
-Jeck6/[/]scale | Eckert VI equal area |
-Jlaea//[/horizon]/scale | Lambert azimuthal equal area |
-Jlcc/////scale | Lambert conic conformal |
-Jmerc/[/[/]]scale | Mercator cylindrical |
-Jmill/[/]scale | Miller cylindrical |
-Jmoll/[/]scale | Mollweide |
-Jnsper///alt/azim/tilt/twist/W/H/scale | General perspective |
-Jomerc///azim/scale | Oblique Mercator, 1: origin and azimuth |
-Jomerc/////scale | Oblique Mercator, 2: two points |
-Jomercp/////scale | Oblique Mercator, 3: origin and pole |
-Jortho//[/horizon]/scale | Azimuthal orthographic |
-Jpolar/[a]scale[/origin] | Polar [azimuthal] () (or cylindrical) |
-Jpoly[/[/]]width | (American) polyconic |
-Jrobin/[/]scale | Robinson |
-Jsinu/[/]scale | Sinusoidal equal area |
-Jstere///[/horizon]/scale | General stereographic |
-Jtmerc//[/]scale | Transverse Mercator |
-Jutm/zone/scale | Universal Transverse Mercator (UTM) |
-Jvandg/[/]scale | Van der Grinten |
-Jwintri/[/]scale | Winkel Tripel |
-Jxy/xscale[lpexpTt][/yscale[lpexpTt]][d] | Linear, log, –, and time |
-K | Append more PS later |
-O | This is an overlay plot |
-P | Select Portrait orientation |
-Rwest/east/south/north[/zmin/zmax][r] | Specify Region of interest |
-U[[just]/dx/dy/][label] | Plot time-stamp on plot |
-V | Run in verbose mode |
-X[acr]off [u] | Shift plot origin in -direction |
-Y[acr]off [u] | Shift plot origin in -direction |
-b[io][csSdD][ncol] | Select binary input or output |
-ccopies | Set number of plot copies [1] |
-f[io]colinfo | Set formatting of ASCII input or output |
-g[+]xXyYdDgap[u] | Segment data by detecting gaps |
-m[io]flag | Set multi-segment data mode |
-:[io] | Expect y/x input rather than x/y |
This section explains features common to all the programs in GMT and summarizes the philosophy behind the system. Some of the features described here may make more sense once you reach the cook-book section where we present actual examples of their use.
GMT programs can accept dimensional quantities in cm, inch, meter, or point (1/72 of an inch)9. There are two ways to ensure that GMT understands which unit you intend to use.
The latter method is less secure as other users may have a different unit set and your script may not work as intended. We therefore recommend you always supply the desired unit explicitly.
There are about 100 parameters which can be adjusted individually to modify the appearance of plots or affect the manipulation of data. When a program is run, it initializes all parameters to the GMT defaults10, then tries to open the file .gmtdefaults4 in the current directory11. If not found, it will look for that file in a sub-directory /̃.gmt of your home directory, and finally in your home directory itself. If successful, the program will read the contents and set the default values to those provided in the file. By editing this file you can affect features such as pen thicknesses used for maps, fonts and font sizes used for annotations and labels, color of the pens, dots-per-inch resolution of the hardcopy device, what type of spline interpolant to use, and many other choices (A complete list of all the parameters and their default values can be found in the gmtdefaults manual pages). Figures 4.2.1, 4.2, and 4.3 show the parameters that affect plots). You may create your own .gmtdefaults4 files by running gmtdefaults and then modify those parameters you want to change. If you want to use the parameter settings in another file you can do so by specifying +<defaultfile> on the command line. This makes it easy to maintain several distinct parameter settings, corresponding perhaps to the unique styles required by different journals or simply reflecting font changes necessary to make readable overheads and slides. Note that any arguments given on the command line (see below) will take precedent over the default values. E.g., if your .gmtdefaults4 file has x offset = 1i as default, the -X1.5i option will override the default and set the offset to 1.5 inches.
There are at least two good reasons why the GMT default options are placed in a separate parameter file:
As mentioned, GMT programs will attempt to open a file named .gmtdefaults4. At times it may be desirable to override that default. There are several ways in which this can be accomplished.
gmtset ANNOT_FONT_PRIMARY Times-Bold ANNOT_FONT_SIZE_PRIMARY 12
These changes will remain in effect until they are overridden.
In addition to those parameters that directly affect the plot there are numerous parameters than modify units, scales, etc. For a complete listing, see the gmtdefaults man pages. We suggest that you go through all the available parameters at least once so that you know what is available to change via one of the described mechanisms.
Note: All examples presented in this document started by copying the file .gmtdefaults4.doc from the directory doc/scripts to .gmtdefaults4. As a result the commands gmtset of other scripts were overall, reverting to a “virgin” of parameters set in .gmtdefaults4.doc. The graphs in Chapter 7 were created using .gmtdefaults4.doc from the directory examples after which the graphs were scaled down by 50%.
Each program requires certain arguments specific to its operation. These are explained in the manual pages and in the usage messages. Most programs are “case-sensitive”; almost all options must start with an upper-case letter. We have tried to choose letters of the alphabet which stand for the argument so that they will be easy to remember. Each argument specification begins with a hyphen (except input file names; see below), followed by a letter, and sometimes a number or character string immediately after the letter. Do not space between the hyphen, letter, and number or string. Do space between options. Example:
pscoast -R0/20/0/20 -Ggray -JM6i -Wthin -B5 -V -P map.ps
Most of the programs take many of the same arguments like those related to setting the data region, the map projection, etc. The 17 switches in Table 4.1 have the same meaning in all the programs (although some programs may not use all of them). These options will be described here as well as in the manual pages, as is vital that you understand how to use these options. We will present these options in order of importance (some are use a lot more than others).
Option | Meaning
|
-B | Defines tickmarks, annotations, and labels for basemaps and axes |
-H | Specifies that input/output tables have header record(s) |
-J | Selects a map projection or coordinate transformation |
-K | Allows more plot code to be appended to this plot later |
-O | Allows this plot code to be appended to an existing plot |
-P | Selects Portrait plot orientation [Default is landscape] |
-R | Defines the extent of the map/plot region |
-U | Plots a time-stamp, by default in the lower left corner of page |
-V | Selects verbose operation; reporting on progress |
-X | Sets the x-coordinate for the plot origin on the page |
-Y | Sets the y-coordinate for the plot origin on the page |
-b | Selects binary input and/or output |
-c | Specifies the number of plot copies |
-f | Specifies the data format on a per column basis |
-g | Identify data gaps based on supplied criteria |
-m | Specifies data in multiple segment format |
-: | Assumes input geographic data are (lat,lon) and not (lon,lat) |
The -R option defines the map region or data domain of interest. It may be specified in one of three ways (Figure 4.4):
For rectilinear projections the first two forms give identical results. Depending on the selected map projection (or the kind of expected input data), the boundary coordinates may take on three different formats:
When used in conjunction with the Cartesian Linear Transformation (-Jx or -JX) —which can be used to map floating point data, geographical coordinates, as well as time coordinates— it is prudent to indicate that you are using geographical coordinates in one of the following ways:
The optional clock string is a 24-hour clock in hh[:mm[:ss[.xxx]]] format. If no clock is given it implies 00:00:00, i.e., the start of the specified day. Note that not all of the specified entities need be present in the data. All calendar date-clock strings are internally represented as double precision seconds since proleptic Gregorian date Mon Jan 1 00:00:00 0001. Proleptic means we assume that the modern calendar can be extrapolated forward and backward; a year zero is used, and Gregory’s reforms12 are extrapolated backward. Note that this is not historical.
This option selects the coordinate transformation or map projection. The general format is
Since GMT version 4.3.0, there is an alternative way to specify the projections: use the same abbreviation as in the mapping package Proj4. The options thus either look like:
The projections available in GMT are presented in Figure 4.5. For details on all GMT projections and the required parameters, see the psbasemap man page. We will also show examples of every projection in the next Chapters, and a quick summary of projection syntax was given in Chapter 3.
This is by far the most complicated option in GMT, but most examples of its usage are actually quite simple. Given as
-B[ps]xinfo[/yinfo[/zinfo]][:."title
string":][Ww][Ee][Ss][Nn][Zz[+]],
this switch specifies map boundaries (or plot axes) to be plotted by using the selected information. The
optional flag following -B selects p(rimary) [Default] or s(econdary) axes information (mostly used for
time axes annotations; see examples below). The components xinfo, yinfo and zinfo are of the form
info[:"axis label":][:="prefix":][:,"unit label":]
where info is one or more concatenated substrings of the form [t]stride[±phase][u]. The t flag sets the axis item of interest; the available items are listed in Table 4.2. By default, all 4 map boundaries (or plot axes) are plotted (denoted W, E, S, N). To change this selection, append the codes for those you want (e.g., WSn). Upper case (e.g., W) will annotate in addition to draw axis/tick-marks. The title, if given, will appear centered above the plot14. Unit label or prefix may start with a leading – to suppress the space between it and the annotation. Normally, equidistant annotations occur at multiples of stride; you can phase-shift this by appending ±phase.
Note that the appearance of certain time annotations (month-, week-, and day-names) may be affected by the TIME_LANGUAGE, TIME_FORMAT_PRIMARY, and TIME_FORMAT_SECONDARY settings.
The unit flag u can take on one of 18 codes; these are listed in Table 4.3. Almost all of these units are time-axis specific. However, the m and c units will be interpreted as arc minutes and arc seconds, respectively, when a map projection is in effect.
Flag | Unit | Description |
Y | year | Plot using all 4 digits |
y | year | Plot using last 2 digits |
O | month | Format annotation using PLOT_DATE_FORMAT |
o | month | Plot as 2-digit integer (1–12) |
U | ISO week | Format annotation using PLOT_DATE_FORMAT |
u | ISO week | Plot as 2-digit integer (1–53) |
r | Gregorian week | 7-day stride from start of week (see TIME_WEEK_START) |
K | ISO weekday | Plot name of weekday in selected language |
k | weekday | Plot number of day in the week (1-7) (see TIME_WEEK_START) |
D | date | Format annotation using PLOT_DATE_FORMAT |
d | day | Plot day of month (1–31) or day of year (1–366) |
(see PLOT_DATE_FORMAT | ||
R | day | Same as d; annotations aligned with week (see TIME_WEEK_START) |
H | hour | Format annotation using PLOT_CLOCK_FORMAT |
h | hour | Plot as 2-digit integer (0–24) |
M | minute | Format annotation using PLOT_CLOCK_FORMAT |
m | minute | Plot as 2-digit integer (0–60) |
C | seconds | Format annotation using PLOT_CLOCK_FORMAT |
c | seconds | Plot as 2-digit integer (0–60) |
There may be two levels of annotations. Here, “primary” refers to the annotation that is closest to the axis (this is the primary annotation), while “secondary” refers to the secondary annotation that is plotted further from the axis. The examples below will clarify what is meant. Note that the terms “primary” and “secondary” do not reflect any hierarchical order of units: The “primary” annotation interval is usually smaller (e.g., days) while the “secondary” annotation interval typically is larger (e.g., months).
Geographic basemaps may differ from regular plot axis in that some projections support a “fancy” form of axis and is selected by the BASEMAP_TYPE setting. The annotations will be formatted according to the PLOT_DEGREE_FORMAT template and DEGREE_SYMBOL setting. A simple example of part of a basemap is shown in Figure 4.6.
The machinery for primary and secondary annotations introduced for time-series axes can also be utilized for geographic basemaps. This may be used to separate degree annotations from minutes- and seconds-annotations. For a more complicated basemap example using several sets of intervals, including different intervals and pen attributes for grid lines and grid crosses, see Figure 4.7.
For non-geographic axes, the BASEMAP_TYPE setting is implicitly set to plain. Other than that, cartesian linear axes are very similar to geographic axes. The annotation format may be controlled with the D_FORMAT parameter. By default, it is set to “%g”, which is a C language format statement for floating point numbers15, and with this setting the various axis routines will automatically determine how many decimal points should be used by inspecting the stride settings. If D_FORMAT is set to another format it will be used directly (.e.g, “%.2f” for a fixed, two decimals format). Note that for these axes you may use the unit setting to add a unit string to each annotation (see Figure 4.8).
Due to the logarithmic nature of annotation spacings, the stride parameter takes on specific meanings. The following concerns are specific to log axes:
Normally, stride will be used to create equidistant (in the user’s unit) annotations or ticks, but because of the exponential nature of the axis, such annotations may converge on each other at one end of the axis. To avoid this problem, you can append p to stride, and the annotation interval is expected to be in transformed units, yet the annotation itself will be plotted as un-transformed units. E.g., if stride = 1 and power = 0.5 (i.e., sqrt), then equidistant annotations labeled 1, 4, 9, ... will appear.
What sets time axis apart from the other kinds of plot axes is the numerous ways in which we may want to tick and annotate the axis. Not only do we have both primary and secondary annotation items but we also have interval annotations versus tickmark annotations, numerous time units, and several ways in which to modify the plot. We will demonstrate this flexibility with a series of examples. While all our examples will only show a single -axis, time-axis is supported for all axes.
Our first example shows a time period of almost two months in Spring 2000. We want to annotate the month intervals as well as the date at the start of each week:
_________________________________________________________________________________
These commands result in Figure 4.11. Note the leading hyphen in the PLOT_DATE_FORMAT removes leading zeros from calendar items (e.g., 02 becomes 2).
The next example shows two different ways to annotate an axis portraying 2 days in July 1969:
_________________________________________________________________________________
The lower example (Figure 4.12) chooses to annotate the weekdays (by specifying a1K) while the upper example choses dates (by specifying a1D). Note how the clock format only selects hours and minutes (no seconds) and the date format selects a month name, followed by one space and a two-digit day-of-month number.
The third example presents two years, annotating both the years and every 3rd month.
_________________________________________________________________________________
Note that while the year annotation is centered on the 1-year interval, the month annotations must be centered on the corresponding month and not the 3-month interval. The PLOT_DATE_FORMAT selects month name only and TIME_FORMAT_PRIMARY selects the 1-character, upper case abbreviation of month names using the current language (selected by TIME_LANGUAGE).
The fourth example (Figure 4.14) only shows a few hours of a day, using relative time by specifying t in the -R option while the TIME_UNIT is d (for days). We select both primary and secondary annotations, ask for a 12-hour clock, and let time go from right to left:
_________________________________________________________________________________
The fifth example shows a few weeks of time (Figure 4.15). The lower axis shows ISO weeks with week numbers and abbreviated names of the weekdays. The upper uses Gregorian weeks (which start at the day chosen by TIME_WEEK_START); they do not have numbers. ______________________________________________________________________________
Our sixth example shows the first five months of 1996, and we have annotated each month with an abbreviated, upper case name and 2-digit year. Only the primary axes information is specified.
Our seventh and final example illustrates annotation of year-days. Unless we specify the formatting with a leading hyphen in PLOT_DATE_FORMAT we get 3-digit integer days. Note that in order to have the two years annotated we need to allow for the annotation of small fractional intervals; normally such truncated interval must be at least half of a full interval. ______________________________________________________________________________________________________________________
The -H[i][n_recs] option lets GMT know that input file(s) have one [Default] or more header records. If there are more than one header record you must specify the number after the -H option, e.g., -H4. The default number of header records if -H is used is one of the many parameters in the .gmtdefaults4 file (N_HEADER_RECS), but can be overridden by -Hn_header_recs. Note that blank lines and records that be start with the character # are automatically skipped. Normally, programs that both read and write tables will output the header records that are found on input. Use -Hi to suppress the writing of header records.
-P selects Portrait plotting mode16. In general, a plot has an x-axis increasing from left to right and a y-axis increasing from bottom to top. If the paper is turned so that the long dimension of the paper is parallel to the x-axis then the plot is said to have Landscape orientation. If the long dimension of the paper parallels the y-axis the orientation is called Portrait (think of taking pictures with a camera and these words make sense). The default Landscape orientation is obtained by translating the origin in the x-direction (by the width of the chosen paper PAPER_MEDIA) and then rotating the coordinate system counterclockwise by 90°. By default the PAPER_MEDIA is set to Letter (or A4 if SI is chosen); this value must be changed when using different media, such as 11" x 17" or large format plotters (Figure 4.18).
The -K and -O options control the generation of PostScript code for multiple overlay plots. All PostScript files must have a header (for initializations), a body (drawing the figure), and a trailer (printing it out) (see Figure 4.19). Thus, when overlaying several GMT plots we must make sure that the first plot call omits the trailer, that all intermediate calls omit both header and trailer, and that the final overlay omits the header. -K omits the trailer which implies that more PostScript code will be appended later [Default terminates the plot system]. -O selects Overlay plot mode and omits the header information [Default initializes a new plot system]. Most unexpected results for multiple overlay plots can be traced to the incorrect use of these options. If you run only one plot program, ignore both the -O and -K options; they are only used when stacking plots.
-U draws UNIX System time stamp. Optionally, append an arbitrary text string (surrounded by double quotes), or the code c, which will plot the current command string (Figure 4.20).
-V selects verbose mode, which will send progress reports to stderr [Default runs “silently”]. The interpretation of this option can be toggled by changing the default VERBOSE.
-X and -Y shift origin of plot by (xoff ,yoff ) inches (Default is (X_ORIGIN, Y_ORIGIN) for new plots17 and (0,0) for overlays (-O)). By default, all translations are relative to the previous origin (see Figure 4.21). Supply offset as c to center the plot in that direction relative to the page margin. Absolute translations (i.e., relative to a fixed point (0,0) at the lower left corner of the paper) can be achieve by prepending “a” to the offsets. Subsequent overlays will be co-registered with the previous plot unless the origin is shifted using these options. The offsets are measured in the current coordinates system (which can be rotated using the initial -P option; subsequent -P options for overlays are ignored).
All GMT programs that accept table data input may read ASCII, native binary, or netCDF data. When using native binary data the user must be aware of the fact that GMT has no way of determining the actual number of columns in the file. You must therefore pass that information to GMT via the binary -bi[s]n option, where n is the actual number of data columns (s indicates single (4 bytes) rather than double (8 bytes) precision). If uppercase S (or D) are used it implies that byte-swapping should be performed just prior to writing (for output) or immediately after reading (for input). Note that n may be larger than m, the number of columns that the GMT program requires to do its task. If n is not given then it defaults to m. If n m an error is generated.
Because of its meta data, reading netCDF tables (i.e., netCDF files containing 1-dimensional arrays) is quite a bit less complex than reading native binary files. When feeding netCDF tables to programs like psxy, the program will automatically recognize the format and read whatever amount of columns are needed for that program. To steer which columns are to be read, the user can either append the suffix ?var1/var2/... to the netCDF file name or add the option -bicvar1/var2/..., where var1, var2, etc. are the names of the variables to be processed. The latter option is particularly practical when more than one file is read: the -bic option will apply to all files.
Currently, netCDF tables can only be input, not output. For more information, see Appendix B.
The -c option specifies the number of plot copies [Default is 1]. This value is embedded in the PostScript file and will make a printer issue the chosen number of copies without respooling.
When map projections are not required we must explicitly state what kind of data each input or output column contains. This is accomplished with the -f option. Following an optional i (for input only) or o (for output only), we append a text string with information about each column (or range of columns) separated by commas. Each string starts with the column number (0 is first column) followed by either x (longitude), y (latitude), T (absolute calendar time) or t (relative time). If several consecutive columns have the same format you may specify a range of columns rather than a single column, i.e., 0-4 for the first 5 columns. For example, if our input file has geographic coordinates (latitude, longitude) with absolute calendar coordinates in the columns 3 and 4, we would specify fi0y,1x,3-4T. All other columns are assumed to have the default, floating point format and need not be set individually. The shorthand -f[io]g means -f[io]0x,1y (geographic coordinates). For more information, see Sections 4.10 and 4.11.
GMT has several mechanisms that can determine line segmentation. Typically, data segments are separated by multiple segment header records (see section 4.4.14 on -m below). However, if key data columns contain a NaN we may also use that information to break lines into multiple segments. This behavior is modified by the parameter NAN_RECORDS which by default is set to skip, meaning such records are considered bad and simply skipped. If you wish such records to indicate a segment boundary then set this parameter to pass. Finally, you may wish to indicate gaps based on the data values themselves. The -g option is used to detect gaps based on one or more criteria (use -g+ if all the criteria must be met; otherwise only one of the specified criteria needs to be met to signify a data gap). Gaps can be based on excessive jumps in the x- or y-coordinates (-gx or -gy), or on the distance between points (-gd). Append the gap distance and optionally a unit for actual distances. For geographic data the optional unit may be meter [Default], kilometer, miles, or nautical miles. For programs that maps data to map coordinates you can optionally specify these criteria to apply to the projected coordinates (by using upper-case -gX, -gY or-gD). In that case, choose from inch, centimeter, meter, or points. [Default unit is controlled by MEASURE_UNIT]. Note: For -gx or -gy with time data the unit is instead controlled by TIME_UNIT.
The -m option states that the input and output table data will contain special records that marks the start of new line (or polygon) segments. These records are identified by their first character, which can be specified as an argument to -m [The default is ]. Append the modifiers i or o if the option should only apply to input or output, respectively. For binary data a multiple segment header is identified as a data record where all fields equal NaN. See section 4.4.14 and Appendix B for more details.
For geographical data, the first column is expected to contain longitudes and the second to contain latitudes. To reverse this expectation you must apply the -: option. Optionally, append i or o to restrict the effect to input or output only. Note that command line arguments that may take geographic coordinates (e.g., -R) always expect longitude before latitude. Also, geographical grids are expected to have the longitude as first (minor) dimension.
GMT programs “remember” the standardized command line options (See Section 4.4) given during their previous invocations and this provides a shorthand notation for complex options. For example, if a basemap was created with an oblique Mercator projection, specified as
-Joc170W/25:30S/33W/56:20N/1:500000
then a subsequent psxy command to plot symbols only needs to state -Jo in order to activate the same projection. In contrast, note that -J by itself will pick the most recently used projection. Previous commands are maintained in the file .gmtcommands4, of which there will be one in each directory you run the programs from. This is handy if you create separate directories for separate projects since chances are that data manipulations and plotting for each project will share many of the same options. Note that an option spelled out on the command line will always override the last entry in the .gmtcommands4 file and, if execution is successful, will replace this entry as the previous option argument in the .gmtcommands4 file. If you call several GMT modules piped together then GMT cannot guarantee that the .gmtcommands4 file is processed in the intended order from left to right. The only guarantee is that the file will not be clobbered since GMT uses advisory file locking. The uncertainty in processing order makes the use of shorthands in pipes unreliable. We therefore recommend that you only use shorthands in single process command lines, and spell out the full command option when using chains of commands connected with pipes.
Each program carries a usage message. If you enter the program name without any arguments, the program will write the complete usage message to standard error (your screen, unless you redirect it). This message explains in detail what all the valid arguments are. If you enter the program name followed by a hyphen (–) only you will get a shorter version which only shows the command line syntax and no detailed explanations. If you incorrectly specify an option or omit a required option, the program will produce syntax errors and explain what the correct syntax for these options should be. If an error occurs during the running of a program, the program will in some cases recognize this and give you an error message. Usually this will also terminate the run. The error messages generally begin with the name of the program in which the error occurred; if you have several programs piped together this tells you where the trouble is.
Most of the programs which expect table data input can read either standard input or input in one or several files. These programs will try to read stdin unless you type the filename(s) on the command line without the above hyphens. (If the program sees a hyphen, it reads the next character as an instruction; if an argument begins without a hyphen, it tries to open this argument as a filename). This feature allows you to connect programs with pipes if you like. If your input is ASCII and has one or more header records, you must use the -H option (see Section 4.4.4). For binary table data no headers are allowed. ASCII files may in many cases also contain sub-headers separating data segments. These are called “multi-segment files” and requires a special option (typically -m); see Appendix B for complete documentation.
If filenames are given for reading, GMT programs will first look for them in the current directory. If the file is not found, the programs will look in two other directories pointed to by environmental parameters (if set). These are GMT_DATADIR and GMT_USERDIR, and they may be set by the user to point to directories that contain data sets of general use. Normally, the first directory (or directories: add multiple paths by separating them with colons (semi-colons under Windows)) will hold data sets of a general nature (tables, grids), although a particular use is to make available large grids accessible via the supplemental programs grdraster or img2grd; see Appendix A for information about these supplemental programs. The second directory may hold miscellaneous data sets more specific to the user; this directory also stores GMT defaults and other configuration files. Data sets that the user finds are often needed may be placed in these directories, thus eliminating the need to specify a full path to the file. Program output is always written to the current directory unless a full path has been specified.
Most of the programs take an optional -V argument which will run the program in the “verbose” mode (see Section 4.4.8). Verbose will write to standard error information about the progress of the operation you are running. Verbose reports things such as counts of points read, names of data files processed, convergence of iterative solutions, and the like. Since these messages are written to stderr, the verbose talk remains separate from your data output.
Most programs write their results, including PostScript plots, to standard output. The exceptions are those which may create binary netCDF grid files such as surface (due to the design of netCDF a filename must be provided; however, alternative binary output formats allowing piping are available; see Section 4.17). With UNIX you can redirect standard output to a file or pipe it into another process. Error messages, usage messages, and verbose comments are written to standard error in all cases. You can use UNIX to redirect standard error as well, if you want to create a log file of what you are doing.
Most of the time, GMT will know what kind of and coordinates it is reading because you have selected a particular coordinate transformation or map projection. However, there may be times when you must explicitly specify what you are providing as input using the -f switch. When binary input data are expected (-bi) they must all be floating point numbers, however for ASCII input there are numerous ways to encode data coordinates (which may be separated by white-space or commas). Valid input data are generally of the same form as the arguments to the -R option (see Section 4.4.1), with additional flexibility for calendar data. Geographical coordinates, for example, can be given in decimal degrees (e.g., -123.45417) or in the [±]ddd[:mm[:ss[.xxx]]][WESN] format (e.g., 123:27:15W).
Because of the widespread use of incompatible and ambiguous formats, the processing of input date components is guided by the template INPUT_DATE_FORMAT in your .gmtdefaults4 file; it is by default set to yyyy-mm-dd. Y2K-challenged input data such as 29/05/89 can be processed by setting INPUT_DATE_FORMAT to dd/mm/yy. A complete description of possible formats is given in the gmtdefaults man page. The clock string is more standardized but issues like 12- or 24-hour clocks complicate matters as well as the presence or absence of delimiters between fields. Thus, the processing of input clock coordinates is guided by the template INPUT_CLOCK_FORMAT which defaults to hh:mm:ss.xxx.
GMT programs that require a map projection argument will implicitly know what kind of data to expect, and the input processing is done accordingly. However, some programs that simply report on minimum and maximum values or just do a reformatting of the data will in general not know what to expect, and furthermore there is no way for the programs to know what kind of data other columns (beyond the leading and columns) contain. In such instances we must explicitly tell GMT that we are feeding it data in the specific geographic or calendar formats (floating point data are assumed by default). We specify the data type via the -f option (which sets both input and output formats; use -fi and -fo to set input and output separately). For instance, to specify that the the first two columns are longitude and latitude, and that the third column (e.g., ) is absolute calendar time, we add -fi0x,1y,2T to the command line. For more details, see the man page for the program you need to use.
The numerical output from GMT programs can be binary (when -bo is used) or ASCII [Default]. In the latter case the issue of formatting becomes important. GMT provides extensive machinery for allowing just about any imaginable format to be used on output. Analogous to the processing of input data, several templates guide the formatting process. These are OUTPUT_DATE_FORMAT and OUTPUT_CLOCK_FORMAT for calendar-time coordinates, OUTPUT_DEGREE_FORMAT for geographical coordinates, and D_FORMAT for generic floating point data. In addition, the user have control over how columns are separated via the FIELD_SEPARATOR parameter. Thus, as an example, it is possible to create limited FORTRAN-style card records by setting D_FORMAT to %7.3lf and FIELD_SEPARATOR to none [Default is tab].
PostScript is a command language for driving graphics devices such as laser printers. It is ASCII text which you can read and edit as you wish (assuming you have some knowledge of the syntax). We prefer this to binary metafile plot systems since such files cannot easily be modified after they have been created. GMT programs also write many comments to the plot file which make it easier for users to orient themselves should they need to edit the file (e.g., % Start of x-axis). All GMT programs create PostScript code by calling the pslib plot library (The user may call these functions from his/her own C or FORTRAN plot programs. See the manual pages for pslib syntax). Although GMT programs can create very individualized plot code, there will always be cases not covered by these programs. Some knowledge of PostScript will enable the user to add such features directly into the plot file. By default, GMT will produce freeform PostScript output with embedded printer directives. To produce Encapsulated PostScript (EPS) that can be imported into graphics programs such as IslandDraw, CorelDraw, Illustrator or Freehand for further embellishment, change the PAPER_MEDIA setting in the .gmtdefaults4 file. See Appendix C and the gmtdefaults man page for more details.
A pen in GMT has three attributes: width, color, and texture. Most programs will accept pen attributes in the form of an option argument, with commas separating the given attributes, e.g.,
-W[width[cipm]],[color],[texture[cipm]]
Table 4.5 contains additional examples of pen specifications suitable for, say, psxy.
Pen example | Comment
|
-W0.5p | Solid black line, 0.5 point thick |
-Wgreen | Solid green line with default width |
-Wthin,red,- | Dashed, thin red line |
-Wfat,. | Fat dotted line [black] |
-W0.1c,120-1-1 | Green (in h-s-v) pen, 1 mm thick |
-Wfaint,100/0/0/0,..- | Very thin, cyan (in c/m/y/k), dot-dot-dashed line |
Many plotting programs will allow the user to draw filled polygons or symbols. The fill specification may take two forms:
-Gfill
-Gpdpi/pattern[:Bcolor[Fcolor]]
Due to PostScript implementation limitations the raster images used with -G must be less than 146 x 146 pixels in size; for larger images see psimage. The format of Sun raster files is outlined in Appendix B. Note that under PostScript Level 1 the patterns are filled by using the polygon as a clip path. Complex clip paths may require more memory than the PostScript interpreter has been assigned. There is therefore the possibility that some PostScript interpreters (especially those supplied with older laserwriters) will run out of memory and abort. Should that occur we recommend that you use a regular grayshade fill instead of the patterns. Installing more memory in your printer may or may not solve the problem!
Table 4.6 contains a few examples of fill specifications.
Fill example | Comment
|
-G128 | Solid gray |
-G127/255/0 | Chartreuse, R/G/B-style |
-G#00ff00 | Green, hexadecimal RGB code |
-G25-0.86-0.82 | Chocolate, h-s-v – style |
-GDarkOliveGreen1 | One of the named colors |
-Gp300/7 | Simple diagonal hachure pattern in b/w at 300 dpi |
-Gp300/7:Bred | Same, but with red lines on white |
-Gp300/7:BredF- | Now the gaps between red lines are transparent |
-Gp100/marble.ras | Using user image of marble as the fill at 100 dpi |
Several programs, such as those which read 2-D gridded data sets and create colored images or shaded reliefs, need to be told what colors to use and over what z-range each color applies. This is the purpose of the color palette table (cpt-file). These files may also be used by psxy and psxyz to plot color-filled symbols. For most applications, you will simply create a cpt-file using the tool makecpt which will take an existing color table and resample it to fit your chosen data range, or use grd2cpt to build a cpt-file based on the data distribution in one or more given grid files. However, in some situations you will need to make a cpt-file by hand or using text tools like awk or perl.
Color palette tables (CPT) comes in two flavors: (1) Those designed to work with categorical data (e.g., data where interpolation of values is undefined) and (2) those designed for regular, continuously-varying data.
Note: This is an experimental component and is only available if you compile GMT with -DGMT_CPT2. Categorical data are information on which normal numerical operations are not defined. As an example, consider various land classifications (desert, forest, glacier, etc.) and it is clear that even if we assigned a numerical value to these categories (e.g., desert = 1, forest = 2, etc) it would be meaningless to compute average values (what would 1.5 mean?). For such data a special format of the CPT files are provided. Here, each category is assigned a unique key, a color or pattern, and an optional label (usually the category name). Keys must be monotonically increasing but do not need to be consecutive. The format is
key | fill | label |
… | ||
key | fill | label |
The fill information follows the format given in Section 4.14. While not always applicable to categorical data, the background color (for key-values key), foreground color (for key-values key), and not-a-number (NaN) color (for key-values = NaN) are all defined in the .gmtdefaults4 file, but can be overridden by the statements
B | R | G | B |
F | R | G | B |
N | R | G | B |
Here, the colors may be specified either in the RGB- (red, green, blue), CMYK- (cyan, magenta, yellow, black), or in the HSV-system (hue, saturation, value, and here the comment # COLOR_MODEL = HSV must be present in the cpt file since there are no other way to distinguish between HSV and RGB). Color names can also be used. Using the RGB system20, the format of the cpt-file is:
z | R | G | B | z | R | G | B | [A] | [;label] |
… | |||||||||
z | R | G | B | z | R | G | B | [A] | [;label] |
Thus, for each “z-slice”, defined as the interval between two boundaries (e.g., z to z), the color can be constant (by letting R = R, G = G, and B = B) or a continuous, linear function of z. The optional flag A is used to indicate annotation of the color scale when plotted using psscale. The optional code A may be L, U, or B to select annotation of the lower, upper, or both limits of the particular -slice. However, the standard -B option can be used by psscale to affect annotation and ticking of color scales. The optional semicolon followed by a text label will make psscale, when used with the -L option, place the supplied label instead of formatted z-values.
As for categorical tables, the background color (for z-values z), foreground color (for z-values z), and not-a-number (NaN) color (for z-values = NaN) are all defined in the .gmtdefaults4 file, but can be overridden by the statements
B | R | G | B |
F | R | G | B |
N | R | G | B |
which can be inserted into the beginning or end of the cpt-file. If you prefer the HSV system, set the .gmtdefaults4 parameter accordingly and replace red, green, blue with hue, saturation, value. Color palette tables that contain grayshades only may replace the r/g/b triplets with a single grayshade in the 0–255 range. For CMYK, give four values in the 0–100 range. Both the min and max color specifications in one z-slice must use the same color system, i.e., you cannot mix “red” and 0/255/100 on the same line.
A few programs (i.e., those that plot polygons such as grdview, psscale, and psxy) can accept pattern fills instead of grayshades. You must specify the pattern as in Section 4.14 (no leading -G of course), and only the first (low ) is used (we cannot interpolate between patterns). Finally, some programs let you skip features whose -slice in the cptfile has grayshades set to –. As an example, consider
30 | p200/16 | 80 | – | ||||
80 | – | 100 | – | ||||
100 | 200 | 0 | 0 | 200 | 255 | 255 | 0 |
200 | yellow | 300 | green | ||||
where slice is painted with pattern # 16 at 200 dpi, slice is skipped, slice is painted in a range of dark red to yellow, whereas the slice will linearly yield colors from yellow to green, depending on the actual value of .
Some programs like grdimage and grdview apply artificial illumination to achieve shaded relief maps. This is typically done by finding the directional gradient in the direction of the artificial light source and scaling the gradients to have approximately a normal distribution on the interval [-1,+1]. These intensities are used to add “white” or “black” to the color as defined by the z-values and the cpt-file. An intensity of zero leaves the color unchanged. Higher values will brighten the color, lower values will darken it, all without changing the original hue of the color (see Appendix I for more details). The illumination is decoupled from the data grid file in that a separate grid file holding intensities in the [-1,+1] range must be provided. Such intensity files can be derived from the data grid using grdgradient and modified with grdhisteq, but could equally well be a separate data set. E.g., some side-scan sonar systems collect both bathymetry and backscatter intensities, and one may want to use the latter information to specify the illumination of the colors defined by the former. Similarly, one could portray magnetic anomalies superimposed on topography by using the former for colors and the latter for shading.
For annotation labels or text strings plotted with pstext, GMT provides several escape sequences that allow the user to temporarily switch to the symbol font, turn on sub- or superscript, etc., within words. These conditions are toggled on/off by the escape sequence @x, where x can be one of several types. The escape sequences recognized in GMT are listed in Table 4.7. Only one level of sub- or superscript is supported. Note that under Windows the percent symbol indicates a batch variable, hence you must use two percent-signs for each one required in the escape sequence for font switching.
Code | Effect
|
@̃ | Turns symbol font on or off |
@+ | Turns superscript on or off |
@- | Turns subscript on or off |
@# | Turns small caps on or off |
@_ | Turns underline on or off |
@%fontno% | Switches to another font; @%% resets to previous font |
@:size: | Switches to another font size; @:: resets to previous size |
@;color; | Switches to another font color; @;; resets to previous color |
@! | Creates one composite character of the next two characters |
@@ | Prints the @ sign itself |
Shorthand notation for a few special European characters has also been added (Table 4.8):
PostScript fonts used in GMT may be re-encoded to include several accented characters used in many European languages. To access these, you must specify the full octal code xxx allowed for your choice of character encodings determined by the CHAR_ENCODING setting described in the gmtdefaults man page. Only the special characters belonging to a particular encoding will be available. Many characters not directly available by using single octal codes may be constructed with the composite character mechanism @!.
Some examples of escape sequences and embedded octal codes in GMT strings using the Standard+ encoding:
2@~p@~r@+2@+h@-0@- E\363tv\363s | = | 2 Eötvös |
10@+-3 @Angstr@om | = | 10 Ångstrøm |
Se@nor Gar@con | = | Señor Garçon |
M@!\305anoa stra@se | = | M | anoa straße |
A@\#cceleration@\# (ms@+-2@+) | = | ACCELERATION (MS) |
The option in pstext to draw a rectangle surrounding the text will not work for strings with escape sequences. A chart of characters and their octal codes is given in Appendix F.
GMT has the ability to read and write grids using more than one grid file format (see Table 4.9 for supported format and their IDs). For reading, GMT will automatically determine the format of grid files, while for writing you will normally have to append =ID to the filename if you want GMT to use a different format than the default.
By default, GMT will create new grid files using the nf format; however, this behavior can be overridden by setting the GRID_FORMAT defaults parameter to any of the other recognized values (or by appending =ID).
GMT can also read netCDF grid files produced by other software packages, provided the grid files satisfy the COARDS and Hadley Centre conventions for netCDF grids. Thus, products created under those conventions (provided the grid is 2-, 3-, 4-, or 5-dimensional) can be read directly by GMT and the netCDF grids written by GMT can be read by other programs that conform to those conventions. Three such programs are ncview, Panoply and ncBrowse; others can be found on the netCDF website.
In addition, users with some C-programming experience may add their own read/write functions and link them with the GMT library to extend the number of predefined formats. Technical information on this topic can be found in the source file gmt_customio.c. Users who are considering this approach should contact the GMT team.
ID | GMT 4 netCDF standard formats
|
nb | GMT netCDF format (byte) (COARDS-compliant) |
ns | GMT netCDF format (short) (COARDS-compliant) |
ni | GMT netCDF format (int) (COARDS-compliant) |
nf | GMT netCDF format (float) (COARDS-compliant) |
nd | GMT netCDF format (double) (COARDS-compliant) |
ID | GMT 3 netCDF legacy formats
|
cb | GMT netCDF format (byte) (depreciated) |
cs | GMT netCDF format (short) (depreciated) |
ci | GMT netCDF format (int) (depreciated) |
cf | GMT netCDF format (float) (depreciated) |
cd | GMT netCDF format (double) (depreciated) |
ID | GMT native binary formats
|
bm | GMT native, C-binary format (bit-mask) |
bb | GMT native, C-binary format (byte) |
bs | GMT native, C-binary format (short) |
bi | GMT native, C-binary format (int) |
bf | GMT native, C-binary format (float) |
bd | GMT native, C-binary format (double) |
ID | Miscellaneous grid formats
|
rb | SUN raster file format (8-bit standard) |
rf | GEODAS grid format GRD98 (NGDC) |
sf | Golden Software Surfer format 6 (float) |
sd | Golden Software Surfer format 7 (double) |
af | Atlantic Geoscience Center AGC (float) |
gd | Read-only via GDAL [EXPERIMENTAL] (float) |
Because some formats have limitations on the range of values they can store it is sometimes necessary to provide more than simply the name of the file and its ID on the command line. For instance, a native short integer file may use a unique value to signify an empty node or NaN, and the data may need translation and scaling prior to use. Therefore, all GMT programs that read or write grid files will decode the given filename as follows:
name[=ID[/scale/offset[/nan]]]
where everything in brackets is optional. If you are reading a grid then no options are needed: just continue to pass the name of the grid file. However, if you write another format you must append the =ID string, where ID is the format code listed above. In addition, should you want to (1) multiply the data by a scale factor, and (2) add a constant offset you must append the /scale/offset modifier. Finally, if you need to indicate that a certain data value should be interpreted as a NaN (not-a-number) you must append the /nan suffix to the scaling string (it cannot go by itself; note the nesting of the brackets!).
Some of the grid formats allow writing to standard output and reading from standard input which means you can connect GMT programs that operate on grid files with pipes, thereby speeding up execution and eliminating the need for large, intermediate grid files. You specify standard input/output by leaving out the filename entirely. That means the “filename” will begin with “=ID ” since no GMT netCDF format allow piping (due to the design of netCDF).
Everything looks clearer after a few examples:
Programs that both read and/or write more than one grid file may specify different formats and/or scaling for the files involved. The only restriction with the embedded grid specification mechanism is that no grid files may actually use the “=” character as part of their name (presumably, a small sacrifice).
One can also define special file suffixes to imply a specific file format; this approach represents a more intuitive and user-friendly way to specify the various file formats. The user may create a file called .gmt_io in the current directory, home directory or in the directory /̃.gmt and define any number of custom formats. The following is an example of a .gmt_io file:
# GMT i/o shorthand file |
# It can have any number of comment lines like this one anywhere |
# suffix | format_id | scale | offset | NaN | Comments |
grd | nf | - | - | - | Default format |
b | bf | - | - | - | Native binary floats |
i2 | bs | - | - | 32767 | 2-byte integers with NaN value |
ras | rb | - | - | - | Sun raster files |
byte | bb | - | - | 255 | Native binary 1-byte grids |
bit | bm | - | - | - | Native binary 0 or 1 grids |
mask | bm | - | - | 0 | Native binary 1 or NaN masks |
faa | bs | 0.1 | - | 32767 | Native binary gravity in 0.1 mGal |
These suffices can be anything that makes sense to the user. To activate this mechanism, set parameter GRIDFILE_SHORTHAND to TRUE in your .gmtdefaults4 file. Then, using the filename stuff.i2 is equivalent to saying stuff.i2=bs/1/0/32767, and the filename wet.mask means wet.mask=bm/1/0/0. For a file intended for masking, i.e., the nodes are either 1 or NaN, the bit or mask format file may be as small as 1/32 the size of the corresponding grid float format file.
When the netCDF file contains more than one 2-dimensional variable, GMT programs will load the first such variable in the file and ignore all others. Alternatively, the user can select the required variable by adding the suffix “?varname” to the file name. For example, to get information on the variable “slp” in file file.nc, use:
Since COARDS-compliant netCDF files are the default, the additional suffix “=nf” can be omitted.
In case the named variable is 3-dimensional, GMT will load the first (bottom) layer. If another layer is required, either add “[index]” or “(level)”, where index is the index of the third (depth) variable (starting at 0 for the first layer) and level is the numerical value of the third (depth) variable associated with the requested layer. To indicate the second layer of the 3-D variable “slp” use as file name: file.nc?slp[1].
When you supply the numerical value for the third variable using “(level)”, GMT will pick the layer closest to that value. No interpolation is performed.
Note that the question mark, brackets and parentheses have special meanings on Unix-based platforms. Therefore, you will need to either escape these characters, by placing a backslash in front of them, or place the whole file name plus modifiers between single quotes or double quotes.
A similar approach is followed for loading 4-dimensional grids. Consider a 4-dimensional grid with the following variables:
To get information on the 1010 grid of pressure at depth 10 and at time 24, one would use:
or (only in case the coordinates increase linearly):
The COARDS conventions set restrictions on the names that can be used for the units of the variables and coordinates. For example, the units of longitude and latitude are “degrees_east” and “degrees_north”, respectively. Here is an example of the header of a COARDS compliant netCDF file (to be obtained using ncdump):
This file contains two grids, which can be plotted separately using the names M2_fes2004.nc?amp and M2_fes2004.nc?pha. The attributes long_name and unit for each variable are combined in GMT to a single unit string. For example, after reading the grid y_unit equals latitude [degrees_north]. The same method can be used in reverse to set the proper variable names and units when writing a grid. However, when the coordinates are set properly as geographical or time axes, GMT will take care of this. The user is, however, still responsible for setting the variable name and unit of the z-coordinate. The default is simply “z”.
For a variety of data processing and plotting tasks there is a need to acknowledge that a data point is missing or unassigned. In the “old days” such information was passed by letting a value like -9999.99 take on the special meaning of “this is not really a value, it is missing”. The problem with this scheme is that -9999.99 (or any other floating point value) may be a perfectly reasonable data value and in such a scenario would be skipped. The solution adopted in GMT is to use the IEEE concept Not-a-Number (NaN) for this purpose. Mathematically, a NaN is what you get if you do an undefined mathematical operation like ; in ASCII data files they appear as the textstring NaN. This value is internally stored with a particular bit pattern defined by IEEE so that special action can be taken when it is encountered by programs. In particular, a library function called isnan is used to test if a floating point is a NaN. GMT uses these tests extensively to determine if a value is suitable for plotting or processing (if a NaN is used in a calculation the result would become NaN as well). Data points whose values equal NaN are not normally plotted (or plotted with the special NaN color given in .gmtdefaults4). Several tools such as xyz2grd, gmtmath, and grdmath can convert user data to NaN and vice versa, thus facilitating arbitrary masking and clipping of data sets. Note that a few computers do not have native IEEE hardware support. At this point, this applies to some of the older Cray super-computers. Users on such machines may have to adopt the old ‘-9999.99” scheme to achieve the desired results.
Data records that contain NaN values for the x or y columns (or the z column for cases when 3-D Cartesian data are expected) are usually skipped during reading. However, the presence of these bad records can be interpreted in two different ways, and this behavior is controlled by the NAN_RECORDS defaults parameter. The default setting (gap) considers such records to indicate a gap in an otherwise continuous series of points (e.g., a line), and programs can act upon this information, e.g., not to draw a line across the gap or to break the line into separate segments. The alternative setting (bad) makes no such interpretation and simply reports back how many bad records were skipped during reading.
GMT relies on several environment parameters, in particular to find data files and program settings.
GMT programs read real-world coordinates and convert them to positions on a plot. This is achieved by selecting one of several coordinate transformations or projections. We distinguish between three sets of such conversions:
The next chapter will be dedicated to GMT map projections in its entirety. Meanwhile, the present chapter will summarize the properties of the Cartesian and Polar coordinate transformations available in GMT, list which parameters define them, and demonstrate how they are used to create simple plot axes. We will mostly be using psbasemap (and occasionally psxy) to demonstrate the various transformations. Our illustrations may differ from those you reproduce with the same commands because of different settings in our .gmtdefaults4 file.) Finally, note that while we will specify dimensions in inches (by appending i), you may want to use cm (c), meters (m), or points (p) as unit instead (see the gmtdefaults man page).
GMT Cartesian coordinate transformations come in three flavors:
These transformations convert input coordinates to locations on a plot. There is no coupling between and (i.e., and ); it is a one-dimensional projection. Hence, we may use separate transformations for the - and -axes (and -axes for 3-D plots). Below, we will use the expression , where is either or (or for 3-D plots). The coefficients in depend on the desired plot size (or scale), the chosen domain, and the nature of itself.
Two subsets of linear will be discussed separately; these are a polar (cylindrical) projection and a linear projection applied to geographic coordinates (with a 360° periodicity in the -coordinate). We will show examples of all of these projections using dummy data sets created with gmtmath, a “Reverse Polish Notation” (RPN) calculator that operates on or creates table data:
_________________________________________________________________________________
There are in fact three different uses of the Cartesian linear transformation, each associated with specific command line options. The different manifestations result from specific properties of three kinds of data:
Selection of the Cartesian linear transformation with regular floating point coordinates will result in a simple linear scaling of the input coordinates. The projection is defined by stating
If the y-scale or y-axis length is different from that of the x-axis (which is most often the case), separate the two scales (or lengths) by a slash, e.g., -Jx0.1i/0.5i or -JX8i/5i. Thus, our data sets will plot as shown in Figure 5.1.
The complete commands given to produce this plot were
_________________________________________________________________________________
Normally, the user’s x-values will increase to the right and the y-values will increase upwards. It should be noted that in many situations it is desirable to have the direction of positive coordinates be reversed. For example, when plotting depth on the y-axis it makes more sense to have the positive direction downwards. All that is required to reverse the sense of positive direction is to supply a negative scale (or axis length). Finally, sometimes it is convenient to specify the width (or height) of a map and let the other dimension be computed based on the implied scale and the range of the other axis. To do this, simply specify the length to be recomputed as 0.
While the Cartesian linear projection is primarily designed for regular floating point x,y data, it is sometimes necessary to plot geographical data in a linear projection. This poses a problem since longitudes have a 360° periodicity. GMT therefore needs to be informed that it has been given geographical coordinates even though a linear transformation has been chosen. We do so by adding a g (for geographical) or d (for degrees) directly after -R or by appending a g or d to the end of the -Jx (or -JX) option. As an example, we want to plot a crude world map centered on 125°E. Our command will be
_________________________________________________________________________________
with the result reproduced in Figure 5.2.
Several particular issues arise when we seek to make linear plots using calendar date/time as the input coordinates. As far as setting up the coordinate transformation we must indicate whether our input data have absolute time coordinates or relative time coordinates. For the former we append T after the axis scale (or width), while for the latter we append t at the end of the -Jx (or -JX) option. However, other command line arguments (like the -R option) may already specify whether the time coordinate is absolute or relative. An absolute time entry must be given as [date]T[clock] (with date given as yyyy[-mm[-dd]], yyyy[-jjj], or yyyy[-Www[-d]], and clock using the hh[:mm[:ss[.xxx]]] 24-hour clock format) whereas the relative time is simply given as the units of time since the epoch followed by t (see TIME_UNIT and TIME_EPOCH for information on specifying the time unit and the epoch). As a simple example, we will make a plot of a school week calendar (Figure 5.3).
_________________________________________________________________________________
When the coordinate ranges provided by the -R option and the projection type given by -JX (including the optional d, g, t or T) conflict, GMT will warn the users about it. In general, the options provided with -JX will prevail.
The log transformation is simply and is selected by appending an l (lower case L) immediately following the scale (or axis length) value. Hence, to produce a plot in which the x-axis is logarithmic (the y-axis remains linear, i.e., a semi-log plot), try
_________________________________________________________________________________
Note that if x- and y-scaling are different and a log-log plot is desired, the l must be appended twice: Once after the x-scale (before the /) and once after the y-scale.
This projection uses and allows us to explore exponential relationships like x versus y. While and can be any values, we will select and which means we will plot versus . We indicate this scaling by appending a p (lower case P) followed by the desired exponent, in our case 0.5. Since we do not need to specify p1 since it is identical to the linear transformation. Thus our command becomes
_________________________________________________________________________________
This transformation converts polar coordinates (angle and radius ) to positions on a plot. Now and , hence it is similar to a regular map projection because and are coupled and (i.e., ) has a 360° periodicity. With input and output points both in the plane it is a two-dimensional projection. The transformation comes in two flavors:
Consequently, the polar transformation is defined by providing
As an example of this projection we will create a gridded data set in polar coordinates using grdmath, a RPN calculator that operates on or creates grid files.
_________________________________________________________________________________
We used grdcontour to make a contour map of this data. Because the data file only contains values with , a donut shaped plot appears in Figure 5.6.
GMT implements more than 30 different projections. They all project the input coordinates longitude and latitude to positions on a map. In general, and , where is implicitly given as the radial vector length to the point on the chosen ellipsoid. The functions and can be quite nasty and we will refrain from presenting details in this document. The interested read is referred to Snyder [1987]21. We will mostly be using the pscoast command to demonstrate each of the projections. GMT map projections are grouped into four categories depending on the nature of the projection. The groups are
Because and are coupled we can only specify one plot-dimensional scale, typically a map scale (for lower-case map projection code) or a map width (for upper-case map projection code). However, in some cases it would be more practical to specify map height instead of width, while in other situations it would be nice to set either the shortest or longest map dimension. Users may select these alternatives by appending a character code to their map dimension. To specify map height, append h to the given dimension; to select the minimum map dimension, append -, whereas you may append + to select the maximum map dimension. Without the modifier the map width is selected by default.
In GMT version 4.3.0 we noticed we ran out of the alphabet for 1-letter (and sometimes 2-letter) projection codes. To allow more flexibility, and to make it easier to remember the codes, we implemented the option to use the abbreviations used by the Proj4 mapping package. Since some of the GMT projections are not in Proj4, we invented some of our own as well. For a full list of both the old 1- and 2-letter codes, as well as the Proj4-equivalents see the quick reference cards in Section 3.2. For example, -JM15c and -JMerc/15c have the same meaning.
This projection, developed by Albers in 1805, is predominantly used to map regions of large east-west extent, in particular the United States. It is a conic, equal-area projection, in which parallels are unequally spaced arcs of concentric circles, more closely spaced at the north and south edges of the map. Meridians, on the other hand, are equally spaced radii about a common center, and cut the parallels at right angles. Distortion in scale and shape vanishes along the two standard parallels. Between them, the scale along parallels is too small; beyond them it is too large. The opposite is true for the scale along meridians. To define the projection in GMT you need to provide the following information:
Note that you must include the “1:” if you choose to specify the scale that way. E.g., you can say 0.5 which means 0.5 inch/degree or 1:200000 which means 1 inch on the map equals 200,000 inches along the standard parallels. The projection center defines the origin of the rectangular map coordinates. As an example we will make a map of the region near Taiwan. We choose the center of the projection to be at 125 °E/20 °N and 25 °N and 45 °N as our two standard parallels. We desire a map that is 5 inches wide. The complete command needed to generate the map below is therefore given by:
_________________________________________________________________________________
The equidistant conic projection was described by the Greek philosopher Claudius Ptolemy about A.D. 150. It is neither conformal or equal-area, but serves as a compromise between them. The scale is true along all meridians and the standard parallels. To select this projection in GMT you must provide the same information as for the other conic projection, i.e.,
The equidistant conic projection is often used for atlases with maps of small countries. As an example, we generate a map of Cuba:
_________________________________________________________________________________
This conic projection was designed by the Alsatian mathematician Johann Heinrich Lambert (1772) and has been used extensively for mapping of regions with predominantly east-west orientation, just like the Albers projection. Unlike the Albers projection, Lambert’s conformal projection is not equal-area. The parallels are arcs of circles with a common origin, and meridians are the equally spaced radii of these circles. As with Albers projection, it is only the two standard parallels that are distortion-free. To select this projection in GMT you must provide the same information as for the Albers projection, i.e.,
The Lambert conformal projection has been used for basemaps for all the 48 contiguous States with the two fixed standard parallels 33°N and 45°N. We will generate a map of the continental USA using these parameters. Note that with all the projections you have the option of selecting a rectangular border rather than one defined by meridians and parallels. Here, we choose the regular WESN region, a “fancy” basemap frame, and use degrees west for longitudes. The generating commands used were
_________________________________________________________________________________
The choice for projection center does not affect the projection but it indicates which meridian (here 100°W) will be vertical on the map. The standard parallels were originally selected by Adams to provide a maximum scale error between latitudes 30.5°N and 47.5°N of 0.5–1%. Some areas, like Florida, experience scale errors of up to 2.5%.
The polyconic projection, in Europe usually referred to as the American polyconic projection, was introduced shortly before 1820 by the Swiss-American cartographer Ferdinand Rodulph Hassler (1770-1843). As head of the Survey of the Coast, he was looking for a projection that would give the least distortion for mapping the coast of the United States. The projection acquired its name from the construction of each parallel, which is achieved by projecting the parallel onto the cone while it is rolled around the globe, along the central meridian, tangent to that parallel. As a consequence, the projection involves many cones rather than a single one used in regular conic projections.
The polyconic projection is neither equal-area, nor conformal. It is true to scale without distortion along the central meridian. Each parallel is true to scale as well, but the meridians are not as they get further away from the central meridian. As a consequence, no parallel is standard because conformity is lost with the lengthening of the meridians.
Below we reproduce the illustration by Snyder [1987], with a gridline every 10° and annotations only every 30° in longitude:
_________________________________________________________________________________
This projection was developed by Lambert in 1772 and is typically used for mapping large regions like continents and hemispheres. It is an azimuthal, equal-area projection, but is not perspective. Distortion is zero at the center of the projection, and increases radially away from this point. To define this projection in GMT you must provide the following information:
Two different types of maps can be made with this projection depending on how the region is specified. We will give examples of both types.
In this mode we define our region by specifying the longitude/latitude of the lower left and upper right corners instead of the usual west, east, south, north boundaries. The reason for specifying our area this way is that for this and many other projections, lines of equal longitude and latitude are not straight lines and are thus poor choices for map boundaries. Instead we require that the map boundaries be rectangular by defining the corners of a rectangular map boundary. Using 0°E/40°S (lower left) and 60°E/10°S (upper right) as our corners we try
_________________________________________________________________________________
Note that an “r” is appended to the -R option to inform GMT that the region has been selected using the rectangle technique, otherwise it would try to decode the values as west, east, south, north and report an error since ’east’ ’west’.
Here, you must specify the world as your region (-Rg or -Rd). E.g., to obtain a hemisphere view that shows the Americas, try
_________________________________________________________________________________
To geologists, the Lambert azimuthal equal-area projection (with origin at 0°/0°) is known as the equal-area (Schmidt) stereonet and used for plotting fold axes, fault planes, and the like. An equal-angle (Wulff) stereonet can be obtained by using the stereographic projection (discussed later). The stereonets produced by these two projections appear below.
This is a conformal, azimuthal projection that dates back to the Greeks. Its main use is for mapping the polar regions. In the polar aspect all meridians are straight lines and parallels are arcs of circles. While this is the most common use it is possible to select any point as the center of projection. The requirements are
A default map scale factor of 0.9996 will be applied by default (although you may change this with MAP_SCALE_FACTOR). However, the setting is ignored when a standard parallel has been specified since the scale is then implicitly given. We will look at two different types of maps.
In our first example we will let the projection center be at the north pole. This means we have a polar stereographic projection and the map boundaries will coincide with lines of constant longitude and latitude. An example is given by
_________________________________________________________________________________
As with Lambert’s azimuthal equal-area projection we have the option to use rectangular boundaries rather than the wedge-shape typically associated with polar projections. This choice is defined by selecting two points as corners in the rectangle and appending an “r” to the -R option. This command produces a map as presented in Figure 6.9:
_________________________________________________________________________________
In terms of usage this projection is identical to the Lambert azimuthal equal-area projection. Thus, one can make both rectangular and hemispheric maps. Our example shows Australia using a projection pole at 130E/30°S. The command used was
_________________________________________________________________________________
By choosing 0°/0°as the pole, we obtain the conformal stereonet presented next to its equal-area cousin in the Section 6.2.1 on the Lambert azimuthal equal-area projection (Figure 6.7).
The perspective projection imitates in 2 dimensions the 3-dimensional view of the earth from space. The implementation in GMT is very flexible, and thus requires many input variables. Those are listed and explained below, with the values used in Figure 6.11 between brackets.
The imagined view of northwest Europe from a Space Shuttle at 230 km looking due east is thus accomplished by the following pscoast command:
_________________________________________________________________________________
The orthographic azimuthal projection is a perspective projection from infinite distance. It is therefore often used to give the appearance of a globe viewed from outer space. As with Lambert’s equal-area and the stereographic projection, only one hemisphere can be viewed at any time. The projection is neither equal-area nor conformal, and much distortion is introduced near the edge of the hemisphere. The directions from the center of projection are true. The projection was known to the Egyptians and Greeks more than 2,000 years ago. Because it is mainly used for pictorial views at a small scale, only the spherical form is necessary.
To specify the orthographic projection the same options -Jg or -JG as the perspective projection are used, but with fewer variables to supply:
Our example of a perspective view centered on 75°W/40°N can therefore be generated by the following pscoast command:
_________________________________________________________________________________
The most noticeable feature of this azimuthal projection is the fact that distances measured from the center are true. Therefore, a circle about the projection center defines the locus of points that are equally far away from the plot origin. Furthermore, directions from the center are also true. The projection, in the polar aspect, is at least several centuries old. It is a useful projection for a global view of locations at various or identical distance from a given point (the map center).
To specify the azimuthal equidistant projection you must supply:
Our example of a global view centered on 100°W/40°N can therefore be generated by the following pscoast command. Note that the antipodal point is 180° away from the center, but in this projection this point plots as the entire map perimeter:
_________________________________________________________________________________
The Gnomonic azimuthal projection is a perspective projection from the center onto a plane tangent to the surface. Its origin goes back to the old Greeks who used it for star maps almost 2500 years ago. The projection is neither equal-area nor conformal, and much distortion is introduced near the edge of the hemisphere; in fact, less than a hemisphere may be shown around a given center. The directions from the center of projection are true. Great circles project onto straight lines. Because it is mainly used for pictorial views at a small scale, only the spherical form is necessary.
To specify the Gnomonic projection you must supply:
Using a horizon of 60°, our example of this projection centered on 120°W/35°N can therefore be generated by the following pscoast command:
_________________________________________________________________________________
Cylindrical projections are easily recognized for its shape: maps are rectangular and meridians and parallels are straight lines crossing at right angles. But that is where similarities between the cylindrical projections supported by GMT (Mercator, transverse Mercator, universal transverse Mercator, oblique Mercator, Cassini, cylindrical equidistant, cylindrical equal-area, Miller, and cylindrical stereographic projections) stops. Each have a different way of spacing the meridians and parallels to obtain certain desirable cartographic properties.
Probably the most famous of the various map projections, the Mercator projection takes its name from the Flemish cartographer Gheert Cremer, better known as Gerardus Mercator, who presented it in 1569. The projection is a cylindrical and conformal, with no distortion along the equator. A major navigational feature of the projection is that a line of constant azimuth is straight. Such a line is called a rhumb line or loxodrome. Thus, to sail from one point to another one only had to connect the points with a straight line, determine the azimuth of the line, and keep this constant course for the entire voyage22. The Mercator projection has been used extensively for world maps in which the distortion towards the polar regions grows rather large, thus incorrectly giving the impression that, for example, Greenland is larger than South America. In reality, the latter is about eight times the size of Greenland. Also, the Former Soviet Union looks much bigger than Africa or South America. One may wonder whether this illusion has had any influence on U.S. foreign policy.
In the regular Mercator projection, the cylinder touches the globe along the equator. Other orientations like vertical and oblique give rise to the Transverse and Oblique Mercator projections, respectively. We will discuss these generalizations following the regular Mercator projection.
The regular Mercator projection requires a minimum of parameters. To use it in GMT programs you supply this information (the first two items are optional and have defaults):
Our example presents a world map at a scale of 0.012 inch pr degree which will give a map 4.32 inch wide. It was created with the command:
_________________________________________________________________________________
While this example is centered on the Dateline, one can easily choose another configuration with the -R option. A map centered on Greenwich would specify the region with -R-180/180/-70/70.
The transverse Mercator was invented by Lambert in 1772. In this projection the cylinder touches a meridian along which there is no distortion. The distortion increases away from the central meridian and goes to infinity at 90° from center. The central meridian, each meridian 90° away from the center, and equator are straight lines; other parallels and meridians are complex curves. The projection is defined by specifying:
The optional latitude of origin defaults to Equator if not specified. Although defaulting to 1, you can change the map scale factor via the MAP_SCALE_FACTOR parameter. Our example shows a transverse Mercator map of south-east Europe and the Middle East with 35°E as the central meridian:
_________________________________________________________________________________
The transverse Mercator can also be used to generate a global map—the equivalent of the 360° Mercator map. Using the command
_________________________________________________________________________________
we made the map illustrated in Figure 6.17. Note that when a world map is given (indicated by -R0/360/s/n), the arguments are interpreted to mean oblique degrees, i.e., the 360° range is understood to mean the extent of the plot along the central meridian, while the “south” and “north” values represent how far from the central longitude we want the plot to extend. These values correspond to latitudes in the regular Mercator projection and must therefore be less than 90°.
A particular subset of the transverse Mercator is the Universal Transverse Mercator (UTM) which was adopted by the US Army for large-scale military maps. Here, the globe is divided into 60 zones between 84°S and 84°N, most of which are 6° wide. Each of these UTM zones have their unique central meridian. Furthermore, each zone is divided into latitude bands but these are not needed to specify the projection for most cases. See Figure 6.18 for all zone designations.
GMT implements both the transverse Mercator and the UTM projection. When selecting UTM you must specify:
In order to minimize the distortion in any given zone, a scale factor of 0.9996 has been factored into the formulae. (although a standard, you can change this with MAP_SCALE_FACTOR). This makes the UTM projection a secant projection and not a tangent projection like the transverse Mercator above. The scale only varies by 1 part in 1,000 from true scale at equator. The ellipsoidal projection expressions are accurate for map areas that extend less than 10° away from the central meridian. For larger regions we use the conformal latitude in the general spherical formulae instead.
Oblique configurations of the cylinder give rise to the oblique Mercator projection. It is particularly useful when mapping regions of large lateral extent in an oblique direction. Both parallels and meridians are complex curves. The projection was developed in the early 1900s by several workers. Several parameters must be provided to define the projection. GMT offers three different definitions:
Our example was produced by the command
_________________________________________________________________________________
It uses definition 3 for an oblique view of some Caribbean islands. Note that we define our region using the rectangular system described earlier. If we do not append an “r” to the -R string then the information provided with the -R option is assumed to be oblique degrees about the projection center rather than the usual geographic coordinates. This interpretation is chosen since in general the parallels and meridians are not very suitable as map boundaries.
This cylindrical projection was developed in 1745 by César-François Cassini de Thury for the survey of France. It is occasionally called Cassini-Soldner since the latter provided the more accurate mathematical analysis that led to the development of the ellipsoidal formulae. The projection is neither conformal nor equal-area, and behaves as a compromise between the two end-members. The distortion is zero along the central meridian. It is best suited for mapping regions of north-south extent. The central meridian, each meridian 90° away, and equator are straight lines; all other meridians and parallels are complex curves. The requirements to define this projection are:
A detailed map of the island of Sardinia centered on the 8°45’E meridian using the Cassini projection can be obtained by running the command:
_________________________________________________________________________________
As with the previous projections, the user can choose between a rectangular boundary (used here) or a geographical (WESN) boundary.
This simple cylindrical projection is really a linear scaling of longitudes and latitudes. The most common form is the Plate Carrée projection, where the scaling of longitudes and latitudes is the same. All meridians and parallels are straight lines. The projection can be defined by:
The first two of these are optional and have defaults. When the standard parallel is defined, the central meridian must be supplied as well.
A world map centered on the dateline using this projection can be obtained by running the command:
_________________________________________________________________________________
Different relative scalings of longitudes and latitudes can be obtained by selecting a standard parallel different from the equator. Some selections for standard parallels have practical properties as shown in Table 6.1.
Projection | Standard parallel
|
Grafarend and Niermann, minimum linear distortion | 61.7° |
Ronald Miller Equirectangular | 50.5° |
Ronald Miller, minimum continental distortion | 43.5° |
Grafarend and Niermann | 42° |
Ronald Miller, minimum overall distortion | 37.5° |
Plate Carrée, Simple Cylindrical, Plain/Plane | 0° |
This cylindrical projection is actually several projections, depending on what latitude is selected as the standard parallel. However, they are all equal area and hence non-conformal. All meridians and parallels are straight lines. The requirements to define this projection are:
While you may choose any value for the standard parallel and obtain your own personal projection, there are seven choices of standard parallels that result in known (or named) projections. These are listed in Table 6.2.
Projection | Standard parallel
|
Balthasart | 50° |
Gall-Peters | 45° |
Hobo-Dyer | 37°30’ (= 37.5°) |
Trystan Edwards | 37°24’ (= 37.4°) |
Caster | 37°04’ (= 37.0666°) |
Behrman | 30° |
Lambert | 0° |
For instance, a world map centered on the 35°E meridian using the Behrman projection (Figure 6.22) can be obtained by running the command:
_________________________________________________________________________________
As one can see there is considerable distortion at high latitudes since the poles map into lines.
This cylindrical projection, presented by Osborn Maitland Miller of the American Geographic Society in 1942, is neither equal nor conformal. All meridians and parallels are straight lines. The projection was designed to be a compromise between Mercator and other cylindrical projections. Specifically, Miller spaced the parallels by using Mercator’s formula with 0.8 times the actual latitude, thus avoiding the singular poles; the result was then divided by 0.8. There is only a spherical form for this projection. Specify the projection by:
For instance, a world map centered on the 90°E meridian at a map scale of 1:400,000,000 (Figure 6.23) can be obtained as follows:
_________________________________________________________________________________
The cylindrical stereographic projections are certainly not as notable as other cylindrical projections, but are still used because of their relative simplicity and their ability to overcome some of the downsides of other cylindrical projections, like extreme distortions of the higher latitudes. The stereographic projections are perspective projections, projecting the sphere onto a cylinder in the direction of the antipodal point on the equator. The cylinder crosses the sphere at two standard parallels, equidistant from the equator. The projections are defined by:
Some of the selections of the standard parallel are named for the cartographer or publication that popularized the projection (Table 6.3).
Projection | Standard parallel
|
Miller’s modified Gall | 66.159467° |
Kamenetskiy’s First | 55° |
Gall’s stereographic | 45° |
Bolshoi Sovietskii Atlas Mira or Kamenetskiy’s Second | 30° |
Braun’s cylindrical | 0° |
A map of the world, centered on the Greenwich meridian, using the Gall’s stereographic projection (standard parallel is 45°, Figure 6.24), is obtained as follows:
_________________________________________________________________________________
GMT supports 8 common projections for global presentation of data or models. These are the Hammer, Mollweide, Winkel Tripel, Robinson, Eckert IV and VI, Sinusoidal, and Van der Grinten projections. Due to the small scale used for global maps these projections all use the spherical approximation rather than more elaborate elliptical formulae.
In all cases, the specification of the central meridian can be skipped. The default is the middle of the longitude range of the plot, specified by the -(R) option.
The equal-area Hammer projection, first presented by the German mathematician Ernst von Hammer in 1892, is also known as Hammer-Aitoff (the Aitoff projection looks similar, but is not equal-area). The border is an ellipse, equator and central meridian are straight lines, while other parallels and meridians are complex curves. The projection is defined by selecting:
A view of the Pacific ocean using the Dateline as central meridian is accomplished thus
_________________________________________________________________________________
This pseudo-cylindrical, equal-area projection was developed by the German mathematician and astronomer Karl Brandan Mollweide in 1805. Parallels are unequally spaced straight lines with the meridians being equally spaced elliptical arcs. The scale is only true along latitudes 40°44’ north and south. The projection is used mainly for global maps showing data distributions. It is occasionally referenced under the name homalographic projection. Like the Hammer projection, outlined above, we need to specify only two parameters to completely define the mapping of longitudes and latitudes into rectangular x/y coordinates:
An example centered on Greenwich can be generated thus:
_________________________________________________________________________________
In 1921, the German mathematician Oswald Winkel a projection that was to strike a compromise between the properties of three elements (area, angle and distance). The German word “tripel” refers to this junction of where each of these elements are least distorted when plotting global maps. The projection was popularized when Bartholomew and Son started to use it in its world-renowned “The Times Atlas of the World” in the mid 20th century. In 1998, the National Geographic Society made the Winkel Tripel as its map projection of choice for global maps.
Naturally, this projection is neither conformal, nor equal-area. Central meridian and equator are straight lines; other parallels and meridians are curved. The projection is obtained by averaging the coordinates of the Equidistant Cylindrical and Aitoff (not Hammer-Aitoff) projections. The poles map into straight lines 0.4 times the length of equator. To use it you must enter
Centered on Greenwich, the example in Figure 6.27 was created by this command:
_________________________________________________________________________________
The Robinson projection, presented by the American geographer and cartographer Arthur H. Robinson in 1963, is a modified cylindrical projection that is neither conformal nor equal-area. Central meridian and all parallels are straight lines; other meridians are curved. It uses lookup tables rather than analytic expressions to make the world map “look” right23. The scale is true along latitudes ±38°. The projection was originally developed for use by Rand McNally and is currently used by the National Geographic Society. To use it you must enter
Again centered on Greenwich, the example below was created by this command:
_________________________________________________________________________________
The Eckert IV and VI projections, presented by the German cartographer Max Eckert-Greiffendorff in 1906, are pseudocylindrical equal-area projections. Central meridian and all parallels are straight lines; other meridians are equally spaced elliptical arcs (IV) or sinusoids (VI). The scale is true along latitudes ±40°30’ (IV) and ±49°16’ (VI). Their main use is in thematic world maps. To select Eckert IV you must use -JKf (f for “four”) while Eckert VI is selected with -JKs (s for “six”). If no modifier is given it defaults to Eckert VI. In addition, you must enter
Centered on the Dateline, the Eckert IV example below was created by this command:
_________________________________________________________________________________
The same script, with s instead of f, yields the Eckert VI map:
The sinusoidal projection is one of the oldest known projections, is equal-area, and has been used since the mid-16th century. It has also been called the “Equal-area Mercator” projection. The central meridian is a straight line; all other meridians are sinusoidal curves. Parallels are all equally spaced straight lines, with scale being true along all parallels (and central meridian). To use it, you need to select:
A simple world map using the sinusoidal projection is therefore obtained by
_________________________________________________________________________________
To reduce distortion of shape the interrupted sinusoidal projection was introduced in 1927. Here, three symmetrical segments are used to cover the entire world. Traditionally, the interruptions are at 160°W, 20°W, and 60°E. To make the interrupted map we must call pscoast for each segment and superpose the results. To produce an interrupted world map (with the traditional boundaries just mentioned) that is 5.04 inches wide we use the scale 5.04/360° = 0.014 and offset the subsequent plots horizontally by their widths (140°0.014 and 80°0.014):
_________________________________________________________________________________
The usefulness of the interrupted sinusoidal projection is basically limited to display of global, discontinuous data distributions like hydrocarbon and mineral resources, etc.
The Van der Grinten projection, presented by Alphons J. van der Grinten in 1904, is neither equal-area nor conformal. Central meridian and Equator are straight lines; other meridians are arcs of circles. The scale is true along the Equator only. Its main use is to show the entire world enclosed in a circle. To use it you must enter
Centered on the Dateline, the example below was created by this command:
_________________________________________________________________________________
In this section we will be giving several examples of typical usage of GMT programs. In general, we will start with a raw data set, manipulate the numbers in various ways, then display the results in diagram or map view. The resulting plots will have in common that they are all made up of simpler plots that have been overlaid to create a complex illustration. We will mostly follow the following format:
A detailed discussion of each command is not given; we refer you to the manual pages for command line syntax, etc. We encourage you to run these scripts for yourself. See Appendix D if you would like an electronic version of all the shell-scripts (both sh and csh scripts are available, as or DOS batch files; only the sh-scripts are discussed here) and support data used below. Note that all examples explicitly specifies the measurement units, so although we use inches you should be able to run these scripts and get the same plots even if you have cm as the default measure unit. The examples are all written to be “quiet”, that is no information is echoed to the screen. Thus, these scripts are well suited for background execution.
Note that we also end each script by cleaning up after ourselves. Because awk is broken as designed on some systems, and nawk is not available on others we refer to $AWK in the scripts below; the do_examples.sh scripts will set this when running all examples.
Finally, be aware that for practical purposes the output PostScript file name is stored as the variable ps.
We want to create two contour maps of the low order geoid using the Hammer equal area projection. Our gridded data file is called osu91a1f_16.nc and contains a global 1° by 1° gridded geoid (we will see how to make gridded files later). We would like to show one map centered on Greenwich and one centered on the dateline. Positive contours should be drawn with a solid pen and negative contours with a dashed pen. Annotations should occur for every 50 m contour level, and both contour maps should show the continents in light gray in the background. Finally, we want a rectangular frame surrounding the two maps. This is how it is done:
_________________________________________________________________________________
The first command draws a box surrounding the maps. This is followed by two sequences of pscoast, grdcontour, grdcontour. They differ in that the first is centered on Greenwich; the second on the dateline. We use the limit option (-L) in grdcontour to select negative contours only and plot those with a dashed pen, then positive contours only and draw with a solid pen [Default]. The -T option causes tickmarks pointing in the downhill direction to be drawn on the innermost, closed contours. For the upper panel we also added - and + to the local lows and highs. You can find this illustration as Figure 7.1.
As our second example we will demonstrate how to make color images from gridded data sets (again, we will defer the actual making of grid files to later examples). We will use the supplemental program grdraster to extract 2-D grid files of bathymetry and Geosat geoid heights and put the two images on the same page. The region of interest is the Hawaiian islands, and due to the oblique trend of the island chain we prefer to rotate our geographical data sets using an oblique Mercator projection defined by the hotspot pole at (68°W, 69°N). We choose the point (190°, 25.5°) to be the center of our projection (e.g., the local origin), and we want to image a rectangular region defined by the longitudes and latitudes of the lower left and upper right corner of region. In our case we choose (160°, 20°) and (220°, 30°) as the corners. We use grdimage to make the illustration:
_________________________________________________________________________________
The first step extracts the 2-D data sets from the local data base using grdraster, which is a supplemental utility program (see Appendix A) that may be adapted to reflect the nature of your data base format. It automatically figures out the required extent of the region given the two corners points and the projection. The extreme meridians and parallels enclosing the oblique region is -R159:50/220:10/3:10/47:35. This is the area extracted by grdraster. For your convenience we have commented out those lines and provided the two extracted files so you do not need grdraster to try this example. By using the embedded grid file format mechanism we saved the topography using kilometers as the data unit. We now have two grid files with bathymetry and geoid heights, respectively. We use makecpt to generate a linear color palette file geoid.cpt for the geoid and use grd2cpt to get a histogram-equalized cpt file topo.cpt for the topography data. To emphasize the structures in the data we calculate the slopes in the north-south direction using grdgradient; these will be used to modulate the color image. Next we run grdimage to create a color-code image of the Geosat geoid heights, and draw a color legend to the right of the image with psscale. Similarly, we run grdimage but specify -Y4.5i to plot above the previous image. Adding scale and label the two plots a) and b) completes the illustration (Figure 7.2).
In this example we will show how to use the GMT programs fitcircle, project, sample1d, spectrum1d, psxy, and pstext. Suppose you have (lon, lat, gravity) along a satellite track in a file called sat.xyg, and (lon, lat, gravity) along a ship track in a file called ship.xyg. You want to make a cross-spectral analysis of these data. First, you will have to get the two data sets into equidistantly sampled time-series form. To do this, it will be convenient to project these along the great circle that best fits the sat track. We must use fitcircle to find this great circle and choose the L estimates of best pole. We project the data using project to find out what their ranges are in the projected coordinate. The minmax utility will report the minimum and maximum values for multi-column ASCII tables. Use this information to select the range of the projected distance coordinate they have in common. The script prompts you for that information after reporting the values. We decide to make a file of equidistant sampling points spaced 1 km apart from -1167 to +1169, and use the UNIX utility $AWK to accomplish this step. We can then resample the projected data, and carry out the cross-spectral calculations, assuming that the ship is the input and the satellite is the output data. There are several intermediate steps that produce helpful plots showing the effect of the various processing steps (example_03[a–f].ps), while the final plot example_03.ps shows the ship and sat power in one diagram and the coherency on another diagram, both on the same page. Note the extended use of pstext and psxy to put labels and legends directly on the plots. For that purpose we often use -Jx1i and specify positions in inches directly. Thus, the complete automated script reads:
_________________________________________________________________________________
The final illustration (Figure 7.3) shows that the ship gravity anomalies have more power than altimetry derived gravity for short wavelengths and that the coherency between the two signals improves dramatically for wavelengths 20 km.
This example will illustrate how to make a fairly complicated composite figure. We need a subset of the ETOPO5 bathymetry24 and Geosat geoid data sets which we will extract from the local data bases using grdraster. We would like to show a 2-layer perspective plot where layer one shows a contour map of the marine geoid with the location of the Hawaiian islands superposed, and a second layer showing the 3-D mesh plot of the topography. We also add an arrow pointing north and some text. This is how to do it:
_________________________________________________________________________________
The purpose of the color palette file zero.cpt is to have the positive topography mesh painted light gray (the remainder is white). The left side of Figure 7.4 shows the complete illustration.
A color version of this figure was used in our first article in EOS Trans. AGU (Oct. 8th, 1991). It was created along similar lines, but instead of a mesh plot we chose a color-coded surface with artificial illumination from a light-source due north. We choose to use the -Qi option in grdview to achieve a high degree of smoothness. Here, we select 100 dpi since that will be the resolution of our final raster (The EOS raster was 300 dpi). We used grdgradient to provide the intensity files. The following script creates the color PostScript file. Note that the size of the resulting output file is directly dependent on the square of the dpi chosen for the scanline conversion. A higher value for dpi in -Qi would have resulted in a much larger output file. The cpt files were taken from Section 7.2.
_________________________________________________________________________________
Instead of a mesh plot we may choose to show 3-D surfaces using artificial illumination. For this example we will use grdmath to make a grid file that contains the surface given by the function , where . The illumination is obtained by passing two grid files to grdview: One with the z-values (the surface) and another with intensity values (which should be in the ±1 range). We use grdgradient to compute the horizontal gradients in the direction of the artificial light source. The gray.cpt file only has one line that states that all z values should have the gray level 128. Thus, variations in shade are entirely due to variations in gradients, or illuminations. We choose to illuminate from the SW and view the surface from SE:
_________________________________________________________________________________
The variations in intensity could be made more dramatic by using grdmath to scale the intensity file before running grdview. For very rough data sets one may improve the smoothness of the intensities by passing the output of grdgradient to grdhisteq. The shell-script above will result in a plot like the one in Figure 7.5.
GMT provides two tools to render histograms: pshistogram and psrose. The former takes care of regular histograms whereas the latter deals with polar histograms (rose diagrams, sector diagrams, and wind rose diagrams). We will show an example that involves both programs. The file fractures.yx contains a compilation of fracture lengths and directions as digitized from geological maps. The file v3206.t contains all the bathymetry measurements from Vema cruise 3206. Our complete figure (Figure 7.6) was made running this script:
_________________________________________________________________________________
Many scientific papers start out by showing a location map of the region of interest. This map will typically also contain certain features and labels. This example will present a location map for the equatorial Atlantic ocean, where fracture zones and mid-ocean ridge segments have been plotted. We also would like to plot earthquake locations and available isochrons. We have obtained one file, quakes.xym, which contains the position and magnitude of available earthquakes in the region. We choose to use magnitude/100 for the symbol-size in inches. The digital fracture zone traces (fz.xy) and isochrons (0 isochron as ridge.xy, the rest as isochrons.xy) were digitized from available maps25. We create the final location map (Figure 7.7) with the following script:
_________________________________________________________________________________
The same figure could equally well be made in color, which could be rasterized and made into a slide for a meeting presentation. The script is similar to the one outlined above, except we would choose a color for land and oceans, and select colored symbols and pens rather than black and white.
The program psxyz allows us to plot three-dimensional symbols, including columnar plots. As a simple demonstration, we will convert a gridded netCDF of bathymetry into an ASCII table and use the height information to draw a 2-D histogram in a 3-D perspective view. Our gridded bathymetry file is called guinea_bay.nc and covers the region from 0 to 5 °E and 0 to 5 °N. Depth ranges from -5000 meter to sea-level. We produce the Figure 7.8 by running this script:
_________________________________________________________________________________
A common application in many scientific disciplines involves plotting one or several time-series as as “wiggles” along tracks. Marine geophysicists often display magnetic anomalies in this manner, and seismologists use the technique when plotting individual seismic traces. In our example we will show how a set of Geosat sea surface slope profiles from the south Pacific can be plotted as “wiggles” using the pswiggle program. We will embellish the plot with track numbers, the location of the Pacific-Antarctic Ridge, recognized fracture zones in the area, and a “wiggle” scale. The Geosat tracks are stored in the files *.xys, the ridge in ridge.xy, and all the fracture zones are stored in the multiple segment file fz.xy. We extract the profile id (which is the first part of the file name for each profile) and the last point in each of the track files to construct an input file for pstext that will label each profile with the track number. We know the profiles trend approximately N40°E so we want the labels to have that same orientation (i.e., the angle with the baseline must be 50°). We do this by extracting the last record from each track, paste this file with the tracks.lis file, and use $AWK to create the format needed for pstext. Note we offset the positions by -0.05 inch with -D in order to have a small gap between the profile and the label:
_________________________________________________________________________________
The output shows the sea-surface slopes along 42 descending Geosat tracks in the Eltanin and Udintsev fracture zone region in a Mercator projection (Figure 7.9).
Our next and perhaps most business-like example presents a three-dimensional bar graph plot showing the geographic distribution of the membership in the American Geophysical Union (AGU). The input data was taken from the January 2008 AGU member directory and added up to give total members per continent. We decide to plot a 3-D column centered on each continent with a height that is proportional to the logarithm of the membership. A log-scale is used since the memberships vary by almost 3 orders of magnitude. We choose a plain linear projection for the basemap and add the columns and text on top. Our script that produces Figure 7.10 reads:
_________________________________________________________________________________
In this example we generate a series of 6 color images, arranged so that they can be cut out and assembled into a 3-D color cube. The six faces of the cube represent the outside of the R-G-B color space. On each face one of the color components is fixed at either 0 or 255 and the other two components vary smoothly across the face from 0 to 255. The cube is configured as a right-handed coordinate system with x-y-z mapping R-G-B. Hence, the 8 corners of the cube represent the primaries red, green, and blue, plus the secondaries cyan, magenta and yellow, plus black and white.
The 6 color faces are generated by feeding grdimage three grids, one for each color component (R, G, and B). In some cases the X or Y axes of a face are reversed by specifying a negative width or height in order to change the variation of the color value in that direction from ascending to descending, or vice versa.
A number of rays emanating from the white and black corners indicate the Hue value (ranging from 0 to 360°). The dashed and dotted lines near the white corner reflect saturation levels, running from 0 to 1 (in black font). On these 3 faces the brightness is a constant value of 1. On the other 3 faces of the cube, around the black corner, the white decimal numbers indicate brightnesses between 0 and 1, with saturation fixed at 1.
Here is the shell script to generate the RGB cube in Figure 7.11:
_________________________________________________________________________________
Our next example (Figure 7.12) operates on a data set of topographic readings non-uniformly distributed in the plane (Table 5.11 in Davis: Statistics and Data Analysis in Geology, J. Wiley). We use triangulate to perform the optimal Delaunay triangulation, then use the output to draw the resulting network. We label the node numbers as well as the node values, and call pscontour to make a contour map and image directly from the raw data. Thus, in this example we do not actually make grid files but still are able to contour and image the data. We use a color palette table topo.cpt (created via minmax and makecpt). The script becomes:
_________________________________________________________________________________
In many areas, such as fluid dynamics and elasticity, it is desirable to plot vector fields of various kinds. GMT provides a way to illustrate 2-component vector fields using the grdvector utility. The two components of the field (Cartesian or polar components) are stored in separate grid files. In this example we use grdmath to generate a surface and to calculate by returning the x- and y-derivatives separately. We superpose the gradient vector field and the surface z and also plot the components of the gradient in separate windows. A pstext call to place a header finishes the plot (Figure 7.13:
_________________________________________________________________________________
This example shows how one goes from randomly spaced data points to an evenly sampled surface. First we plot the distribution and values of our raw data set (same as in Section 7.12). We choose an equidistant grid and run blockmean which preprocesses the data to avoid aliasing. The dashed lines indicate the logical blocks used by blockmean; all points inside a given bin will be averaged. The logical blocks are drawn from a temporary file we make on the fly within the shell script. The processed data is then gridded with the surface program and contoured every 25 units. A most important point here is that blockmean, blockmedian, or blockmode should always be run prior to running surface, and both of these steps must use the same grid interval. We use grdtrend to fit a bicubic trend surface to the gridded data, contour it as well, and sample both grid files along a diagonal transect using grdtrack. The bottom panel compares the gridded (solid line) and bicubic trend (dashed line) along the transect using psxy (Figure 7.14):
_________________________________________________________________________________
This example (Figure 7.15) demonstrates some off the different ways one can use to grid data in GMT, and how to deal with unconstrained areas. We first convert a large ASCII file to binary with gmtconvert since the binary file will read and process much faster. Our lower left plot illustrates the results of gridding using a nearest neighbor technique (nearneighbor) which is a local method: No output is given where there are no data. Next (lower right), we use a minimum curvature technique (surface) which is a global method. Hence, the contours cover the entire map although the data are only available for portions of the area (indicated by the gray areas plotted using psmask). The top left scenario illustrates how we can create a clip path (using psmask) based on the data coverage to eliminate contours outside the constrained area. Finally (top right) we simply employ pscoast to overlay gray land masses to cover up the unwanted contours, and end by plotting a star at the deepest point on the map with psxy. This point was extracted from the grid files using grdinfo.
_________________________________________________________________________________
pscontour (for contouring) and triangulate (for gridding) use the simplest method of interpolating data: a Delaunay triangulation (see Section 7.12) which forms as a union of planar triangular facets. One advantage of this method is that it will not extrapolate beyond the convex hull of the input (x, y) data. Another is that it will not estimate a z value above or below the local bounds on any triangle. A disadvantage is that the surface is not differentiable, but has sharp kinks at triangle edges and thus also along contours. This may not look physically reasonable, but it can be filtered later (last panel below). surface can be used to generate a higher-order (smooth and differentiable) interpolation of onto a grid, after which the grid may be illustrated (grdcontour, grdimage, grdview). surface will interpolate to all (x, y) points in a rectangular region, and thus will extrapolate beyond the convex hull of the data. However, this can be masked out in various ways (see Section 7.15).
A more serious objection is that surface may estimate z values outside the local range of the data (note area near x = 0.8, y = 5.3). This commonly happens when the default tension value of zero is used to create a “minimum curvature” (most smooth) interpolant. surface can be used with non-zero tension to partially overcome this problem. The limiting value should approximate the triangulation, while a value between 0 and 1 may yield a good compromise between the above two cases. A value of 0.5 is shown here (Figure 7.16). A side effect of the tension is that it tends to make the contours turn near the edges of the domain so that they approach the edge from a perpendicular direction. A solution is to use surface in a larger area and then use grdcut to cut out the desired smaller area. Another way to achieve a compromise is to interpolate the data to a grid and then filter the grid using grdfft or grdfilter. The latter can handle grids containing “NaN” values and it can do median and mode filters as well as convolutions. Shown here is triangulate followed by grdfilter. Note that the filter has done some extrapolation beyond the convex hull of the original x, y values. The “best” smooth approximation of depends on the errors in the data and the physical laws obeyed by z. GMT cannot always do the “best” thing but it offers great flexibility through its combinations of tools. We illustrate all four solutions using a cpt file that contains color fills, predefined patterns for interval (900,925) and NaN, an image pattern for interval (875,900), and a “skip slice” request for interval (700,725).
_________________________________________________________________________________
This example demonstrates how pscoast can be used to set up clip paths based on coastlines. This approach is well suited when different gridded data sets are to be merged on a plot using different color palette files. Merging the files themselves may not be doable since they may represent different data sets, as we show in this example. Here, we lay down a color map of the geoid field near India with grdimage, use pscoast to set up land clip paths, and then overlay topography from the ETOPO5 data set with another call to grdimage. We finally undo the clippath with a second call to pscoast with the option -Q (Figure 7.17):
_________________________________________________________________________________
We also plot a color legend on top of the land. So here we basically have three layers of “paint” stacked on top of each other: the underlaying geoid map, the land mask, and finally the color legend. This legend makes clear how grd2cpt distributed the colors over the range: they are not of equal length put are associated with equal amounts of area in the plot. Since the high amounts (in red) are not very prevalent, that color spans a long range.
For this image it is appropriate to use the -I option in psscale so the legend gets shaded, similar to the geoid grid. See Appendix M to learn more about color palettes and ways to draw color legends.
To demonstrate potential usage of the new programs grdvolume and gmtselect we extract a subset of the Sandwell & Smith altimetric gravity field26 for the northern Pacific and decide to isolate all seamounts that (1) exceed 50 mGal in amplitude and (2) are within 200 km of the Pratt seamount. We do this by dumping the 50 mGal contours to disk, then making a simple $AWK script center.awk that returns the mean location of the points making up each closed polygon, and then pass these locations to gmtselect which retains only the points within 200 km of Pratt. We then mask out all the data outside this radius and use grdvolume to determine the combined area and volumes of the chosen seamounts. Our illustration is presented in Figure 7.18.
_________________________________________________________________________________
GMT 3.1 introduced color patterns and this examples give a few cases of how to use this new feature. We make a phony poster that advertises an international conference on GMT in Honolulu. We use grdmath, makecpt, and grdimage to draw pleasing color backgrounds on maps, and overlay pscoast clip paths to have the patterns change at the coastlines. The middle panel demonstrates a simple pscoast call where the built-in pattern # 86 is drawn at 100 dpi but with the black and white pixels replaced with color combinations. At the same time the ocean is filled with a repeating image of a circuit board (provides in Sun raster format). The text GMT in the center is an off-line PostScript file that was overlaid using psimage. The final panel repeats the top panel except that the land and sea images have changed places (Figure 7.19).
_________________________________________________________________________________
One is often required to make special maps that shows the distribution of certain features but one would prefer to use a custom symbol instead of the built-in circles, squares, triangles, etc. in the GMT plotting programs psxy and psxyz. Here we demonstrate one approach that allows for a fair bit of flexibility in designing ones own symbols. The following recipe is used when designing a new symbol.
M [ -Gfill ] [ -Wpen ] | Start new element at , |
D | Draw straight line from current point to , around (, ) |
A | Draw arc segment of radius from angle to |
We also add a few stand-alone circles (for other symbols, see psxy man page):
c [ -Gfill ] [ -Wpen ] | Draw single circle of radius around , |
The optional -G and -W can be used to hardwire the color fill and pen for segments (use – to disallow fill or line for any specific feature). By default the segments are painted based on the values of the command line settings.
Manually applying these rules to our volcano symbol results in a definition file volcano.def:
______________________________________________________________________________
Without much further discussion we also make a definition file bullseye.def for a multi-colored bulls eye symbol. Note that the symbol can be created beyond the -0.5 to +0.5 range, as shown by the red lines. There is no limit in GMT to the size of the symbols. The center, however, will always be at (0,0). This is the point to which the coordinates in psxy refers.
______________________________________________________________________________
The values refer to positions and dimensions illustrated in the Figure above.
We are now ready to give it a try. Based on the hotspot locations in the file hotspots.d (with a 3rd column giving the desired symbol sizes in inches) we lay down a world map and overlay red volcano symbols using our custom-built volcano symbol and psxy. We do something similar with the bulls eye symbols. Without the -G option, however, they get the colors defined in bullseye.def.
Here is our final map script that produces Figure 7.20:
_________________________________________________________________________________
Given these guidelines you can easily make your own symbols. Symbols with more than one color can be obtained by making several symbol components. E.g., to have yellow smoke coming out of red volcanoes we would make two symbols: one with just the cone and caldera and the other with the bubbles. These would be plotted consecutively using the desired colors. Alternatively, like in bullseye.def, we may specify colors directly for the various segments. Note that the custom symbols (Appendix N), unlike the built-in symbols in GMT, can be used with the built-in patterns (Appendix E). Other approaches are also possible, of course.
As discussed in Section 4.4.3, the annotation of time-series is generally more complicated due to the extra degrees of freedom afforded by the dual annotation system. In this example we will display the trend of the stock price of RedHat (RHAT) from their initial public offering until late 2006. The data file is a comma-separated table and the records look like this:
Hence, we have a single header record and various prices in USD for each day of business. We will plot the trend of the opening price as a red line superimposed on a yellow envelope representing the low-to-high fluctuation during each day. We also indicate when and at what cost Paul Wessel bought a few shares, and zoom in on the developments since 2004; in the inset we label the time-axis in Finnish in honor of Linus Thorvalds. Because the time coordinates are Y2K-challenged and the order is backwards (big units of years come after smaller units like days) we must change the default input/output formats used by GMT. Finally, we want to prefix prices with the $ symbol to indicate the currency. Here is how it all comes out:
_________________________________________________________________________________
which produces the plot in Figure 7.21, suggesting Wessel has missed a few trains if he had hoped to cash in on the Internet bubble...
The next example uses the command-line tool wget to obtain a data file from a specified URL27. In the example script this line is commented out so the example will run even if you do not have wget (we use the supplied neic_quakes.d which normally would be created by wget); remove the comment to get the actual current seismicity plot using the live data. The main purpose of this script is not to show how to plot a map background and a few circles, but rather demonstrate how a map legend may be composed using the new tool pslegend. Some scripting is used to pull out information from the data file that is later used in the legend. The legend will normally have the email address of the script owner; here that command is commented out and the user is hardwired to “GMT guru”. The USGS logo, taken from their web page and converted to a Sun raster file, is used to spice up the legend.
_________________________________________________________________________________
The script produces the plot in Figure 7.22, giving the URL where these and similar data can be obtained.
While motorists recently have started to question the old saying “all roads lead to Rome”, aircraft pilots have known from the start that only one great-circle path connects the points of departure and arrival28. This provides the inspiration for our next example which uses grdmath to calculate distances from Rome to anywhere on Earth and grdcontour to contour these distances. We pick five cities that we connect to Rome with great circle arcs, and label these cities with their names and distances (in km) from Rome, all laid down on top of a beautiful world map. Note that we specify that contour labels only be placed along the straight map-line connecting Rome to its antipode, and request curved labels that follows the shape of the contours.
_________________________________________________________________________________
The script produces the plot in Figure 7.23; note how interesting the path to Seattle appears in this particular projection (Hammer). We also note that Rome’s antipode lies somewhere near the Chatham plateau (antipodes will be revisited in Section 7.25).
Although we are not seismologists, we have yet another example involving seismicity. We use seismicity data for the Australia/New Zealand region to demonstrate how we can extract subsets of data using geospatial criteria. In particular, we wish to plot the epicenters given in the file oz_quakes.d as red or green circles. Green circles should only be used for epicenters that satisfy the following three criteria:
All remaining earthquakes should be plotted in red. Rather that doing the selection process twice we simply plot all quakes as red circles and then replot those that pass our criteria. Most of the work here is done by gmtselect; the rest is carried out by the usual pscoast and psxy workhorses. Note for our purposes the Dateline is just a line along the 180° meridian.
_________________________________________________________________________________
The script produces the plot in Figure 7.24. Note that the horizontal distance from the dateline seems to increase as we go south; however that is just the projected distance (Mercator distortion) and not the actual distance which remains constant at 1000 km.
As promised in Section 7.23, we will study antipodes. The antipode of a point at is the point at . We seek an answer to the question that has plagued so many for so long: Given the distribution of land and ocean, how often is the antipode of a point on land also on land? And what about marine antipodes? We use grdlandmask and grdmath to map these distributions and calculate the area of the Earth (in percent) that goes with each of the three possibilities. To make sense of our grdmath equations below, note that we first calculate a grid that is +1 when a point and its antipode is on land, -1 if both are in the ocean, and 0 elsewhere. We then seek to calculate the area distribution of dry antipodes by only pulling out the nodes that equal +1. As each point represent an area approximated by where the term’s actual dimension depends on , we need to allow for that shrinkage, normalize our sum to that of the whole area of the Earth, and finally convert that ratio to percent. Since the , terms appear twice in these expressions they cancel out, leaving the somewhat intractable expressions below where the sum of for all is known to equal :
_________________________________________________________________________________
In the end we obtain a funny-looking map depicting the antipodal distribution as well as displaying in legend form the requested percentages (Figure 7.25). Note that the script is set to evaluate a global 30 minute grid for expediency (), hence several smaller land masses that do have terrestrial antipodes do not show up. If you want a more accurate map you can set the parameter to a smaller increment (try 5 and wait a few minutes).
The call to grdimage includes the –Sn to suspend interpolation and only return the value of the nearest neighbor. This option is particularly practical for plotting categorical data, like these, that should not be interpolated.
Next, we present a recent extension to the -JG projection option which allows the user to specify a particular altitude (this was always at infinity before), as well as several further parameters to limit the view from the chosen vantage point. In this example we show a view of the eastern continental US from a height of 160 km. Below we add a view with a specific tilt of 55° and azimuth 210°; here we have chosen a boresight twist of 45°. We view the land from New York towards Washington, D.C.
_________________________________________________________________________________
At this point the full projection has not been properly optimized and the map annotations will need additional work. Also, note that the projection is only implemented in pscoast and grdimage. We hope to refine this further and extend the availability of the full projection to all of the GMT mapping programs.
Next, we show how to plot a data grid that is distributed in projected form. The gravity and predicted bathymetry grids produced by David Sandwell and Walter H. F. Smith are not geographical grids but instead given on a spherical Mercator grid. The GMT supplement imgsrc has tools to extract subsets of these large grids. If you need to make a non-Mercator map then you must extract a geographic grid using img2grd and then plot it using your desired map projection. However, if you want to make a Mercator map then you can save time and preserve data quality by avoiding to re-project the data set twice since it is already in a Mercator projection. This example shows how this is accomplished. We use the -M option in img2grd29 to pull out the grid in Mercator units (i.e., do not invert the Mercator projection) and then simply plot the grid using a linear projection with a suitable scale (here 0.25 inches per degrees of longitude). To overlay basemaps and features that has geographic longitude/latitude coordinates we must remember two key issues:
_________________________________________________________________________________
This map of the Tasman Sea shows the marine gravity anomalies with land painted black. A color scale bar was then added to complete the illustration.
Next, we present a similar case: We wish to plot a data set given in UTM coordinates and want it to be properly registered with overlying geographic data, such as coastlines or data points. The mistake many GMT rookies make is to specify the UTM projection with their UTM data. However, that data have already been projected and is now in linear meters. The only sensible way to plot such data is with a linear projection, yielding a UTM map. In this step one can choose to annotate or tick the map in UTM meters as well. To plot geographic (lon/lat) data on the same map there are a few things you must consider:
_________________________________________________________________________________
Our script illustrates how we would plot a UTM grid of elevations near Kilauea volcano on the Big Island of Hawaii. Given we are in UTM zone 5Q, the script determines the geographic coordinates of the lower left and upper right corner of the UTM grid, then uses that region when overlaying the coastline and light blue ocean. We place a scale bar and label Kilauea crater to complete the figure.
Next, we demonstrate how gridding on a spherical surface can be accomplished using Green’s functions of surface splines, with or without tension. Global gridding does not work particularly well in Cartesian coordinates hence the chosen approach. We use greenspline to produce a crude topography grid for Mars based on radii estimates from the Mariner 9 and Viking Orbiter spacecrafts. This data comes from Smith and Zuber [Science, 1996] and is used here as a small (N = 370) data set we can use to demonstrate spherical surface gridding. Since greenspline must solve a N by N matrix system your system memory may impose limits on how large data sets you can handle; also note that the spherical surface spline in tension is particularly slow to compute.
_________________________________________________________________________________
Our script must first estimate the ellipsoidal shape of Mars from the parameters given by Smith and Zuber so that we can remove this reference surface from the gridded radii. We run the gridding twice: First with no tension using Parker’s [1990] method and then with tension using the Wessel and Becker [2008] method. The grids are then imaged with grdimage and grdcontour and a color scale is placed between them.
Finally, we end with a simple mathematical illustration of sine and cosine, highlighting the it graph mode for linear projections and the new curved vectors for angles.
_________________________________________________________________________________
The script simply draws a graph basemap, computes sine and cosine and plots them as lines, then indicates on a circle that these quantities are simply the projections of an unit vector on the x- and y-axis, at the given angle.
Unlike the previous chapter, in this chapter we will explore what is involved in creating animations (i.e., movies). Of course, an animation is nothing more than a series of individual images played back in an orderly fashion. Here, these images will have been created with GMT. To ensure a smooth transition from frame to frame we will be following some general guidelines when writing our scripts. Since there is no “movie” mode in GMT we must take care of all the book-keeping in our script. Thus, animations may require a bit of planning and may use more advanced scripting than the previous static examples. Note: This is a new chapter introduced with the 4.4.0 version and should be considered work in progress.
Most, if not all, animation scripts must deal with several specific phases of movie making:
We will discuss these phases in more detail before showing our first example.
Note that in our examples below the animation scripts only produce a PostScript plot of the first frame, then exits. This is to allow our install script to just produce the frame for inclusion in our documentation, and to avoid a lengthy movie production when installing GMTṪo actually build these animation examples you must run the scripts with an argument (any argument, it does not matter).
Our first animation is not very ambitious: We wish to plot the sine function from 0–360° and take snap shots every 20°. To get a smooth curve we must sample the function much more frequently; we settle on 10 times more frequently than the frame spacing. We place a bright red circle at the leading edge of the curve, and as we move forward in time (here, angles) we dim the older circles to a dark red color. We add a label that indicates the current angle value. Once the 18 frames are completed we convert them to a single animated GIF file and write a plain HTML wrapper with a simple legend. Opening the HTML page anim01.html in the browser will display the animation.
_________________________________________________________________________________
Make sure you understand the purpose of all the steps in our script. In this case we did some trial-and-error to determine the exact values to use for the map projection, the region, the spacing around the frame, etc. so that the final result gave a reasonable layout. Do this planning on a single PostScript plot before running a lengthy animation script.
Our next animation uses a gridded topography for parts of Colorado (US); the file is distributed with the tutorial examples. Here, we want to use grdimage to generate a shaded-relief image sequence in which we sweep the illumination azimuth around the entire horizon. The resulting animation illustrates how changing the illumination azimuth can bring out subtle features (or artifacts) in the gridded data. The red arrow points in the direction of the illumination.
_________________________________________________________________________________
As you can see, these sorts of animations are not terribly difficult to put together, especially since our vantage point is fixed. In the next example we will move the “camera” around and must therefore deal with how to frame perspective views.
Our third animation keeps a fixed gridded data set but moves the camera angle around the full 360°. We use grdview to generate a shaded-relief image sequence using the new enhanced -E option. No additional information is plotted on the image. As before we produce an animated GIF image and a simple HTML wrapper for it.
_________________________________________________________________________________
Our next animation simulates what an imaginary satellite might see as it passes in a great circle from New York to Miami at an altitude of 160 km. We use the general perspective view projection with grdimage and use project to create a great circle path between the two cities, sampled every 5 km. The main part of the script will make the DVD-quality frames from different view points, draw the path on the ground, and add frame numbers to each frame. As this animation generates 355 frames we can use 3rd party tools to turn the image sequence into a MPEG-4 movie31. Note: At the moment, grdview cannot use general perspective view projection to allow “fly-through” animations like Fledermaus; we expect to add this functionality in a future version.
_________________________________________________________________________________
Most public-domain (and even commercial) software comes with bugs, and the speed with which such bugs are detected and removed depends to a large degree on the willingness of the user community to report these to us in a useful manner. When your car breaks down, simply telling the mechanic that it doesn’t work will hardly speed up the repair or cut back costs! Therefore, we ask that if you detect a bug, first make sure that it in fact is a bug and not a user error. Then, send us email about the problem. Be sure to include all the information necessary for us to recreate the situation in which the bug occurred. This will include the full command line used and, if possible, the data file used by the program. Send the bug-reports to gmt-help@lists.hawaii.edu. We will try to fix bugs as soon as our schedules permit and inform users about the bug and availability of updated code (See Appendix D).
Two electronic mailing lists are available to which users may subscribe. gmt-group@lists.hawaii.edu and is primarily a way for us to notify the users when bugs have been fixed or when new updates have been installed in the ftp directory (See Appendix D). We also maintain another list (gmt-help@lists.hawaii.edu) which interested users may subscribe to. It basically provides a forum for GMT users to exchange ideas and ask questions about GMT usage, installation and portability, etc. Please use this utility rather than sending questions directly to us personally. We hope you appreciate that we simply do not have time to be everybody’s personal GMT tutor.
The electronic mailing lists are maintained automatically by a program. To subscribe to one or both of the lists, send a message to listserv@lists.hawaii.edu containing the command(s):
subscribe gmt-group your
full name, not email address
subscribe gmt-help your
full name, not email address
(Do not type the angular brackets ). You may also register electronically via the GMT home web page (gmt.soest.hawaii.edu). For information on what commands you may send, send a message containing the word help. You must interact with the listserver to be added to or removed from the mailing lists! We strongly recommend that you at least subscribe to gmt-group since this is how we can notify you of future updates and bug-fixes. Most new users will also benefit from having the other forum (gmt-help) as they struggle to realign their sense of logic with that of GMT. While anybody may post messages to gmt-help, access to gmt-group is restricted to minimize net traffic. Any message sent to gmt-group will be intercepted by the GMT manager who will determine if the message is important enough to cause thousands of mailtools to go BEEP. Communication with other GMT users should go via gmt-help. Finally, all GMT information is provided online at the main GMT home page in Hawaii, i.e., gmt.soest.hawaii.edu. Changes to GMT will also be posted on this page. The main GMT page has links to the official GMT ftp sites around the world.
These packages are for the most part written and supported by us, but there are some exceptions. They provide extensions to GMT that are needed for particular rather than general applications. The packages are provided with GMT and are installed by default unless they requires non-GMT libraries; see the main GMT configuration process. Questions or bug reports for this software should be addressed to the person(s) listed in the README file associated with the particular program. It is not guaranteed that these programs are fully ANSI-C, Y2K, or POSIX compliant, or that they necessarily will install smoothly on all platforms, but most do. Note that the data sets some of these programs work on are not distributed with these packages; they must be obtained separately. The contents of the supplemental archive may change without notice; at this writing it contains these directories:
This package contains grdraster which you can use to extract data from global gridded data sets such as those available from NGDC. We have used it to prepare some of the grids in the examples (Chapter 6). You can also customize it to read your own data sets. The package is maintained by the GMT developers.
This package contains gshhg which you can use to extract shoreline polygons from the Global Self-consistent Hierarchical High-resolution Geography (GSHHG) available separately from NGDC or the GSHHG home page (GSHHG contains both GSHHS and WDBII; GSHHS is the polygon data base from which the GMT coastlines derive while WDBII provides rivers and boundaries). It also contains gshhg_dp for cleverly decimating a shoreline, and gshhgtograss to convert shoreline segments to the GRASS database format; the latter program is maintained by Simon Cox32. The package is maintained by Paul Wessel.
This package consists of the program img2mercgrd to extract subsets of the global gravity and predicted topography solutions derived from satellite altimetry33. The package is maintained by Walter Smith34.
This package contains the programs pscoupe, psmeca, pspolar, and psvelo which are used by seismologists and geodesists for plotting focal mechanisms (including cross-sections and polarities), error ellipses, velocity arrows, rotational wedges, and more. The package is maintained by Kurt Feigl35 and Genevieve Patau36.
Here you will find the mex files grdinfo, grdread, and grdwrite, which can be used in Matlab or Octave to read and write grid files. The package originated with David Sandwell, UCSD, and was subsequently modified by Paul Wessel and Phil Sharfstein, UCSB. It is now maintained by Paul Wessel.
This package currently holds the programs mgd77convert, mgd77header, mgd77info, mgd77list, mgd77magref, mgd77manage, mgd77path, mgd77sniffer, and mgd77track which can be used to extract information or data values from or plot marine geophysical data files in the ASCII MGD77 or netCDF MGD77+ formats37). We expect this package eventually to replace the mgg package. The package is maintained by Paul Wessel.
This package holds the legacy programs binlegs, dat2gmt, gmt2dat, gmtinfo, gmtlegs, gmtlist, gmtpath, gmttrack, and mgd77togmt, which can be used to maintain, access, extract data from, and plot marine geophysical data files converted from the MGD77 format to the .gmt format38). The package is maintained by the GMT developers.
At the moment, this package contains the programs dimfilter, which is an extension of grdfilter in that it allows for spatial directional filtering, psmegaplot which you can use to make large posters using a simple laserwriter, makepattern which generates raster patterns from GMT 3.0 icon files, gmt2kml which converts GMT table data to Google Earth’s KML format, gmtdigitize which provides a GMT interface to a digitizing tablet via a serial port, gmtstitch which can be used to assemble pieces digitized lines into complete lines or polygons, gmtdp which performs line reduction using the Douglas-Peucker algorithm, kml2gmt which extracts GMT table data from Google Earth KML files, and nc2xy which can extract data from column-oriented netCDF files. The package is maintained by Paul Wessel. The increasingly popular utility ps2raster, which simplifies the rasterization of GMTPostScript to raster formats (see Appendix C), was moved to the general tools starting with GMT 4.2.0.
This package contains programs to plot SEGY seismic data files using the GMT mapping transformations and postscript library. pssegy generates a 2-D plot (x:location and y:time/depth) while pssegyz generates a 3-D plot (x and y: location coordinates, z: time/depth). Locations may be read from predefined or arbitrary portions of each trace header. Finally, segy2grd can convert SEGY data to a GMT grid file. The package is maintained by Tim Henstock39.
This package contains the main programs sphtriangulate, which you can use to generate data for Delaunay or Voronoi diagrams, sphdistance which calculates distances from lines to grid nodes using Voronoi decomposition of the data, and sphinterpolate which performs gridding under tension on a sphere. These programs passes the heavy work onto the two Fortran-77 packages SSRFPACK and STRIPACK by Robert Renka; here they have been translated to C with assistance from f2c. The package is maintained by Paul Wessel.
This package contains the plate tectonic programs backtracker, which you can use to move geologic markers forward or backward in time, grdrotater which rotates entire grids using a finite rotation, hotspotter which generates CVA grids based on seamount locations and a set of absolute plate motion stage poles (grdspotter does the same using a bathymetry grid instead of seamount locations), originator, which associates seamounts with the most likely hotspot origins, and rotconverter which does various operations involving finite rotations on a sphere. The package is maintained by Paul Wessel.
This package contains the tools x2sys_datalist, which allows you to extract data from almost any binary or ASCII data file, and x2sys_cross which determines crossover locations and errors generated by one or several geospatial tracks. Newly added are the tools x2sys_init, x2sys_binlist, x2sys_get, x2sys_list, x2sys_put, x2sys_report, x2sys_solve and x2sys_merge which extends the track-management system employed by the mgg supplement to generic track data of any format. This package represents a new generation of tools intended to replace the old “X_SYSTEM” crossover tools (below). The package is maintained by Paul Wessel.
This package contains the tools x_edit, x_init, x_list, x_over, x_remove, x_report, x_setup, x_solve_dc_drift, and x_update. Collectively, they make up the old “XSYSTEM” crossover tools. This package with remain in the GMT supplemental archive until x2sys is complete. The package is maintained by Paul Wessel.
The package contains an X11 editor (xgridedit) for visual editing of grid files. It was originally developed by Hugh Fisher, CRES, in March 1992 but is now maintained by Lloyd Parkes40.
These files have N records which have M fields each. Most programs can read multicolumn files, but require that the x [and y] variable(s) be stored in the 1st [and 2nd] column (There are, however, some exceptions to this rule, such as filter1d and sample1d). GMT can read both ASCII and binary table data.
The first data record may be preceded by 1 or more header records. When using such files, make sure to use the -H option and set the parameter N_HEADER_RECS in the .gmtdefaults4 file (System default is 1 header record if -H is set; you may also use -Hnrecs directly). Fields within a record must be separated by spaces, tabs, or commas. Each field can be an integer or floating-point number or a geographic coordinate string using the [+-]dd[:mm[:ss]][WSNEwsne] format. Thus, 12:30:44.5W, 17.5S, 1:00:05, and 200:45E are all valid input strings.
When dealing with time- or (x,y)-series it is usually convenient to have each profile in separate files. However, this may sometimes prove impractical due to large numbers of profiles. An example is files of digitized lineations where the number of individual features may range into the thousands. One file per feature would in this case be unreasonable and furthermore clog up the directory. GMT provides a mechanism for keeping more than one profile in a file. Such files are called multiple segment files and are identical to the ones just outlined except that they have subheaders interspersed with data records that signal the start of a segment. The subheaders may be of any format, but all must have the same character in the first column. When using such files, you must specify the -m option. The unique character is by default ’’, but you can override that by appending your chosen character to the M option. E.g., -mH will look for subheaders starting with H, whereas -m’*’ will check for asterisks (The quotes are necessary since * has special meaning to UNIX). Some programs such as psxy will examine the subheaders to see if they contain -W and -G options for specifying pen and fill attributes for individual segments, -Z to change color via a cpt-file, or -L for label specifications. These settings (and occasionally others) will override the corresponding command line options.
GMT programs also support native binary tables to speed up input-output for i/o-intensive tasks like gridding and preprocessing. Files may have no header (hence the -H option cannot be used) and all data must either be single or double precision (no mixing allowed). Multiple segment files are allowed (-m) and the segment headers are assumed to be records where all the fields equal NaN. Flags appended to -m are ignored. The format and number of fields are specified with the -b option. Thus, for input you may set -bi[s][n], where s designates single precision (default is d for double) and n is the number of fields. For output, use -bo[s] (the programs know how many columns to write, unless you use -m in which case we need to know the number of output columns up front). If you need to swap the byte-order on either input or output you must use upper case S or D instead.
More and more programs are now producing binary data in the netCDF format, and so GMT programs started to support tabular netCDF data (files containing one or more 1-dimensional arrays) starting with GMT version 4.3.0. Because of the meta data contained in those files, reading them is much less complex than reading native binary tables, and even than ASCII tables. GMT programs will read as many 1-dimensional columns as are needed by the program, starting with the first 1-dimensional it can find in the file. To specifically specify which variables are to be read, append the suffix ?var1/var2/... to the netCDF file name or add the option -bicvar1/var2/..., where var1, var2, etc. are the names of the variables to be processed. The latter option is particularly practical when more than one file is read: the -bic option will apply to all files. Currently, GMT only reads, but does not write, netCDF tabular data.
By default, GMT stores 2-D grids as COARDS-compliant netCDF files. COARDS (which stands for Cooperative Ocean/Atmosphere Research Data Service) is a convention used by many agencies distributing gridded data for ocean and atmosphere research. Sticking to this convention allows GMT to read gridded data provided by other institutes and other programs. Conversely, other general domain programs will be able to read grids created by GMT. COARDS is a subset of a more extensive convention for netCDF data called CF-1.0 (Climate and Forecast, version 1.0). Hence, GMT grids are also automatically CF-1.0-compliant. However, since CF-1.0 has more general application than COARDS, not all CF-1.0 compliant netCDF files can be read by GMT.
The netCDF grid file in GMT has several attributes (See Table B.1) to describe the content. The routine that deals with netCDF grid files is sufficiently flexible so that grid files slightly deviating from the standards used by GMT can also be read.
Attribute | Description
|
Global attributes | |
Conventions | COARDS, CF-1.0 (optional) |
title | Title (optional) |
source | How file was created (optional) |
node_offset | 0 for gridline node registration (default), 1 for pixel registration |
- and -variable attributes
| |
long_name | Coordinate name (default: “Longitude” and “Latitude”) |
units | Unit of the coordinate (default: “degrees_east” and “degrees_north”) |
actual_range | Minimum and maximum and of region; if absent |
(or valid_range) | the first and last - and -values are queried |
-variable attributes
| |
long_name | Name of the variable (default: “z”) |
units | Unit of the variable (no default) |
scale_factor | Factor to multiply with (default: 1) |
add_offset | Offset to add to scaled (default: 0) |
actual_range | Minimum and maximum (optional) |
_FillValue | Value associated with missing data points; if absent an |
(or missing_value) | appropriate default value is assumed, depending on data type. |
By default, the first 2-dimensional variable in a netCDF file will by read as the variable and the coordinate axes and will be determined from the dimensions of the variable. GMT will recognize whether the (latitude) variable increases or decreases. Both forms of data storage are handled appropriately.
For more information on the use of COARDS-compliant netCDF files, and on how to load multi-dimensional grids, read Section 4.18.
GMT also allows other formats to be read. In addition to the default netCDF format it can use binary floating points, short integers, bytes, and bits, as well as 8-bit Sun raster files (colormap ignored). Additional formats may be used by supplying read/write functions and linking these with the GMT libraries. The source file gmt_customio.c has the information that programmers will need to augment GMT to read custom grid files. We anticipate that the number of pre-programmed formats will increase as enterprising users implement what they need. See Section 4.17 for more information.
Scanline format means that the data are stored in rows (y = constant) going from the “top” ( (north)) to the “bottom” ( (south)). Data within each row are ordered from “left” ( (west)) to “right” ( (east)). The node_offset signals how the nodes are laid out. The grid is always defined as the intersections of all x ( ) and y ( ) lines. The two scenarios differ in which area each data point represents. The default node registration in GMT is gridline node registration. Most programs can handle both types, and for some programs like grdimage a pixel registered file makes more sense. Utility programs like grdsample and grdproject will allow you to convert from one format to the other; grdedit can make changes to the grid header and convert a pixel- to a gridline-registred grid, or vice versa.
In this registration, the nodes are centered on the grid line intersections and the data points represent the average value in a cell
of dimensions ()
centered on the nodes (Figure B.1). In the case of grid line registration the number of nodes are related to region
and grid spacing by
which for the example in Figure B.1 yields .
Here, the nodes are centered in the grid cells, i.e., the areas between grid lines, and the data points represent the
average values within each cell (Figure B.2. In the case of pixel registration the number of nodes are related to
region and grid spacing by
Thus, given the same region (-R), the pixel node registered grids have one less column and one less row than the grid line registered grids; here we find .
GMT has the option to specify boundary conditions in some programs that operate on grids (grdsample -L; grdgradient -L; grdtrack -L; nearneighbor -L; grdview -L). The boundary conditions come into play when interpolating or computing derivatives near the limits of the region covered by the grid. The default boundary conditions used are those which are “natural” for the boundary of a minimum curvature interpolating surface. If the user knows that the data are periodic in x (and/or y), or that the data cover a sphere with x,y representing longitude,latitude, then there are better choices for the boundary conditions. Periodic conditions on x (and/or y) are chosen by specifying x (and/or y) as the boundary condition flags; global spherical cases are specified using the g (geographical) flag. Behavior of these conditions is as follows:
“Pole conditions” use a 180° phase-shift of the data, requiring 180 modulo .
on the boundary, where is represented by the values in the grid file, and is the derivative in the direction normal to a boundary, and
is the two-dimensional Laplacian operator.
The old style native grid file format that was common in earlier version of GMT is still supported, although the use of netCDF files is strongly recommended. The file starts with a header of 892 bytes containing a number of attributes defining the content. The grdedit utility program will allow you to edit parts of the header of an existing grid file. The attributes listed in Table B.2 are contained within the header record in the order given (except the -array which is not part of the header structure, but makes up the rest of the file). As this header was designed long before 64-bit architectures became available, the jump from the first three integers to the subsequent doubles in the structure does not occur on a 16-byte alignment. While GMT handles the reading of these structures correctly, enterprising programmers must take care to read this header correctly (see our code for details).
Parameter | Description
|
int nx | Number of nodes in the -dimension |
int ny | Number of nodes in the y-dimension |
int node_offset | 0 for grid line registration, 1 for pixel registration |
double x_min | Minimum -value of region |
double x_max | Maximum -value of region |
double y_min | Minimum -value of region |
double y_max | Maximum -value of region |
double z_min | Minimum -value in data set |
double z_max | Maximum -value in data set |
double x_inc | Node spacing in -dimension |
double y_inc | Node spacing in -dimension |
double z_scale_factor | Factor to multiply -values after read |
double z_add_offset | Offset to add to scaled -values |
char x_units[80] | Units of the -dimension |
char y_units[80] | Units of the -dimension |
char z_units[80] | Units of the -dimension |
char title[80] | Descriptive title of the data set |
char command[320] | Command line that produced the grid file |
char remark[160] | Any additional comments |
TYPE z[nx*ny] | 1-D array with -values in scanline format |
The Sun raster file format consists of a header followed by a series of unsigned 1-byte integers that represents the bit-pattern. Bits are scanline oriented, and each row must contain an even number of bytes. The predefined 1-bit patterns in GMT have dimensions of 64 by 64, but other sizes will be accepted when using the -GpP option. The Sun header structure is outline in Table B.3.
Parameter | Description
|
int ras_magic | Magic number |
int ras_width | Width (pixels) of image |
int ras_height | Height (pixels) of image |
int ras_depth | Depth (1, 8, 24, 32 bits) of pixel |
int ras_length | Length (bytes) of image |
int ras_type | Type of file; see RT_* below |
int ras_maptype | Type of colormap; see RMT_* below |
int ras_maplength | Length (bytes) of following map |
After the header, the color map (if ras_maptype is not RMT_NONE) follows for ras_maplength bytes, followed by an image of ras_length bytes. Some related definitions are given in Table B.4.
Macro name | Description
|
RAS_MAGIC | 0x59a66a95 |
RT_STANDARD | 1 (Raw pixrect image in 68000 byte order) |
RT_BYTE_ENCODED | 2 (Run-length compression of bytes) |
RT_FORMAT_RGB | 3 ([X]RGB instead of [X]BGR) |
RMT_NONE | 0 (ras_maplength is expected to be 0) |
RMT_EQUAL_RGB | 1 (red[ras_maplength/3],green[],blue[]) |
Numerous public-domain programs exist, such as xv and convert (in the ImageMagick package), that will translate between various raster file formats such as tiff, gif, jpeg, and Sun raster. Raster patterns may be created with GMT plotting tools by generating PostScript plots that can be rasterized by ghostscript and translated into the right raster format.
Now that you made some nice graphics with GMT, it is time to add them to a document, an article, a report, your dissertation, a poster, a web page, or a presentation. Of course, you could try the old-fashioned scissors and glue stick. More likely, you want to incorporate your graphics electronically into the document. Depending on the application, the GMT PostScript file will need to be converted to Encapsulated PostScript (EPS), Portable Document Format (PDF), or some raster format (e.g., JPEG, PNG, or TIFF) in order to incorporate them into the document.
A large number of questions to the GMT-Help mailing list are related to these rendering issues, showing that something as seemingly straightforward as incorporating a PostScript file into a document is a far from trivial exercise. This Chapter will show how to include GMT graphics into documents and how to achieve the best quality results.
GMT can produce both freeform PostScript files and the more restricted Encapsulated PostScript files (EPS). The former is intended to be sent to a printer or PostScript previewer, while the latter is intended to be included in another document (but should also be able to print and preview). You control what kind of PostScript that GMT produces by manipulating the PAPER_MEDIA parameter (see the gmtdefaults man page for how this is accomplished). Note that a freeform PostScript file may contain special operators (such as Setpagedevice) that is specific to printers (e.g., selection of paper tray). Some previewers (among them, Sun’s pageview) do not understand these valid instructions and may fail to image the file. Also, embedding freeform PostScript with such instructions in it into a larger document can create printing to fail. While you could choose another viewer (we recommend gv (ghostview)) to view single plots prepared by GMT, it is generally wiser anyhow to select EPS output when you are creating a plot intended for inclusion into a larger document. Some programs (and some publishers as well) do not allow the use of instructions like Setpagedevice as part of embedded graphics.
An EPS file that is to be placed into another document needs to have correct bounding box parameters. These are found in the PostScript Document Comment %%BoundingBox. Applications that generate EPS files should set these parameters correctly. Because GMT makes the PostScript files on the fly, often with several overlays, it is not possible to do so accurately. However, GMT does make an effort to ensure that the BoundingBox is large enough to contain the entire composite plot41. Therefore, if you need a “tight” BoundingBox you need to post-process your PostScript file. There are several ways in which this can be accomplished.
should give a tight BoundingBox; epstool assumes the plot is page size and not a huge poster.
should also do the trick. The downside is that this program adds an “image” of the plot in the preamble of the EPS file, thus increasing the file size significantly. This image is a rough rendering of your PostScript graphics that some programs will show on screen while you are editing your document. This image is basically a placeholder for the PostScript graphics that will actually be printed.
will convert the PostScript file myplot.ps into an encapsulated PostScript file myplot.eps which is exactly cropped to the tightest possible BoundingBox.
If you do not want to modify your illustration but just include it in a text document: many word processors (such as Microsoft Word, Corel WordPerfect, and Apple Pages) will let you include a PostScript file that you may place but not edit. Newer versions of those programs also allow you to include PDF versions of your graphics. Except for Pages, you will not be able to view the figure on-screen, but it will print correctly.
Since Adobe’s PDF (Portable Document Format) seems to become the de facto standard for vector graphics, you are often well off converting GMT produced PostScript files to PDF. Being both vector formats (i.e., they basically describe all objects, text and graphics as lines and curves), such conversion sounds awfully straightforward and not worth a full section in this document. But experience has shown differently, since most converters cut corners by using the same tool (Aladdin’s ghostscript) with basic default options that are not devised to produce the best quality PDF files.
For some applications it is practical or even essential that you convert your PostScript file into a raster format, such as GIF (Graphics Interchange Format), TIFF (Tagged Image File Format), PNG (Portable Network Graphics), or JPEG (Joint Photographic Experts Group). A web page is better served with a raster image that will immediately show on a web browser, than with a PostScript file that needs to be downloaded to view, despite the better printing quality of the PostScript image. A less obvious reason to convert your image to a raster format is to by-pass PowerPoint’s rendering engine in case you want to embed the image into a presentation.
The are a number of programs that will convert PostScript files to PDF or raster formats, like Aladdin’s pstopdf, pbmplus’ pstoimg, or ImageMagick’s convert, most of which run ghostscript behind the scenes. The same is true for viewers like gv and Apple’s Preview. So a lot of the times when people report that their PostScript plot does not look right but prints fine, it is the way ghostscript is used with its most basic settings that is to blame.
Here are some notorious pitfalls with ghostscript (and other rendering programs for that matter).
Clearly, this can lead to some unwanted results. First, all edges and lines get blurry and second, the assumption of a white background causes the gray shades to stand out when transferring the image to background with a different color (like the popular sleep-inducing blue in PowerPoint presentations). A more surprising effect of anti-aliasing is that the seams between tiles that make up the land mask when using pscoast will become visible. The anti-aliasing somehow decides to blur the edges of all polygons, even when they are seamlessly connected to other polygons.
It is therefore wise to overrule the default anti-aliasing option and over-sample the image yourself by choosing a higher resolution.
The remedy to all the problems mentioned in the previous section is readily available to you in the form of the GMT utility ps2raster. It is designed to provide the best quality PDF and raster files using ghostscript as a rendering engine. The program ps2raster avoids anti-aliasing and lossy compression techniques that are default to ghostscript and includes the fonts into the resulting PDF file to ensure portability. By default the fonts are rendered at 720 dots-per-inch in a PDF file and images are sampled to 300 dpi, but that can be changed with the -E option. Simply run
to convert all PostScript files to PDF while cropping it to the smallest possible BoundingBox. Or use the -Tg option to convert your files to PNG.
The -P option of ps2raster may also come in handy. When you have not supplied the -P option in your first GMT plot command, your plot will be in Landscape mode. That means that the plot will be rotated 90 degrees (anti-clockwise) to fit on a Portrait mode page when coming out of the printer. The -P option of ps2raster will undo that rotation, so that you do not have to do so within your document. This will only affect Landscape plots; Portrait plots will not be rotated.
Nearly all illustrations in this GMT documentation were GMT-produced PostScript files. They were converted to PDF files using ps2raster and then included into a LATEX document that was processed with pdflatex to create the PDF document you are reading.
To add the graphics into the LATEX document we use the \includegraphics command supplied by the graphicx package. In the preamble of your LATEX document you will need to include the line
The inclusion of the graphics will probably be inside a floating figure environment; something like this
Note that the \includegraphics command does not require you to add the suffix .pdf to the file name. If you run pdflatex, it will look automatically for myplot.pdf. If you run latex, it will use myplot.eps instead.
You can scale your plot using the options width=, height=, or scale=. In addition, if your original graphics was produced in Landscape mode (i.e., you did not use GMT’s -P option: not while plotting, nor in ps2raster), you will need to rotate the plot as well. For example,
will rotate the image 90° clockwise and scale it such that its width (after rotation) will be 80% of the width of the text column.
In Figure C.1 we have attempted to include Figure 7.20 into a PowerPoint presentation. First the PostScript file was converted to PDF (using ps2raster), then loaded into PowerPoint and the white background color was made transparent using the formatting toolbar (shown on the left side of Figure C.1). Clearly, when we let PowerPoint do the rendering, we do not get the best result:
On the central column of Figure C.1 we have included PNG versions of a portion of the same example. This shows the workings of anti-aliasing and different resolutions. All samples were obtained with convert. The one on the top uses all default settings, resulting in an anti-aliased image at 72 dpi resolution (very much like the PDF included directly into PowerPoint).
Just switching anti-aliasing off (middle) is clearly not an option either. It is true that we got rid of the gray blurring and the seams between the tiles, but without anti-aliasing the image becomes very blocky. The solution is to render the image at a higher resolution (e.g., 300 dpi) without anti-aliasing and then shrink the image to the appropriate size (bottom of the central column in Figure C.1). The scaling, rotation as well as the selection of the transparent color can be accomplished through the “Formatting” tool bar and the “Format Picture” dialogue box of PowerPoint (Figure C.2), which can be found by double clicking the included image (or selecting and right-clicking or control-clicking on a one-button mouse).
These examples do not constitute endorsements of the products mentioned above; they only represent our limited experience with adding PostScript to various types of documents. For other solutions and further help, please post messages to gmt-help@lists.hawaii.edu.
All the source code, support data, PDF and HTML versions of all documentation (including UNIX manual pages) can be obtained by anonymous ftp from several mirror sites. We also maintain a GMT page on the World Wide Web (http://gmt.soest.hawaii.edu); see this page for installation directions which allow for a simplified, automatic install procedure (for the purchase of CD-R and DVD-R media, see http://www.geoware-online.com.)
The GMT compressed tar archives requires bzip2 to expand. If this utility is not installed on your system, you must obtained it by your system’s package manager or install it separately42. The GMT archives are as follows:
The netCDF library that makes up the backbone of the grid file i/o operations can be obtained from Unidata by downloading he file netcdf.tar.Z from the anonymous FTP directory of unidata.ucar.edu.
For Windows users who just want executables we have three Windows installers available. Choose one of the first two and optionally the third:
Usually, only one of the GMT 32- or 64-bit installers will be needed.
GMT provides 90 different bit and hachure patterns that can be selected with the -Gp or -GP option in most plotting programs. The left side of each image was created using -Gp, the right side shows the inverted version using -GP. These patterns are reproduced below at 300 dpi.
The characters and their octal codes in the Standard and ISOLatin1 encoded fonts are shown in Figure F.1. Dark gray areas signify codes reserved for control characters. In order to use all the extended characters (shown in the light gray boxes) you need to set CHAR_ENCODING to Standard+ or ISOLatin1+ in your .gmtdefaults4 file43.
The chart for the Symbol character set (GMT font number 12) and Pifont ZapfDingbats character set (font number 34) are presented in Figure F.2 below. The octal code is obtained by appending the column value to the ?? value, e.g., is 266 in the Symbol font. The euro currency symbol is 240 in the Symbol font and will print if your printer supports it (older printer’s firmware will not know about the euro).
GMT uses the standard 35 fonts that come with most PostScript laserwriters. If your printer does not support
some of these fonts, it will automatically substitute the default font (which is usually Courier). The following is a
list of the GMT fonts:
For the special fonts Symbol (12) and ZapfDingbats (34), see the octal charts in Appendix F. When specifying fonts in GMT, you can either give the entire font name or just the font number listed in this table. To change the fonts used in plotting basemap frames, see the man page for gmtdefaults. For direct plotting of text-strings, see the man page for pstext. To add additional fonts that you may have purchased or that are available at your institution, see instructions in the CUSTOM_font_info.d under the share/pslib directory.
GMT creates valid (so far as we know) Adobe PostScript Level 2. It does not use operators specific to Level 3 and should therefore produce output that should print on all PostScript printers44. Sometimes unexpected things happen when GMT output is sent to certain printers or displays. This section lists some things we have learned from experience, and some work-arounds. Note that many of these lessons are now rather old so hopefully these workarounds no longer apply to anybody...
When you try to display a PostScript file on a device, such as a printer or your screen, then a program called a PostScript device driver has to compute which device pixels should receive which colors (black or white in the case of a simple laser printer) in order to display the file. At this stage, certain device-dependent things may happen. These are not limitations of GMT or PostScript, but of the particular display device. The following bugs are known to us based on our experiences:
/setgray {dup dup setrgbcolor} def
immediately following the first line in the file (starts with %!PS.)
The parameter DOTS_PR_INCH can be set by the user through the .gmtdefaults4 file or gmtset. By default it is equal to the value in the gmt_defaults.h file, which is supplied with 300 when you get GMT from us. This seems a good size for most applications, but should ideally reflect the resolution of your hardcopy device (most laserwriters have at least 300 dpi, hence our default value). GMT computes what the plot should look like in double precision floating point coordinates, and then converts these to integer coordinates at DOTS_PR_INCH resolution. This helps us find out that certain points in a path lie on top of other points, and we can remove these, making smaller paths. Small paths are important for the laserwriter bugs above, and also to make fill operations compute faster. Some users have set their DOTS_PR_INCH to very large numbers. This only makes the PostScript output bigger without affecting the appearance of the plot. However, if you want to make a plot which fits on a page at first, and then later magnify this same PostScript file to a huge size, the higher DPI is important. Your data may not have the higher resolution but on certain devices the edges of fonts will not look crisp if they are not drawn with an effective resolution of 300 dpi or so. Beware of making an excessively large path. Note that if you change dpi the linewidths produced by your -W options will change, unless you have appended p for linewidth in points.
Note for users of pageview in Sun OpenWindows: GMT now offers some octal escape sequences to load European alphabet characters in text strings (see Section 4.16). When this feature is enabled, the header on GMT PostScript output includes a section defining special fonts. The definition is added to the header whether or not your plot actually uses the fonts.
Users who view their GMT PostScript output using pageview in OpenWindows on Sun computers or user older laserwriters may have difficulties with the European font definition. If your installation of OpenWindows followed a space-saving suggestion of Sun, you may have excluded the European fonts, in which case pageview will fail to render your plot.
Ask your system administrator about this, or run this simple test: (1) View a GMT PostScript file with
pageview. If it comes up OK, you will be fine. If it comes up blank, open the “Edit PostScript” button and examine
the lower window for error messages. (The European font problem generates lots of error messages in this
window). (2) Verify that the PostScript file is OK, by sending it to a laserwriter and making sure it comes out. (3)
If the PostScript file is OK but it chokes pageview, then edit the PostScript file, cutting out everything between
the lines:
%%%%% START OF EUROPEAN FONT DEFINITION %%%%%
bunch of
definitions
%%%%% END OF EUROPEAN FONT DEFINITION %%%%%
Now try pageview on the edited version. If it now comes up, you have a limited subset of OpenWindows installed. If you discover that these fonts cause you trouble, then you can edit your .gmtdefaults4 file to set CHAR_ENCODING = Standard, which will suppress the printing of this definition in the GMT PostScript header. You can make output which will be viewable in pageview without any editing. However, you would have to reset this to TRUE before attempting to use European fonts, and then the output will become un-pageview-able again. If you try to concatenate segments of GMT PostScript made with and without the European fonts enabled, then you may find that you have problems, either with the definition, or because you ask for something not defined.
When making images and perspective views of large amounts of data, the GMT programs can take some time to run, the resulting PostScript files can be very large, and the time to display the plot can be long. Fine tuning a plot script can take lots of trial and error. We recommend using grdsample to make a low resolution version of the data files you are plotting, and practice with that, so it is faster; when the script is perfect, use the full-resolution data files. We often begin building a script using only psbasemap or pscoast to get the various plots oriented correctly on the page; once this works we replace the psbasemap calls with the actually desired GMT programs.
If you want to make color shaded relief images and you haven’t had much experience with it, here is a good first cut at the problem: Set your COLOR_MODEL to HSV using gmtset. Use makecpt or grd2cpt to make a continuous color palette spanning the range of your data. Use the -Nt option on grdgradient. Try the result, and then play with the tuning of the .gmtdefaults4, the cpt file, and the gradient file.
In this Appendix, we are going to try to explain the relationship between the RGB, CMYK, and HSV color systems so as to (hopefully) make them more intuitive. GMT allows users to specify colors in cpt files in either of these three systems. Interpolation between colors is performed in either RGB or HSV, depending on the specification in the cpt file. Below, we will explain why this all matters.
Remember your (parents’) first color television set? Likely it had three little bright colored squares on it: red, green, and blue. And that is exactly what each color on the tube is made of: varying levels of red, green and blue light. Switch all of them off, , then you have black. All of them at maximum, , creates white. Your computer screen works the same way.
A mix of levels of red, green, and blue creates basically any color imaginable. In GMT each color can be represented by the triplet //. For example, 127/255/0 (half red, full green, and no blue) creates a color called chartreuse. The color sliders in the graphics program GIMP are an excellent way to experiment with colors, since they show you in advance how moving one of the color sliders will change the color. As Figure I.1a shows: increase the red and you will get a more yellow color, while lowering the blue level will turn it into brown.
Is chocolate your favorite color, but you do not know the RGB equivalent values? Then look them up in Figure I.2 or type man gmtcolors for a full list. It’s 210/105/30. But GMT makes it easy on you: you can specify pen, fill, and palette colors by any of the more than 500 unique colors found in that file.
Are you very web-savvy and work best with hexadecimal color codes as they are used in HTML? Even that is allowed in GMT. Just start with a hash mark (#) and follow with the 2 hexadecimal characters for red, green, and blue. For example, you can use #79ff00 for chartreuse, #D2691E for chocolate.
If you have played around with RGB color sliders, you will have noticed that it is not intuitive to make a chosen color lighter or darker, more saturated or more gray. It would involve changing three sliders. To make it easier to manipulate colors in terms of lightness and saturation, another coordinate system was invented: HSV (hue, saturation, value). Those terms can be made clear best by looking at the color sliders in Figure I.1a. Hue (running from 0° to 360°) gives you the full spectrum of saturated colors. Saturation (from 0 to 1, or 100%) tells you how ‘full’ your color is: reduce it to zero and you only have gray scales. Value (from 0 to 1, or 100%) will bring you from black to a fully saturated color. Note that “value” is not the same as “intensity”, or “lightness”, used in other color geometries. “Brilliance” may be the best alternative word to describe “value”. Apple calls it as “brightness”, and hence refers to HSB for this color space.
Want more chartreuse or chocolate? You can specify them in GMT as 90-1-1 and 25-0.86-0.82, respectively.
We are going to try to give you a geometric picture of color mixing in RGB and HSV by means of a tour of the RGB cube depicted in Figure 7.11. The geometric picture is most helpful, we think, since HSV are not orthogonal coordinates and not found from RGB by a simple algebraic transformation. So here goes: Look at the cube face with black, red, magenta, and blue corners. This is the = 0 face. Orient the cube so that you are looking at this face with black in the lower left corner. Now imagine a right-handed cartesian (,,) coordinate system with origin at the black point; you are looking at the plane with increasing to your right, increasing away from you, and increasing up. Keep this sense of (,,) as you look at the cube.
Now tip the cube such that the black corner faces down and the white corner up. When looking from the top, you can see the hue, contoured in gray solid lines, running around in 360° counter-clockwise. It starts with shades of red (0°), then goes through green (120°) and blue (240°), back to red.
On the three faces that are now on the lower side (with the white print) one of (,,) is equal to 0. These three faces meet at the black corner, where . On these three faces the colors are fully saturated: = 1. The dashed white lines indicate different levels of , ranging from 0 to 1 with contours every 0.1.
On the upper three faces (with the black print), one of (,,) is equal to the maximum value. These three faces meet at the white corner, where . On these three faces value is at its maximum: = 1 (or 100%). The dashed black lines indicate varying levels of saturation: ranges from 0 to 1 with contours every 0.1.
Now turn the cube around on its vertical axis (running from the black to the white corner). Along the six edges that zigzag around the “equator”, both saturation and value are maximum, so . Twirling the cube around and tracing the zigzag, you will visit six of the eight corners of the cube, with changing hue (): red (0°), yellow (60°), green (120°), cyan (180°), blue (240°), and magenta (300°). Three of these are the RGB colors; the other three are the CMY colors which are the complement of RGB and are used in many color hardcopy devices (see below). The only cube corners you did not visit on this path are the black and white corners. They lie on the vertical axis where hue is undefined and . Any point on this axis is a shade of gray.
Let us call the points where (points along the RYGCBM path described above) the “pure” colors. If we start at a pure color and we want to whiten it, we can keep constant and while decreasing ; this will move us along one of the cube faces toward the white point. If we start at a pure color and we want to blacken it, we can keep constant and while decreasing ; this will move us along one of the cube faces toward the black point. Any point in (,,) space which can be thought of as a mixture of pure color + white, or pure color + black, is on a face of the cube.
The points in the interior of the cube are a little harder to describe. The definition for above works at all points in (non-gray) (,,) space, but so far we have only looked at (, ) on the cube faces, not inside it. At interior points, none of (,,) is equal to either 0 or 255. Choose such a point, not on the gray axis. Now draw a line through your point so that the line intersects the gray axis and also intersects the RYGCBM path of edges somewhere. It is always possible to construct this line, and all points on this line have the same hue. This construction shows that any point in RGB space can be thought of as a mixture of a pure color plus a shade of gray. If we move along this line away from the gray axis toward the pure color, we are “purifying” the color by “removing gray”; this move increases the color’s saturation. When we get to the point where we cannot remove any more gray, at least one of (,,) will have become zero and the color is now fully saturated; . Conversely, any point on the gray axis is completely undersaturated, so that there. Now we see that the black point is special, is both 0 and 1 at the same time. In other words, at the black point saturation in undefined (and so is hue). The convention is to use at this point.
It remains to define value. To do so, try this: Take your point in RGB space and construct a line through it so that this line goes through the black point; produce this line from black past your point until it hits a face on which . All points on this line have the same hue. Note that this line and the line we made in the previous paragraph are both contained in the plane whose hue is constant. These two lines meet at some arbitrary angle which varies depending on which point you chose. Thus HSV is not an orthogonal coordinate system. If the line you made in the previous paragraph happened to touch the gray axis at the black point, then these two lines are the same line, which is why the black point is special. Now, the line we made in this paragraph illustrates the following: If your chosen point is not already at the end of the line, where , then it is possible to move along the line in that direction so as to increase (,,) while keeping the same hue. The effect this has on a color monitor is to make the color more “brilliant”, your hue will become “stronger”; if you are already on a plane where at least one of (,,) = 255, then you cannot get a stronger version of the same hue. Thus, measures brilliance or strength. Note that it is not quite true to say that measures distance away from the black point, because is not equal to .
Another representation of the HSV space is the color cone illustrated in Figure I.3.
“Pure” colors are around the edge of the circular surface at
the top. Hue runs counter-clockwise. Saturation decreases
to the center. Value increases from zero (black) at the
bottom to 1 at the top. Gray shades are along the vertical
axis.
From studying the RGB cube, we hope you will have understood that there are different routes to follow between two colors, depending whether you are in the RGB or HSV system. Suppose you would make an interpolation between blue and red. In the RGB system you would follow a path diagonally across a face of the cube, from 0/0/255 (blue) via 127/0/127 (purple) to 255/0/0 (red). In the HSV system, you would trace two edges, from 240-1-1 (blue) via 300-1-1 (magenta) to 360-1-1 (red). That is even assuming software would be smart enough to go the shorter route. More likely, red will be recorded as 0-1-1, so hue will be interpolated the other way around, reducing hue from 240° to 0°, via cyan, green, and yellow.
Depending on the design of your color palette, you may want to have it either way. By default, GMT interpolates in RGB space, even when the original color palette is in the HSV system. However, when you add the line #COLOR_MODEL=+HSV (with the leading ‘+’ sign) in the header of the color palette file, GMT will not only read the color representation as HSV values, but also interpolate colors in the HSV system. That means that H, S, and V values are interpolated linearly between two colors, instead of their respective R, G, and B values.
The top row in Figure I.4 illustrates two examples: a blue-white-red scale (the polar palette in Appendix M) interpolated in RGB and the rainbow palette interpolated in HSV. The bottom row of the Figure demonstrates how things can go terribly wrong when you do the interpolation in the other system.
GMT uses the HSV system to achieve artificial illumination of colored images (e.g., -I option in grdimage) by changing the saturation s and value v coordinates of the color. When the intensity is zero (flat illumination), the data are colored according to the cpt file. If the intensity is non-zero, the color is either lightened or darkened depending on the illumination. The color is first converted to HSV (if necessary) and then darkened by moving (,) toward (HSV_MIN_SATURATION, HSV_MIN_VALUE) if the intensity is negative, or lightened by sliding (,) toward (HSV_MAX_SATURATION, HSV_MAX_VALUE) if the illumination is positive. The extremes of the and are defined in the .gmtdefaults4 file and are usually chosen so the corresponding points are nearly black ( = 1, = 0) and white ( = 0, = 1). The reason this works is that the HSV system allows movements in color space which correspond more closely to what we mean by “tint” and “shade”; an instruction like “add white” is easy in HSV and not so obvious in RGB.
The RGB system is understandable because it is cartesian, and we all learned cartesian coordinates in school. But it doesn’t help us create a tint or shade of a color; we cannot say, “We want orange, and a lighter shade of orange, or a less vivid orange”. With HSV we can do this, by saying, “Orange must be between red and yellow, so its hue is about = 30°; a less vivid orange has a lesser , a darker orange has a lesser ”. On the other hand, the HSV system is a peculiar geometric construction, more like a cone (Figure I.3). It is not an orthogonal coordinate system, and it is not found by a matrix transformation of RGB; these make it difficult in some cases too. Note that a move toward black or a move toward white will change both and , in the general case of an interior point in the cube. The HSV system also doesn’t behave well for very dark colors, where the gray point is near black and the two lines we constructed above are almost parallel. If you are trying to create nice colors for drawing chocolates, for example, you may be better off guessing in RGB coordinates.
Finally, you can imagine that printers work in a different way: they mix different paints to make a color. The more paint, the darker the color, which is the reverse of adding more light. Also, mixing more colored paints does not give you true black, so that means that you really need four colors to do it right. Open up your color printer and you’ll probably find four cartridges: cyan, magenta, yellow (often these are combined into one), and black. They form the CMYK system of colors, each value running from 0 to 1 (or 100%). In GMT CMYK color coding can be achieved using /// quadruplets.
Obviously, there is no unique way to go from the 3-dimensional RGB system to the 4-dimensional CMYK system. So, again, there is a lot of hand waving applied in the transformation. Strikingly, CMYK actually covers a smaller color space than RGB. We will not try to explain you the details behind it, just know that there is a transformation needed to go from the colors on your screen to the colors on your printer. It might explain why what you see is not necessarily what you get. If you are really concerned about how your color plots will show up in your PhD thesis, for example, it might be worth trying to save and print all your color plots using the CMYK system. Letting GMT do the conversion to CMYK may avoid some nasty surprises when it comes down to printing. To specify the color space of your PostScript file, set PS_COLOR in the .gmtdefaults4 file to RGB, HSV, or CMYK.
The GMT programs filter1d (for tables of data indexed to one independent variable) and grdfilter (for data given as 2-dimensional grids) allow filtering of data by a moving-window process. (To filter a grid by Fourier transform use grdfft.) Both programs use an argument -Ftypewidth to specify the type of process and the window’s width (in 1-d) or diameter (in 2-d). (In filter1d the width is a length of the time or space ordinate axis, while in grdfilter it is the diameter of a circular area whose distance unit is related to the grid mesh via the -D option). If the process is a median, mode, or extreme value estimator then the window output cannot be written as a convolution and the filtering operation is not a linear operator. If the process is a weighted average, as in the boxcar, cosine, and gaussian filter types, then linear operator theory applies to the filtering process. These three filters can be described as convolutions with an impulse response function, and their transfer functions can be used to describe how they alter components in the input as a function of wavelength.
Impulse responses are shown here for the boxcar, cosine, and gaussian filters. Only the relative amplitudes of the filter weights shown; the values in the center of the window have been fixed equal to 1 for ease of plotting. In this way the same graph can serve to illustrate both the 1-d and 2-d impulse responses; in the 2-d case this plot is a diametrical cross-section through the filter weights (Figure J.1).
Although the impulse responses look the same in 1-d and 2-d, this is not true of the transfer functions; in 1-d the transfer function is the Fourier transform of the impulse response, while in 2-d it is the Hankel transform of the impulse response. These are shown in Figures J.2 and J.3, respectively. Note that in 1-d the boxcar transfer function has its first zero crossing at , while in 2-d it is around . The 1-d cosine transfer function has its first zero crossing at ; so a cosine filter needs to be twice as wide as a boxcar filter in order to zero the same lowest frequency. As a general rule, the cosine and gaussian filters are “better” in the sense that they do not have the “side lobes” (large-amplitude oscillations in the transfer function) that the boxcar filter has. However, they are correspondingly “worse” in the sense that they require more work (doubling the width to achieve the same cut-off wavelength).
One of the nice things about the gaussian filter is that its transfer functions are the same in 1-d and 2-d. Another nice property is that it has no negative side lobes. There are many definitions of the gaussian filter in the literature (see page 7 of Bracewell45). We define equal to 1/6 of the filter width, and the impulse response proportional to . With this definition, the transfer function is and the wavelength at which the transfer function equals 0.5 is about 5.34 , or about 0.89 of the filter width.
Starting with version 3.0, GMT use a completely new coastline database and the pscoast utility was been completely rewritten to handle the new file format. Many users have asked us why it has taken so long for GMT to use a high-resolution coastline database; after all, such data have been available in the public domain for years. To answer such questions we will take you along the road that starts with these public domain data sets and ends up with the database used by GMT.
There are two well-known public-domain data sets that could be used for this purpose. Once is known as the World Data Bank II or CIA Data Bank (WDB) and contains coastlines, lakes, political boundaries, and rivers. The other, the World Vector Shoreline (WVS) only contains shorelines between saltwater and land (i.e., no lakes). It turns out that the WVS data is far superior to the WDB data as far as data quality goes, but as noted it lacks lakes, not to mention rivers and borders. We decided to use the WVS whenever possible and supplement it with WDB data. We got these data over the Internet; they are also available on CD-ROM from the National Geophysical Data Center in Boulder, Colorado46.
In order to paint continents or oceans it is necessary that the coastline data be organized in polygons that may be filled. Simple line segments can be used to draw the coastline, but for painting polygons are required. Both the WVS and WDB data consists of unsorted line segments: there is no information included that tells you which segments belong to the same polygon (e.g., Australia should be one large polygon). In addition, polygons enclosing land must be differentiated from polygons enclosing lakes since they will need different paint. Finally, we want pscoast to be flexible enough that it can paint the land or the oceans or both. If just land (or oceans) is selected we do not want to paint those areas that are not land (or oceans) since previous plot programs may have drawn in those areas. Thus, we will need to combine polygons into new polygons that lend themselves to fill land (or oceans) only (Note that older versions of pscoast always painted lakes and wiped out whatever was plotted beneath).
The WVS and WDB together represent more than 100 Mb of binary data and something like 20 million data points. Hence, it becomes obvious that any manipulation of these data must be automated. For instance, the reasonable requirement that no coastline should cross another coastline becomes a complicated processing step.
GMT uses a polygon-assembly routine that carries out these tasks on the fly.
We will demonstrate the power of the new database by starting with a regional hemisphere map centered near Papua New Guinea and zoom in on a specified point. The map regions will be specified in projected km from the projection center, e.g., we may want the map to go from -2000 km to +2000 km in the longitudinal and the latitudinal direction. However, GMT programs expects degrees in the -R option that specifies the desired region. Given the chosen map projection we can automate this process by using a simple shell function that we call getbox:
_________________________________________________________________________________
Also, as we zoom in on the projection center we want to draw the outline of the next map region on the plot. To do that we need the geographical coordinates of the four corners of the region rectangle. Again, we automate this task by using our simple function getrect:
_________________________________________________________________________________
We begin with an azimuthal equidistant map of the hemisphere centered on 130°21’E, 0°12’S, which is slightly west of New Guinea, near the Strait of Dampier. The edges of the map are all 9000 km true distance from the projection center. At this scale (and for global maps) the crude resolution data will usually be adequate to capture the main geographic features. To avoid cluttering the map with insignificant detail we only plot features (i.e., polygons) that exceed 500 km in area. Smaller features would only occupy a few pixels on the plot and make the map look “dirty”. We also add national borders to the plot. The crude database is heavily decimated and simplified by the DP-routine: The total file size of the coastlines, rivers, and borders database is only 283 kbytes. The plot is produced by the script:
_________________________________________________________________________________
Here, we use the OBLIQUE_ANNOTATION bit flags to achieve horizontal annotations and set ANNOT_MIN_SPACING to suppress some longitudinal annotations near the S pole that otherwise would overprint. The box indicates the outline of the next map.
We have now reduced the map area by zooming in on the map center. Now, the edges of the map are all 2000 km true distance from the projection center. At this scale we choose the low resolution data that faithfully reproduce the dominant geographic features in the region. We cut back on minor features less than 100 km in area. We still add national borders to the plot. The low database is less decimated and simplified by the DP-routine: The total file size of the coastlines, rivers, and borders combined grows to 907 kbytes; it is the default resolution in GMT. The plot is generated by the script:
_________________________________________________________________________________
We continue to zoom in on the map center. In this map, the edges of the map are all 500 km true distance from the projection center. We abandon the low resolution data set as it would look too jagged at this scale and instead employ the intermediate resolution data that faithfully reproduce the dominant geographic features in the region. This time, we ignore features less than 20 km in area. Although the script still asks for national borders none exist within our region. The intermediate database is moderately decimated and simplified by the DP-routine: The combined file size of the coastlines, rivers, and borders now exceeds 3.35 Mbytes. The plot is generated by the script:
_________________________________________________________________________________
The relentless zooming continues! Now, the edges of the map are all 100 km true distance from the projection center. We step up to the high resolution data set as it is needed to accurately portray the detailed geographic features within the region. Because of the small scale we only ignore features less than 1 km in area. The high resolution database has undergone minor decimation and simplification by the DP-routine: The combined file size of the coastlines, rivers, and borders now swells to 12.3 Mbytes. The map and the final outline box are generated by these commands:
_________________________________________________________________________________
We now arrive at our final plot, which shows a detailed view of the western side of the small island of Waigeo. The map area is approximately 40 by 40 km. We call upon the full resolution data set to portray the richness of geographic detail within this region; no features are ignored. The full resolution has undergone no decimation and it shows: The combined file size of the coastlines, rivers, and borders totals a (once considered hefty) 55.9 Mbytes. Our final map is reproduced by the single command:
_________________________________________________________________________________
We hope you will study these examples to enable you to make efficient and wise use of this vast data set.
While GMT can be ported to non-UNIX systems such as Windows, it is also true that one of the strengths of GMT lies its symbiotic relationship with UNIX. We therefore recommend that GMT be installed in a POSIX-compliant UNIX environment such as traditional UNIX-systems, Linux, or Mac OS X. If abandoning your non-UNIX operating system is not an option, consider one of these solutions:
Because GMT works best in conjugation with UNIX tools we suggest you install GMT using the Cygwin product from Cygnus (now assimilated by Redhat, Inc.). This free version works on any Windows version and it comes with both the Bourne Again shell bash and the tcsh. You also have access to most standard GNU development tools such as compilers and text processing tools (awk, grep, sed, etc.). Note that executables prepared for Windows will also run under Cygwin.
Follow the instructions on the Cygwin page51 on how to install the package; note you must explicitly add all the development tool packages (e.g., gcc etc) as the basic installation does not include them by default. Once you are up and running under Cygwin, you may install GMT the same way you do under any other UNIX platform by either running the automated install via install_gmt4.sh or manually running configure first, then type make all. If you install via the web form, make sure you save the parameter file without DOS CR/LF endings. Use dos2unix to get rid of those if need be.
Finally, from Cygwin’s User Guide: By default, no Cygwin program can allocate more than 384 MB of memory (program and data). You should not need to change this default in most circumstances. However, if you need to use more real or virtual memory in your machine you may add an entry in either the HKEY_LOCAL_MACHINE (to change the limit for all users) or HKEY_CURRENT_USER (for just the current user) section of the registry. Add the DWORD value heap_chunk_in_mb and set it to the desired memory limit in decimal Mb. It is preferred to do this in Cygwin using the regtool program included in the Cygwin package. (For more information about regtool or the other Cygwin utilities, see the Section called Cygwin Utilities in Chapter 3 of the Cygwin’s User Guide or use the help option of each utility.) You should always be careful when using regtool since damaging your system registry can result in an unusable system. This example sets the local machine memory limit to 1024 Mb:
For more installation details see the general README file.
SFU52 is also similar to Cygwin in that it provides precompiled UNIX tools for DOS/WIN32, including the sh and csh shells.
DJGPP53 is similar to Cygwin in that it provides precompiled UNIX tools for DOS/WIN32, including the bash shell. At the time of this writing we have not been successful in compiling netCDF in this environment. This is fully due to our limited understanding of the innards of the netCDF installation whose configure script did not work for us. As soon as this problem is overcome we expect a smooth install similar to that of Cygwin.
GMT will compile and install using the Microsoft Visual C/C++ compiler. We expect other WIN32 C compilers to give similar results. Since configure cannot be run you must manually rename gmt_notposix.h.in to gmt_notposix.h. The netCDF home page gives full information on how to compile and install netCDF; precompiled libraries are also available. At present we simply have a lame gmtinstall.bat file that compiles the entire GMT package, and gmtsuppl.bat which compiles most of the supplemental programs. If you just need to run GMT and do not want to mess with compilations, get the precompiled binaries from the GMT ftp sites.
GMT has been ported to OS/2 by Allen Cogbill, Los Alamos National Laboratory. One must have EMX installed in order to use the executables. All features that are present in the UNIX version of GMT are available in the OS/2 version. All executables may be obtained using links in the following document, which provides more detail on the port.
GMT will install directly under Mac OS X since it is fully Unix compliant.
Figure M.1 shows each of the 22 built-in color palettes, stored in so-called cpt tables. The programs makecpt and grd2cpt are used to access these master cpt tables and translate/scale them to fit the user’s range of -values. The top half of the color bars in the Figure shows the original color scale, which can be either discrete or continuous, though some (like globe) are a mix of the two. The bottom half the color bar are built by using makecpt -T-1/1/0.25, thus splitting the color scale into 8 discrete colors.
The use of color legends has already been introduced in Chapter 7 (examples 2, 16, and 17). Things become a bit more complicated when you want to label the legend with names for certain intervals (like geological time periods in the example below). To accomplish that, one should add a semi-colon and the label name at the end of a line in the cpt table and add the -L option to the psscale command that draws the color legend. This option also makes all intervals in the legend of equal length, even it the numerical values are not equally spaced.
Normally, the name labels are plotted at the lower end of the intervals. But by adding a gap amount (even when zero) to the -L option, they are centered. The example below also shows how to annotate ranges using -Li (in which case no name labels should appear in the cpt file), and how to switch the color bar around (by using a negative length).
_________________________________________________________________________________
GMT comes with several custom plot symbols ready to go. They are used in psxy and psxyz using the -Sk
option. To make your own custom plot symbol, please follow the instructions given in the man pages of those two
programs. The following is a plot of each symbol. Note that we only show the symbol outline and not any fill. Be
aware that some symbols may have a hardwired fill or no-fill component. Also note that some symbols,
in particular the geometric ones, duplicate what is already available as standard built-in symbols.
The GMT programs grdcontour (for data given as 2-dimensional grids) and pscontour (for x,y,z tables) allow for contouring of data sets, while psxy and psxyz can plot lines based on x,y- and x,y,z-tables, respectively. In both cases it may be necessary to attach labels to these lines. Clever or optimal placements of labels is a very difficult topic, and GMT provides several algorithms for this placement as well as complete freedom in specifying the attributes of the labels. Because of the richness of these choices we present this Appendix which summarizes the various options and gives several examples of their use.
While the previous GMT versions 1–3 allowed for a single algorithm that determined where labels would be placed, GMT 4 allows for five different algorithms. Furthermore, a new “symbol” option (-Sq for “quoted line”) has been added to psxy and psxyz and hence the new label placement mechanisms apply to those programs as well. The contouring programs expect the algorithm to be specified as arguments to -G while the line plotting programs expect the same arguments to follow -Sq. The information appended to these options is the same in both cases and is of the form [code]info. The five algorithms correspond to the five codes below (some codes will appear in both upper and lower case; they share the same algorithm but differ in some other ways). In what follows, the phrase “line segment” is taken to mean either a contour or a line to be labeled. The codes are:
Only one algorithm can be specified at any given time.
Determining where to place labels is half the battle. The other half is to specify exactly what are the attributes of the labels. It turns out that there are quite a few possible attributes that we may want to control, hence understanding how to specify these attributes becomes important. In the contouring programs, one or more attributes may be appended to the -A option using the format +code[args] for each attribute, whereas for the line plotting programs these attributes are appended to the -Sq option following a colon (:) that separates the label codes from the placement algorithm. Several of the attributes do not apply to contours so we start off with listing those that apply universally. These codes are:
For contours, the label will be the value of the contour (possibly modified by +u or +=). However, for quoted lines other options apply:
We will demonstrate the use of these options with a few simple examples. First, we will contour a subset of the global geoid data used in GMT Example 01; the region selected encompasses the world’s strongest “geoid dipole”: the Indian Low and the New Guinea High.
Our first example uses the default placement algorithm. Because of the size of the map we request contour labels every 1.5 inches along the lines:
_________________________________________________________________________________
As seen in Figure O.1, the contours are placed rather arbitrary. The string of contours for to align well but that is a fortuitous consequence of reaching the 1.5 inch distance from the start at the bottom of the map.
We now exercise the option for specifying exactly how many labels each contour line should have:
_________________________________________________________________________________
By selecting only one label per contour and requiring that labels only be placed on contour lines whose length exceed 1 inch, we achieve the effect shown in Figure O.2.
Here, we specify four points where we would like contour labels to be placed. Our points are not exactly on the contour lines so we give a nonzero “slop” to be used in the distance calculations: The point on the contour closest to our fixed points and within the given maximum distance will host the label.
_________________________________________________________________________________
The angle of the label is evaluated from the contour line geometry, and the final result is shown in Figure O.3.
To aid in understanding the algorithm we chose to specify “debug” mode (+d) which placed a small circle at each of the fixed points.
Often, it will suffice to place contours at the imaginary intersections between the contour lines and a well-placed straight line segment. The -Gl or -GL algorithms work well in those cases:
_________________________________________________________________________________
The obvious choice in this example is to specify a great circle between the high and the low, thus placing all labels between these extrema.
The thin debug line in Figure O.4 shows the great circle and the intersections where labels are plotted. Note that any number of such lines could be specified; here we are content with just one.
If (1) the number of intersecting straight line segments needed to pick the desired label positions becomes too large to be given conveniently on the command line, or (2) we have another data set or lines whose intersections we wish to use, the general crossing algorithm makes more sense:
_________________________________________________________________________________
In this case, we have created three strands of lines whose intersections with the contours define the label placements, presented in Figure O.5.
We will now demonstrate some of the ways to play with the label attributes. To do so we will use psxy on a great-circle line connecting the geoid extrema, along which we have sampled the ETOPO5 relief data set. The file transect.d thus contains lon, lat, dist, geoid, relief, with distances in km.
This example will change the orientation of labels from along-track to across-track, and surrounds the labels with an opaque, outlined textbox so that the label is more readable. We choose the place the labels every 1000 km along the line and use that distance as the label. The labels are placed normal to the line:
_________________________________________________________________________________
The composite illustration in Figure O.6 shows the new effects. Note that the line connecting the extrema does not end exactly at the ‘-’ and ‘+’ symbols. This is because the placements of those symbols are based on the mean coordinates of the contour and not the locations of the (local or global) extrema.
A small variation on this theme is to place the labels parallel to the line, use spherical degrees for placement, append the degree symbol as a unit for the labels, choose a rounded rectangular textbox, and inverse-video the label:
_________________________________________________________________________________
The output is presented as Figure O.7.
In the next example we will use the bathymetry values along the transect as our label, with placement determined by the distance along track. We choose to place labels every 1500 km. To do this we need to pull out those records whose distances are multiples of 1500 km and create a “fixed points” file that can be used to place labels and specify the labels. This is done with awk.
_________________________________________________________________________________
The output is presented as Figure O.8.
Finally, we will make a more complex composite illustration that uses several of the label placement and label attribute settings discussed in the previous sections. We make a map showing the tsunami travel times (in hours) from a hypothetical catastrophic landslide in the Canary Islands55. We lay down a color map based on the travel times and the shape of the seafloor, and travel time contours with curved labels as well as a few quoted lines. The final script is
_________________________________________________________________________________
with the complete illustration presented as Figure O.9.
In Chapter 4 it is described how GMT creates several (temporary) files to communicate between the different commands that make up the script that finally creates a plot. Among those files are:
A cure to all these woes is the isolation mode introduced in GMT version 4.2.2. This mode allows you to run a GMT script without leaving any traces other than the resulting PostScript or data files, and not altering the .gmtdefaults4 or .gmtcommands4 files. Those files will be placed in a temporary directory instead. And if properly set up, this temporary directory will only be used by a single script, even if another GMT script is running simultaneously. This also provides the opportunity to create any other temporary files that the script might create in the same directory.
The example below shows how isolation mode works.
_________________________________________________________________________________
The files .gmtdefaults4 and .gmtcommands4 are automatically created in the temporary directory $GMT_TMPDIR. The script is also adjusted such that the temporary grid file lat.nc and colormap lat.cpt are created in that directory as well. To make things even more easy, GMT now provides a set of handy shell functions in gmt_shell_functions.sh: simply include that file in the script and the creation and the removal of the temporary directory is reduced to a single command.
_________________________________________________________________________________
We encourage all GMT users to start using version 4 immediately; it has been tested extensively by the GMT team and has benefitted from bug reports for the 3.4.x versions. Users who still worry about the new version breaking things may install GMT 3.4.x versions and 4.x and use our utility gmtswitch to select their current version should the need to switch arises. You will find gmtswitch in the top-level GMT4.x directory; install as explained below.
Because GMT 4.x is backwards compatible with the 3.4.x series yet maintains its parameters and history in separate hidden files (e.g., .gmtdefaults4 versus .gmtdefaults) it is possible to install and use both versions on the same workstation. To simplify such setups we supply the utility gmtswitch which simplifies switching back and forth between any number of installed GMT 3-versions and GMT 4.x. Place the gmtswitch Bourne shell script in your general executable path (not in one of the GMT bin directories) and run it after you have finished installing all GMT versions of interest. The first time you run gmtswitch it will try to find all the available versions installed on your file system. The versions found will be listed in the file .gmtversions in your home directory; each line is the full path to a GMT root directory (e.g., /usr/local/GMT3.4.2). You may manually add or remove entries there at any time. You are then instructed to make two changes to your environment (the details are shell-dependent but explained by gmtswitch):
Make those edits, logout, and log and back in again. The next time you run gmtswitch you will be able to switch between versions. Typing gmtswitch with no argument will list the available versions in a numerical menu and prompt you to choose one, whereas gmtswitch version will immediately switch to that version (version must be a piece of unique text making up the full path to a version, e.g., 3.4.2). If you use tcsh or csh you may have to type “rehash” to initiate the path changes.
1Version 1.0 was then informally released at the Lamont-Doherty Earth Observatory.
2Use standard UNIX tools such as awk or perl to reformat files should your date and clock components reside in separate columns.
3Used to color the background, foreground, and Not-a-Number areas.
4See GNU General Public License (www.gnu.org/copyleft/gpl.html) for terms on redistribution and modifications.
5The tools can also be installed on other platforms (see Appendix L).
6One public-domain RIP is ghostscript, available from www.gnu.org.
7Programs now also allow for fast, binary multicolumn file i/o.
8While the netCDF format is the default, other formats are also possible, including user-defined formats.
9PostScript definition. In the typesetting industry a slightly different definition of point (1/72.27 inch) is used.
10Choose between SI and US default units by modifying gmt_setup.conf in the GMT share directory.
11To remain backwards compatible with GMT 3.4.x we will also look for .gmtdefaults but only if .gmtdefaults4 cannot be found.
12The Gregorian Calendar is a revision of the Julian Calendar which was instituted in a papal bull by Pope Gregory XIII in 1582. The reason for the calendar change was to correct for drift in the dates of significant religious observations (primarily Easter) and to prevent further drift in the dates. The important effects of the change were (a) Drop 10 days from October 1582 to realign the Vernal Equinox with 21 March, (b) change leap year selection so that not all years ending in “00” are leap years, and (c) change the beginning of the year to 1 January from 25 March. Adoption of the new calendar was essentially immediate within Catholic countries. In the Protestant countries, where papal authority was neither recognized not appreciated, adoption came more slowly. England finally adopted the new calendar in 1752, with eleven days removed from September. The additional day came because the old and new calendars disagreed on whether 1700 was a leap year, so the Julian calendar had to be adjusted by one more day.
13While UTM coordinates clearly refer to points on the Earth, in this context they are considered “other”. Thus, when we refer to “geographical” coordinates herein we imply longitude, latitude.
14However, it is suppressed when a 3-D view is selected.
15Please consult the man page for printf or any book on C .
16For historical reasons, the GMT Default is Landscape, see gmtdefaults to change this.
17Ensures that boundary annotations do not fall off the page.
18For an overview of color systems such as HSV, see Appendix I.
19Convert other graphics formats to Sun ras format using ImageMagick’s convert program.
20For CMYK the format obviously involves two extra columns.
21Snyder, J. P., 1987, Map Projections A Working Manual, U.S. Geological Survey Prof. Paper 1395.
22This is, however, not the shortest distance. It is given by the great circle connecting the two points.
23Robinson provided a table of -coordinates for latitudes every 5°. To project values for intermediate latitudes one must interpolate the table. Different interpolants may result in slightly different maps. GMT uses the interpolant selected by the parameter INTERPOLANT in the .gmtdefaults4 file.
24These data are available on CD-ROM from NGDC (www.ngdc.noaa.gov).
25These data are available on CD-ROM from NGDC (www.ngdc.noaa.gov).
26See http://topex.ucsd.edu/marine_grav/mar_grav.html.
27You can also use the utility curl
28Pedants who wish to argue about the “other” arc going the long way should consider using it.
29You could also use img2mercgrd directly – your only option under DOS
30While QuickTime is free you must upgrade to QuickTime Pro (USD 30) to use the authoring functions.
31QuickTime Pro can do this, as can most video-editing programs.
32Simon.Cox@csiro.au
33For data bases, see http://topex.ucsd.edu/marine_grav/mar_grav.html.
34Walter.HF.Smith@noaa.gov
35Kurt.Feigl@cnes.fr
36patau@ipgp.jussieu.fr
37The ASCII MGD77 data are available on CD-ROM from NGDC (www.ngdc.noaa.gov).
38These data are available on CD-ROM from NGDC (www.ngdc.noaa.gov).
39Timothy.J.Henstock@soc.soton.ac.uk
40lloyd@must-have-coffee.gen.nz
41In contrast, regular GMT PostScript files simply have a %%BoundingBox that equal the size of the chosen paper.
42http://www.bzip.org
43If you chose SI units during the installation then the default encoding is ISOLatin1+, otherwise it is Standard+.
44Note, however, that the -Q option in grdimage will exercise a PostScript Level 3 feature called colormasking.
45R. Bracewell, The Fourier Transform and its Applications, McGraw-Hill, London, 444p., 1965.
46www.ngdc.noaa.gov
47Douglas, D.H., and T. K. Peucker, 1973, Algorithms for the reduction of the number of points required to represent a digitized line or its caricature, Canadian Cartographer, 10, 112–122.
48The full and high resolution files are in separate archives because of their size. Not all users may need these files as the intermediate data set is better than the data provided with version 2.1.4.
49If you need complete polygons in a simpler format, see the article on GSHHS (Wessel, P., and W. H. F. Smith, 1996, A Global, self-consistent, hierarchical, high-resolution shoreline database, J. Geophys. Res. 101, 8741–8743).
50Microsoft Services for UNIX is formerly known as Interix, in the distant past known as OpenNT.
51cygwin.com
52See www.microsoft.com/technet/interopmigration/unix/sfu for details.
53See www.gnu.org for details.
54or whatever planet we are dealing with.
55Travel times were calculated using Geoware’s travel time calculator, ttt; see (http://www.geoware-online.com)