GMT Notes

Some Pointers on Plotting in GMT

Plotting time series
Example 1
This example plots from year-month-day and time of day:

# set default time formats
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm # this could also be hh:mm:ss
gmtset PLOT_DATE_FORMAT dd-o #dd is day, o is month as text

set PROJ = -JX5iT/5i # Append a T to denote a time axis
# Time coordinates are in the format you gave above. Clock is optional, in which case the 00:00 and 23:59 would be removed.
set COORD = -R1998-07-15T00:00/1998-08-15T23:59/0/30
set INCREMX = 1U # In this case, I use the ISO week format to annotate every one week following PLOT_DATE_FORMAT
set INCREMY = 5

# My input file is of the format yyyy-mm-ddThh:mm data (i.e. one row might be 1998-12-31T112:30 10)
# if using a 12 hour clock, instead of military, I would need to append AM or PM after my time
psxy IN $COORD $PROJ -Y5i -W4,255/0/0 -Ba$INCREMX\:"Day":/$INCREMY\:"Temperature, deg. C":WSne -K -P >! $OUT

Example 2
This example reads in data in the format day of year:

gmtset INPUT_DATE_FORMAT yyyy-jjj # input format year-day_of_year
gmtset PLOT_DATE_FORMAT dd-o # output format again day-month as text
gmtset TIME_FORMAT_PRIMARY abbreviated # this time abbreviate the month

set PROJ = -JX5iT/5i # include T for time axis
set COORD = -R1998-083T/1998-272T/0/30 # Now the time coordinates are year-day_of_year
set INCREMX = 4U # I'm plotting at every 4 weeks
set INCREMY = 5

# This time my input data is year day_of_year data
awk '{if(NR>2)printf("%d-%03d %f\n",$1,$2,$3)}' $IN | psxy $COORD $PROJ -Y5i -W4,255/0/0 -Ba$INCREMX\:"Day":/$INCREMY\:"Temperature, deg. C":WSne -K -P >! $OUT

Check GMT website under the psbasemap man page for other PLOT_DATE_FORMAT options

Useful for plotting polygons (include -M -L in psxy)
// reverse the line order of a file
awk ' { line[NR] = $0 }
END { for (i=NR;i>0;i--)
print line[i] }' listlogs

Histograms
In addition to the pshistogram option, I like to use the psxy command with the -Sb1 option. For instance: awk '{print $1,$2}' $IN | psxy $PROJ $COORD -G255/0/0 -W0/0/0 -Sb1 -K -P -O >> $OUT which will fill (-G) in bars that are 1 point wide with red and outline (-W) with black. This requires that you have the appropriate x, y values.

Plotting vector lines
To plot vector lines that are color-coded, use psxy with a multiple segment file for your input. The header for each line of segment is -Zval, where val=the value corresponding to your color in the CPT file. Use -M to denote that you are plotting multiple segments. An example segment file looks like this:
> -Z24.200000
-113.690000 44.560000
-113.690000 44.440000
-113.690000 44.440000
-113.560000 44.440000

Issues with CPT files?
Jenny Adam sent the following email awhile back:
If anybody is having trouble plotting in GMT - and the error has to do with the cpt file, e.g. the error message is something like:

grdimage: GMT Fatal Error: Error when decoding diff.cpt - aborts!
psscale: GMT Fatal Error: Error when decoding diff.cpt - aborts!

I think Fengge was this first person to report this problem last summer sometime. This means that you are using a GMT script that has been handed down from person to person (I believe it was originally written by Ed Maurer) that has an awk statement following the makecpt statement that looks similar to this:

makecpt -T$BAR -Crainbow -Z | \
awk '{if(NR==4){a=$2;b=$3;c=$4} \
if($1!="B" && $1!="F" && $1!="N"){d=$6;e=$7;f=$8} \
if($1=="B"){$2=a;$3=b;$4=c} \
if($1=="F"){$2=d;$3=e;$4=f} \
if($1=="N"){$2=255;$3=255;$4=255} \
print $0 }' >! $CPT

What the awk statement does is it sets the colors of the values that either are smaller than or larger than the defined color range ($BAR) - it sets the colors to be the same as the min and max colors (e.g. "B" and "F"). This is handy because you don't have odd-looking white and black clumps in your spatial plot. It also sets the color of nodata values ("N") to be white (255/255/255). The problem now is that GMT seems to be working only in a HSV rather than a RGB color scheme, and the HSV color range only accepts values from 0 to 1 for the saturation and value components. Therefore, what you need to do is change the second to last line of this script to look like this:

makecpt -T$BAR -Crainbow -Z | \
awk '{if(NR==4){a=$2;b=$3;c=$4} \
if($1!="B" && $1!="F" && $1!="N"){d=$6;e=$7;f=$8} \
if($1=="B"){$2=a;$3=b;$4=c} \
if($1=="F"){$2=d;$3=e;$4=f} \
if($1=="N"){$2=0;$3=0;$4=1} \
print $0 }' >! $CPT

In HSV, 0/0/1 gives white and 0/0/0 gives black. This should solve the problem.

Miscellaneous
From Ben Livneh:
In order to widen the gmt frame, you go into the .gmtdefaults file and on the line that reads "PAPER_MEDIA" you can set it to a custom size as follows:
PAPER_MEDIA = Custom_x
e.g. For a wide plot
PAPER_MEDIA = Custom_5000x492
Now, I believe you may be able to accomplish this from within your gmt script by typing something like
gmtset PAPER_MEDIA Custom_5000x492

Converting figures when you are done!
There are many ways to do this:

convert -trim -antialias -density 300 FILENAME.ps FILENAME.png
Use density 300 for web posting, 600 for publication.
Use gsview.
- gsview <.ps file>
- in menu FILE select PS TO EPS
- DO NOT select "automatically calculate bounding box" because the problem from ps2eps.scr will be the same here, some part of your figure will be cropped, then OK
- Manually with the mouse, click on each of the four boundaries ( left, bottom, right, top)
- rename the file to .eps
- DONE
- check your file. Note that the size between the .ps and the .eps should be about the same.
Use gsview to create a pdf file .... follow these steps
Get a hold of the in-house ps2eps.scr script (ask around the group) that sets BoundingBox and setpagedevice of the .ps file, which is the only difference between ps and eps. There have been some issues lately, so use with caution.

Last updated 7/23/2010 by Elizabeth Clark