Maps binbash psbasemap R0100010 JX5i3i Ba20g10a2f1g2WSne X2 Y14 K P gt plotps sample1d I1 ltlt END nawk print 1 sqrt1 gt inputdat 0 ID: 480674
Download Presentation The PPT/PDF document "Generic Mapping Tools (GMT)" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Slide1
Generic Mapping Tools (GMT)
MapsSlide2
#!/bin/bash
psbasemap
-R0/100/0/10 -JX5i/3i -Ba20g10/a2f1g2WSne \ -X2 -Y14 -K -P > plot.pssample1d -I1 << END |\nawk '{print $1, sqrt($1)}' > input.dat0100ENDpsxy input.dat -R -JX -W5t15_15:0 -K -O >> plot.pssample1d input.dat -I10 |\psxy -R -JX -St0.5 -G255/0/0 -W5/0/255/0 -O >> plot.psgs plot.ps
Final script for Example 1Slide3
Correction on the –W flag
psxy
input.dat –R –JX -W5t15_15:0 -K –O >> plot.psThe –W command can be quite complicated and is specified in order –W[width],[color],[texture]
In this example, the first 5 represents line width,
the
t
signifies a texture follows, the 15_15 specifies the dash and space width,
and
the :0 specifies that a dash is used.
Alternately,
-W5,0,15_15:-
does the same thing.Slide4
Common command options on first, and possibly subsequent,
calls
Need on all calls-R Define region for plot – will need on first call and at least “–R” on subsequent-J Define
projection
for plot – will need this on all calls if need to define regionSlide5
Common command options on first, and possibly subsequent,
calls
(Generally) Need on first call only-B Borders -- annotation, frame, grid. Only need on first (or a single) call.-P Switch
between landscape and portrait
modes
-X
Shift
X
axis
-Y
Shift
Y axisSlide6
Common command options on first, and possibly subsequent,
calls.
Need when needed.-K Don’t close PostScript (showpage), use when more will follow need
on all but last GMT
call
-O
Don’t
initialize PostScript, use when appending to pre-existing file
need
on all but first GMT call
use
both
–K
and
–O
when putting a large number of GMT call outputs togetherSlide7
Common command options on first, and possibly subsequent,
calls.
Need when needed.-V Verbose (prints out stuff to standard error for user).-H Header
records (tells GMT to skip first H lines of
ascii
input file)Slide8
GMT Defaults
There are about 100 parameters which can be adjusted individually to modify the appearance of plots or affect the manipulation of data. Each as a default value.
GMT defaults are kept in a file called ~/.gmtdefaults4. There are tons of them and you can find out what they are and what the mean reading the man page for gmtdefaults.Slide9
When a program is run, it initializes all parameters to the GMT defaults, then tries to open the file .gmtdefaults4 in the
current directory
. If not found, it will look for that file in a sub-directory ~/.gmt, 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. To view your current gmtdefault settings%gmtdefaults –LTo view the list of options for each default parameter %man gmtdefaults Slide10
Plotting Defaults
example of start of .gmtdefaults4
# GMT-SYSTEM 4.2.1 Defaults file#-------- Plot Media Parameters --PAGE_COLOR = 255/255/255PAGE_ORIENTATION = landscapePAPER_MEDIA = letter#--- Basemap Annotation Parameters --ANNOT_MIN_ANGLE = 20ANNOT_MIN_SPACING = 0ANNOT_FONT_PRIMARY = HelveticaANNOT_FONT_SIZE = 14pANNOT_OFFSET_PRIMARY = 0.075iSlide11
Changing the defaults
You can edit your local copy of .gmtdefaults4 using
nedit or vim You can explicitly reset a default within a script using the command gmtset#!/bin/shgmtset PAPER_MEDIA letter MEASURE_UNIT cmgmtset OUTPUT_DEGREE_FORMAT +D gmtset PLOT_DEGREE_FORMAT +D Slide12
How about
making pretty
MAPS?Slide13
Map projections available in GMTSlide14
List of “standard” command line options
.
The –J option sets the “projection”One has to look at the man page for each one as “different things vary”Slide15
pscoast
-R-90/-70/0/20 -JM6i -P -B5g5 -G180/120/60 > map1.ps
“All” gmt programs plot “maps” through the projection command line option or switch (even the x-y plot).Slide16
pscoast
-R-90/-70/0/20 -JM6i -P -B5g5 -G180/120/60 > map1.ps
All projections give you two selections for specifying the scale (note GMT takes the mapmakers attitude that a map has to have a predetermined/known scale – nicely filling the page does not cut it.)Slide17
pscoast
-R-90/-70/0/20 -JM6i -P -B5g5 -G180/120/60 > map1.ps
-Jmparameters (Mercator [C]). Specify one of:
-
Jm
scale
or
-
JM
width
Give scale along
equator
(
1:xxxx or UNIT/degree).
Slide18
pscoast
-R-90/-70/0/20 -JM6i -P -B5g5 -G180/120/60 > map1.ps
-Jmlon0/lat0/scale
or
-JM
lon0
/
lat0
/
width
Give
central meridian, standard latitude and scale along
parallel
(
1:xxxx or UNIT/degree, UNIT = number inches or
cms
). Slide19
Mercator Projection:
One way to address plotting sphere on a plane (which is whole ‘
nother subject)Conformal (maintains shapes)Cylindrical projectionSlide20
pscoast
-R-130/-70/24/52 -JB-100/35/33/45/6i -B10g5:."
Conic\ Projection": -N1/2p -N2/0.25p -A500 -G200 -W0.25p -P >! map.ps
Region is “rectangle” on the spherical earth
.
-N
for boundaries
(international, US/Canadian/
Mexican
state boundaries “built in”)
, rivers
.Slide21
pscoast
-R-130/-70/24/52 -JB-100/35/33/45/6i -B10g5:."
Conic\ Projection": -N1/2p -N2/0.25p -A500 -G200 -W0.25p -P >! map.ps
-
A
to get rid of small water/island features
Albers
projection
(
b
/
B
) – need to know something (center and standard parallels).Slide22
pscoast
-R-130/-70/24/52 -JB-100/35/33/45/6i -B10g5:."
Conic\ Projection”: -N1/2p -N2/0.25p -A500 -G200 -W0.25p -P >! map.ps -Jb
lon0
/
lat0
/
lat1
/
lat2
/
scale
or
-JB
lon0
/
lat0
/
lat1
/
lat2
/
width
Albers - Give
projection center, two standard parallels, and scale
(1:xxxx or UNIT/degree). Slide23
Albers
Also conformal (maintains/conserves shape)
Conical projectionSlide24
pscoast
-R0/360/-90/90 -JG280/30/6i -Bg30/g15 -Dc -A5000 \
-G255/255/255 -S150/50/150 -P >! map.ps azimuthal Orthographic
projection
(
g
/
G
) mimics
looking at earth from infinite
distance.Slide25
pscoast
-R0/360/-90/90 -JG280/30/6i -Bg30/g15 -Dc -A5000 \
-G255/255/255 -S150/50/150 -P >! map.ps New option
-Dc
Controls resolution of coastline
f
full
h
high
l
low
c
crude
Helps manage file sizes.Slide26
Example 2
Simple
location mapwith focalmechanismsSlide27
gmtset
: To change individual GMT default parameters
pscoast: Plot coastlines, filled continents, rivers, and political borders psxy: Plot symbols, polygons, and lines in 2-Dpstext: Plot text strings psmeca: Plot focal mechanisms on mapsSlide28
#!/bin/sh
###Create a
basemap for the Sumatra-Andaman Island region#Set global variablesPlots=../BasicInfoCPT=/usr/local/GMT4.2.0/share/cptLATMIN=-6.2LATMAX=16LONMIN=90LONMAX=106OUTFILE=$0.ps #What does the $0 signify on this line?#Explicitly set some GMT default valuesgmtset MEASURE_UNIT cm ANOT_FONT_SIZE 12 LABEL_FONT_SIZE 12
gmtset
BASEMAP_TYPE fancySlide29
#Map
# Create the basic
basemap using GMT coastline datapscoast -JM15 -R$LONMIN/$LONMAX/$LATMIN/$LATMAX -Ba2f1WNes -Y5 -X3 -K -P \ -Dh -A100 -N1 -W1 -G155 > $OUTFILE -D Selects the resolution of the data set ((f)ull, (h)igh, (i)ntermediate, (l)ow
, (
c)rude
)
-A Features with an area smaller than
min_area
in km^2
-W and –G Control the fill color and line color of the land regions
-N Plots political boundaries. 1 is for national boundaries, 2 is for states in the AmericasSlide30
#Plot basic tectonic information, including subduction zones and volcanoes
psxy
$Plots/trench.right.gmt -R -JM -M -W5/0 -Sf0.3i/0.08irt -G0 -O -K >> $OUTFILEpsxy $Plots/trench.left.gmt -R -JM -M -W5/0 -Sf0.3i/0.08ilt -G0 -O -K >> $OUTFILEpsxy $Plots/trench.other.gmt -R -JM -M -Wt30_30:1p5/0 -O -K >>
$OUTFILE
psxy
$Plots
/
transform.gmt
-R -JM
-M
-Wt10_10:1p5/0 -O -K >>
$OUTFILE
psxy
$Plots
/
ridge.gmt
-R -JM
-M
-W1/0 -O -K >>
$OUTFILE
psxy
$Plots
/
volcanoes.simkin+siebert.gmt
-R -JM -St.35 -W1/0 -G0 -O -K >>
$OUTFILE
-M Multiple segment file. Segments are separated by a record whose first character is flag. [Default is '>']. Example of
trench.right.gmt
>
92.278 6.948
91.903 7.663
>Slide31
#Information for the 2004 Sumatra event
awk '{print $7, $6, $9/30}'
$Plots/122604_032705.pde |\ #prints lat, lon,mag/30 psxy -R -JM -Sc -W1/0 -G255/255/100 -O -K >> $OUTFILE When no value is specified on the –S flag (ie 0.20), then the size of the circle (denoted by the c) is controlled by the third column of info sent to the psxy command. So we are scaling the circle size by what earthquake parameter in this example? Example of 122604_032705.pde from http://neic.usgs.gov/neis/epic/epic_rect.html PDE-W 2004 12 26 005853.45 3.30 95.98 30 9.00
MwHRV
9C M STS...M
PDE-W 2004 12 26 011710.33 4.94 94.27 30 5.50
mb
GS .. . ....... Slide32
psxy
-R -JM -Sa.89 -W1/0 -G255/255/0 -: -O -K << END >> $OUTFILE3.3 95.98END -Sa creates a star -: indicates that the data is in y-x format rather than x-y format. This is useful because we speak of earthquakes in terms of latitude (y) and longitude (x) and therefore tend to write geographic coordinates in this order as well.Slide33
#Global CMT solutions for events above Mw ~6.6
psmeca
-R -JM -Sm0.5 -G255/255/100 -T -O -K << END >> $OUTFILE#lon lat depth mrr mtt mpp mrt mrp mtp iexp name94.26 3.09 29 1.04 -0.43 -0.61 2.98 -2.40 0.43 29 X Y 122604A
92.79 6.61 14 5.26 -0.84 -4.41 3.95 -2.91 2.10 26 X Y 122604B
92.45 8.58 12 0.94 -0.12 -0.81 0.02 -0.27 0.32 26 X Y 122604C
92.20 4.99 12 -0.02 -0.64 0.66 -0.11 -0.16 -0.80 26 X Y 010105A
95.38 2.84 12 0.26 -0.17 -0.09 0.90 -0.91 0.16 26 X Y 022605A
END
psmeca creates lower hemisphere projections, aka
beachballs
. It is a rather complicated little command.
-S Selects the meaning of the columns in the data file.
-Sa Focal mechanisms in Aki and Richard convention.
-Sc Focal mechanisms in Harvard CMT convention.
-
Sm
Seismic moment tensor (Harvard CMT, with zero trace).
-
Sd
plot only the double couple part of moment tensor.
-
Sz
plot anisotropic part of moment tensor (zero trace).
http://
www.globalcmt.org/CMTsearch.html
-Sp Focal mechanisms given with partial data on both planes.
-
Sx
Principal axis.
-G Controls fill of the compressive quadrant of the moment tensor
-T Plots the nodal planesSlide34
pstext
-R -JM -: -O << END >> $OUTFILE-3.5 90.5 14 0 0 0 Yellow: 2004 Mw 9.2 series-5.75 90.5 14 0 0 0 CMT solutions for Mw >= 6.6END format is x, y, size, angle, font#, justify, text %pstext –L will tell you the font # and font name you can use [T|M|B][L|C|R] to specify justification(top/middle/bottom/left/center/right)gs $OUTFILE #gs indicates ghostscript
if you are using vim to edit this script
:!
script.gmt
will run the script, open the resulting file using
gs
, and you can see what you have done. If you don’t like it, type quit, and continue editing the file.Slide35
nawk
commands as input to GMT
nawk '{print $9, $8, $11}' EBH.HDFYou can put this into GMT several waysIf this is the only file you want to plot – this would worknawk '{print $9, $8, $11}' EBH.HDF |
pxsy
…Slide36
If you had a number of files that needed conversion you could do it this way (only need one
psxy
call)psxy … << END ……`nawk '{print $9, $8, $11}' EBH.HDF`
…
END
Converting each file on the fly.Slide37
If you want to do the same thing to a list of files
filelist
=“$SAMDATA/eq-rupt-1995.dat $DEM/eq-rupt-1960.dat”for FILE in $filelistdopsxy -R -$PROJ$SCALE -M$ -: $CONTINUE -W$LINETHICK/$PURPLE $FILE \$VBSE >> $OUTPUTFILE
done
Other ways to make list
(notice the different kinds of quotes: “,’ and `)
filelist
=`
ls
-1 $ROOT/dem/topocontours/andes_3000_*`
contourlist
=‘1 2 3 4
’Slide38
Some
other
nawk tricks – doing math and passing variables to nawk (quote heaven)SCALE=`echo $STNDTMLON | \nawk ‘{print ($1>=0?$1:360+$1)”/”’${jTRESCALE}_1’*’$FACTOR’}’ `