Intro to shells, SAC, and GMT

1. Intro to shells, SAC, and GMT Rob Porritt UC Berkeley PhD candidate rob@seismo.berkeley.edu

2. Outline • Shells • History • Scripting • SAC • Getting SAC • Interactive • Within a script • GMT • Example script

3. Working on a terminal

4. 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

5. Using shells • Anytime you use a terminal your are explicitly using a shell

6. Scripts • But you can also put a bunch of commands into a script to run repeatedly!

7. 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

8. 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)

9. Getting SAC • http://www.iris.edu/software/sac/sac.request.htm

10. Using SAC • Command line based Seismic Analysis Code • Set of tools to process data. Frequency filtering, decimation, fft, correlation, integration, differentiation, etc… Ground Velocity (m/s) Raw data (counts)

11. Doing it automatically • Make a script!

12. 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

13. 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)

14. SRTM • 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 From http://spaceplace.jpl.nasa.gov/en/kids/srtm_make2.shtml

15. 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

16. Anatomy…. Make a variable for the output plot name 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 case the projection is Mercator and the scale is such that the x-axis of the map is 5.5 inches long Draws the basemap psbasemap and the four variables misc variable contains other information: -V be verbose, -P portrait mode, -X2i shift by 2 inches in x Return to the terminal and type chmod +x plot_socorro_topo.gmt Then run the command with ./plot_socorro_topo.gmt View the plot with a viewer such as ggv (ie ggv socorro.ps)

17. 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

18. 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 grdimagewill draw a color image of the topography in socorro.grd using the color palette (-C) socorro.cpt -K tells the interpreter there is more commands to be appended >> indicates postscript code is being appended to the file created by the previous command -O tells the interpreter to overlay this information on a previous map Run the c-shell

19. 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

20. Changing the Color Palette In this case I have used the haxby color palette – experiment with other palettes

21. The grd2cpt command has been commented out – this does not need to be run any more Add Illumination grdimage has had the paramter –I added to it. This indicates that the resultant image will be illuminated with the file socorro.grad 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

22. The Map…. There is now a lot more information in the map – you can start to see structures and river channels By changing the –Aoption in grdgradient, try recreating the image with different illumination angles. Does one angle help us to visualize structures better than another?

23. Add Contours The grdgrdient command has been commented out – this does not need to be run any more Note that we are now using –O and –K. This indicates that this is neither the top of bottom layer in the postscript file grdcontour creates a layer in the postcript file that will display contours every 100 m (-C100)

24. 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

25. Add A Color Scale The bar is labeled every 500 m and annotated with “meters” 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 psscale draws a color bar using the color palette socorro2.cpt

26. 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?

27. 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 C-Shell . 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) A star (-Sa) with a 0.4 inch diameter is drawn This is the longitude and latitude of the library

28. 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),

29. Add Some Text to the Map The text will be black 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 C-Shell 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 This is the longitude and latitude at which the text will be plotted

30. 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….

31. Add an Inset to the Map We are adding more layers after psscale, so we have added a -K 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 We are redefining the area of the map – the inset will include the contiguous states. It will be 2 inches across A black star is plotted at Socorro pscoast is used to plot a low-resolution (-Dl) coastline is plotted. It has a line thickness of 1 point and is black

32. 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.

33. 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

34. With profile