Transcript Intro to shells, SAC, and GMT

Intro to shells, SAC, and
GMT
Rob Porritt
UC Berkeley PhD candidate
[email protected]
Outline
 Shells
 History
 Scripting
 SAC
 Getting SAC
 Interactive
 Within a script
 GMT
 Example script
Working on a terminal
History of shells
 1977 by Stephen Bourne at AT&T Bell Labs – the
bourne shell (/bin/sh)
 1978 team at Berkeley – the C-shell (/bin/csh)
 Bourne shell updated to the bourne again shell
(/bin/bash)
 C-shell updated to the TCSH shell (/bin/tcsh)
 Many more, but these are the most common
Using shells
 Anytime you use a terminal your are explicitly using a
shell
Scripts
 But you can also put a bunch of commands into a script
to run repeatedly!
Making a script
 First line:
#!/bin/tcsh
Tells the default shell to use the interpreter located at
/bin/tcsh
 Write a series of commands and #defines comments
 chmod +x filename
 ./filename
Useful commands
 ssh – login to a remote computer
 set – make a variable which can be returned with $var
 foreach (tcsh)
 screen – run a process not connected to a terminal (ctrl
+ a + d to detach. screen –r to reattach)
Getting SAC
 http://www.iris.edu/software/sac/sac.request.htm
Using SAC
 Command line
based Seismic
Analysis Code
 Set of tools to
process data.
Frequency
Ground Velocity
(m/s)
filtering,
decimation, fft,
correlation,
integration,
Raw datadifferentiation,
(counts)
etc…
Doing it automatically
 Make a script!
GMT
 Paint by numbers
 Start with example scripts and edit for new purposes
 Rest of the time we’ll spend making a script to make a
map of Socorro
 Original tutorial by Andy Goodliffe at the University of
Alabama – editted by Rob
 You can use this script for future plots – just copy,
paste, edit
Where did the data come from?
 Space Shuttle Radar Topography Mission (SRTM)
data
 1 Arc-second (~30 m) grids
 Downloaded from the USGS Seamless server
(http://seamless.usgs.gov/) and converted to a GMT
format (NetCDF) grid
 Re-sampled (using grdfilter) to a 0.0005 degree grid
 Also available from NASA
(ftp://e0srp01u.ecs.nasa.gov)
SRTM
From http://spaceplace.jpl.nasa.gov/en/kids/srtm_make2.shtml
•A transmit antenna illuminates the
terrain with a radar beam which is
scattered by the surface. Both
arrays receive the reflected signal
•The signal coming to one antenna
may have traveled slightly further
than that arriving at the other
•This translates into a phase
difference
•By measuring the phase
difference we can determine the
angle from which the signal came
•When combined with travel time,
we can determine the distance to
that point
The Basic C-Shell
 Open your favorite text editor on the UNIX system
and create a file named plot_socorro_topo.gmt
 File names are arbitrary – Make them reasonable so
you know what they do
 If you’re a glutton for punishment like me you can
use VIM by typing vi plot_socorro_topo.gmt
 Or you can use xemacs by typing xemacs
plot_socorro_topo.gmt &
 The & at the end puts the xemacs command in the
background so that you get your prompt back
Make a variable for the output plot name
Anatomy….
Indicates that the script is a TC-Shell
Sets the format of the map annotation
Sets a variable for the area of the map
Sets a variable for the tick marks on
the map – annotations every 0.2
minutes, tick marks every 10 minutes
Sets a variable for
the map projection
and scale – in this
Draws the
misc variable contains other information:
case the projection
basemap
-V be verbose, -P portrait mode, -X2i
is Mercator and
psbasemap
the scale is such
and the four shift by 2 inches in x
that the x-axis of
variables
the map is 5.5
Return to the terminal and type chmod +x plot_socorro_topo.gmt
inches long
Then run the command with ./plot_socorro_topo.gmt
View the plot with a viewer such as ggv (ie ggv socorro.ps)
The Map…
Hopefully your map looks something like
this (depending on settings in your
.gmtdefaults file)
Edit your C-Shell so that
set boundaries = “-Ba0.2f10m”
Now reads
set boundaries = “-BSWnea0.2f10m”
This will make it so that only the
southern and western axes are labeled
Note: to find other options in
psbasemap, simply type man psbasemap
– this will display the help pages
Add Topography
You have each been given a two files – socorro.grd and socorro.cpt – these are
the SRTM grid files and color palette file respectively
grdimage will
draw a color
image of the
topography in
socorro.grd
using the color
palette (-C)
socorro.cpt
-O tells the
interpreter to
overlay this
information on a
previous map
-K tells the
interpreter
there is more
commands to be
appended
Run the c-shell
>> indicates
postscript code
is being
appended to the
file created by
the previous
command
The Map….
White is the highest
topography, brown is the lowest.
We can all agree that this color
palette is ugly – lets take a
moment to create a new one
using grd2cpt
Type man grd2cpt to see how to
do this…. Or simply type
grd2cpt
Changing the Color Palette
In this case I have used the
haxby color palette –
experiment with other palettes
The grd2cpt
command has
been commented
out – this does
not need to be
run any more
Add Illumination
grdgradient creates an illumination file with simulated
illumination from the north (-A0), dimensionless gradient (-M)
normalized to 1 (-Nt). The output illumination file is
socorro.grad
grdimage has
had the
paramter –I
added to it.
This indicates
that the
resultant
image will be
illuminated
with the file
socorro.grad
The Map….
There is now a lot more
information in the map – you can
start to see structures and
river channels
By changing the –A option in
grdgradient, try recreating the
image with different
illumination angles. Does one
angle help us to visualize
structures better than
another?
The grdgrdient
command has
been commented
out – this does
not need to be
run any more
grdcontour
creates a
layer in the
postcript
file that will
display
contours
every 100 m
(-C100)
Add Contours
Note that we
are now using
–O and –K.
This indicates
that this is
neither the
top of
bottom layer
in the
postscript
file
The Map….
We now have an idea of the size
of the mountains around us – we
know from running the script last
time that the minimum and
maximum contours are 1300 m
and 3200 m respectively
By changing the –C option in
grdcontour, change the contour
interval to 200 m. By using the –A
option, add annotations every
1000 m
Add A Color Scale
The bar is
labeled every
500 m and
annotated with
“meters”
psscale draws
a color bar
using the color
palette
socorro2.cpt
This is the location (x
– 1.25 inches, y = -0.5
inches) and size (4 x
0.25 inches). h
indicates that the bar
is horizontal
The Map….
The uninitiated reader now
everything they need to know
most of the map details.
Why have I not added a regular
scale bar?
Add
the
Location
of
the
Library
psxy plots
points on maps
using a number
of different
methods (see
man psxy).
Typically we
would read
these points
from a file, but
since we are
only plotting
one point, we
will read it
from the CShell . psxy can
also be used to
plot a regular
graph
The color fill is the
star is white (255)
As input, psxy reads
until it reaches END
The line
thickness of
the star is 1
point (1p) and
the line is black
(0)
This is the longitude and
latitude of the library
A star (-Sa) with
a 0.4 inch
diameter is drawn
The Map….
Colors in GMT are specified in
RGB format – thus 255/0/0 is
red, 0/255/0 is green, etc. You
can mix colors to create the color
of your choice. By changing the –G
option in psxy, change the color of
the star (for example, -G255/0/0
would make the star red),
Add Some Text to the Map
ptext plots
text on maps
Typically we
would read
text strings
from a file,
but since we
are only
plotting one
text string we
will read it
from the CShell
This is the longitude and latitude
at which the text will be plotted
The text will be
black
The text will be
centered on the
point (CM)
Font number 1
(helvetica) will be
used
The text will be
plotted horizontally
(0 degrees)
The font size is
14 points
The Map….
The text is plotted on top of the
star. By changing the text
justification (currently CM)
and/or the text position, move the
text so that it is under the star.
You might want to also change the
text color….
Add an Inset to the Map
We are
adding more
layers after
psscale, so we
have added a
-K
We are
redefining
the area of
the map – the
inset will
include the
contiguous
states. It will
be 2 inches
across
pscoast is used to plot a low-resolution (-Dl)
coastline is plotted. It has a line thickness of 1 point
and is black
This creates the
basemap for the
inset. The
background is
white (-G255), axis
are drawn without
label or ticks (Bnesw), and the
origin is offset
horizontally by 0.2
inches (-X0.2) and
vertically by 7.65
(-Y7.65) inches
relative to the
previous basemap
A black star
is plotted at
Socorro
The Map….
What can you do with this map
now? The postscript format image
can be read by a number of other
packages – Adobe Illustrator
being one example….
This map was the end point for
last year’s tutorial. If there’s time
and willingness I can also show
some other tools to extract data
from a grid and use psxy for an xy
data plot.
The full script with lots of
comments is available on the
resource website.
Code to make the profile
Extract xyz file info from the
topo grid
If the latitude is right, spit out
the longitude and elevation
to psxy
With profile