-
Notifications
You must be signed in to change notification settings - Fork 4
/
cascadia_seis_wcomments.com
executable file
·241 lines (169 loc) · 10.1 KB
/
cascadia_seis_wcomments.com
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#!/bin/bash
# This file is for the first class on data visualization with GMT
# Written by Gina Schmalzle
# This code is a bash shell script (http://en.wikipedia.org/wiki/Bash_(Unix_shell)) that calls GMT files.
# awk is used to manipulate data files (http://en.wikipedia.org/wiki/AWK).
# A unix shell is a "command line interpreter" that allows definition of variables,
# and execution of commands that could be done via the command line. A nice shell script
# will make it easier to change variables, and easier to implement several lines of code.
# Shell scripts are commonly used with GMT programs.
# Most if not all of the commands in this document have associated 'man' (manual) pages. To access them type:
# man whatever_your_command_is
# If you cannot access your man pages through your command prompt, an alternative would be to type man command in
# google.
# To make this file executable, you will have to change the mode of the file (ie, read, write and/or execute)
# In your directory you will need to type:
# chmod u+x ./cascadia_seis.com
# The '#' marks to the left indicates a comment. Anything written after them is not read when the file is executed.
# This file will create a map of Cascadia that includes a grid of topography data from ETOPO1 (ETOPO1_Bed_g_gmt4.grd) and
# seismicity data. These data will be applied in "layers", very similar to how GIS packages have layers. The layers may be
# turned off or on by commenting/uncommenting lines.
# The grd file is already in gmt format. Generating and using grids is another class in itself, but here I will introduce
# you to using GMT formated grid files.
# Also included on the map are earthquakes locations color coded by depth from the ANSS catalog for 2000 to 2012 (anss_eq_2000_2012.dat)
#MAKE A MAP!
# Define the names of the input and output files
out=cascadia_seis.ps # This will be the name of your map generated by this file
seis_data=anss_eq_2000_2012.dat # ANSS earthquake catalog
topo=./ETOPO1_Bed_g_gmt4.grd # ETOPO1 topography grid
# Define map characteristics
# Define your area
north=50
south=40
east=-118
west=-132
# Define your map boundary annotation
# Here we define tick marks every 2 degrees and we print the degree on the West and South sides of the plot
# and keep the ticks (but don't label) on the east and north sides
tick='-B2/2WSen'
# Define Map Projection
# Here we define a Mercator Projection of size = 15
proj='-JM15'
#Start with GMT commands with embedded definitions....
# Help with any of these commands can be obtained by looking at the 'man' files. Simply type at the command line: man gmt_command
# If the man files are not properly installed you can also type in man gmt_command (e.g., man psbasemap) in google and it will come up.
#This line sets up the 'basemap' meaning here you will define the region, boundary annotations and projections.
#You can accomplish this also with other commands (including psxy, pscoast, etc...), but it is good many times to start with psbasemap.
psbasemap -R$west/$east/$south/$north $proj $tick -P -Y12 -K > $out
# This is your first line of GMT Code!!! Whoo-hoo! In long hand this line would look like this:
#
# psbasemap -R-132/-118/40/50 -JM15 -P -Y12 -K > cascadia_seis.ps
# What the options mean:
# psbasemap = plots postscript basemaps
# -R -- defines the area of your map (note that we defined north, south, east and west above and they are inserted into the -R option.
# The Projection (-JM) and tick marks (-B) were defined above.
# Note that when you call a defined variable, you must include a '$' before the variable name
# -P Sets the figure to "Portrait" mode. No -P is landscape.
# -Y Orients the figure vertically (-X orients it horizontally).
# -K means that there will be more 'stuff' appended to the postscript file.
# '>' means that the command output, which would normally print to screen will be directed into your new file (cascadia_seis.ps, shown here as $out)
# In addition it means that it believes cascadia_seis.ps is a new file. If it is not, it will erase all existing info in the file and re-write it with
# the new information.
#plot grid
# We would like the topography to be the map background, so it needs to be the first layer. Hence, we get started with a hard part...
# Helpful hint...
#
# use grdinfo your_grd_file.grd
# to find info about your grid file, such as the min and max values
#
# You will need to make some color palettes. These are files that tell what colors certain properties are displayed.
# For example, your ETOPO grid has a latitude, longitude and a elevation, and you want to color code the topography
# by elevation. The following lines will tell you how to do that...
# First, Make a color palatte
# Typing: makecpt
# at the command line will give you information on pre-existing color schemes
# This will make a color pallete of typical, pre-defined topography colors:
makecpt -Crelief -T-8000/8000/500 -Z > topo.cpt
#makecpt = makes GMT color palette tables
#-C tells GMT what pre-defined color palette to use
#-T defines the range and increment
#-Z states that the colors will change continuously (rather than discretely)
#topo.cpt is a new file containing your color pallete information that will be used later.
#This next line is not necessary, but may be used to make the image appear sharper.
#grdgradient helps to illuminate ridges in the topography from a specified angle.
#grdgradient $topo -A135 -Ne0.8 -Gshadow.grd
#grdgradient=Makes illumination shadow
#-A is the angle from which the light is shown
#-N normalizes the shadow according to equations stated in man grdgradient
#-G lists the name of your output grid
# Overlay the grid onto your map
# Here you are adding the grid as a layer to your postscript file
# This command includes a shadow grid file:
# grdimage $topo -R -J -O -K -Ctopo.cpt -Ishadow.grd >> $out
# This command omits the shadow file:
grdimage $topo -R -J -O -K -Ctopo.cpt >> $out
#grdimage = creates an image from a 2D netcdf grid file
#-R = Sets the region. Notice here I don't have to state the min and max values again.
#-J = Sets the projection. Again the type and size don't have to be restated.
#-O = Overlay. The output for this line is being appended to a previous postscript code
# i.e., you are adding another layer
#-K = You will be appending another layer
#-C = You will be using the color pallete topo.cpt
# Now, back to the easy stuff..
# Add coastlines
pscoast -R -J -O -K -W2 -Df -Na -Ia -Lf-130.8/46/10/200+lkm >> $out
#pscoast = adds coastlines
#-W = Sets the line width and color. Default color = black = 0 and does not have to be explicitly stated.
#-Df = What is the resolution of the coasline dataset? f = fine
#-Na = Draws politcal boundaries, a = draw all the boundaries, see man pscoast for more options
#-Ia = Draw Rivers, a = draw all rivers, see man pscoast for more options
#-Lf = Draw a fancy map scale, f = fancy, centered on -130.8, 46 degrees. +200 = length, +lkm = kilometers
# Add seismic locations and color code them by depth
# Make color palette
# Ahh, another color pallete...
# This time, let's make it rainbow colored and call is seis.cpt
makecpt -Crainbow -T0/50/10 -Z > seis.cpt
# Columns 4, 3 and 5 of the data file are the longitude, latitude and depth, respectively. This is the order
# your data need to be in for psxy (see man file)
awk '{print($4,$3,$5)}' $seis_data | psxy -R -J -O -K -W.1 -Sc.1 -Cseis.cpt -H15 >> $out
# psxy = Plot 2D lines, polygons and symbols on a map. Fun fact -- psxyz plots in 3D.
# -W.1 = Draws the black outline of the circles.
# -Sc.1 = Defines the shape and size; c = circle, size = 0.1
# -H = Header. The first 15 lines of the file contain header information and will not be read.
# -C = defines the color palette to be used for the depth. We could also make all the circles one color.
# In this case, remove the -C option and use -G instead. -G defines the color of the circle in either white-black
# or red/green/blue format. Example colors: -G0 (black); -G255 (White); -G255/0/0 (Red)
# GMT has made this a little easier. You could also say -Gblack or -Gred, but there are a limited amount of colors
# you could do that with.
# Add a scale
psscale -D0/3.2/6/1 -B10:Depth:/:km: -Cseis.cpt -O -K >> $out
# pscale = Adds a scale to go with your color palette
# -D = set the position of the scale
# -B = set and annotate the scale tick marks and lables.
# -C = specify your color palette
# Now, let's take a subset of seismic data and project them onto a line....
#First, let's view the transect line
#Plot transect line
psxy center.dat -R -J -O -K -W1 -Sc.3 -G255/0/0 >> $out
psxy center.dat -R -J -O -K -W5/255/0/0 >> $out
# You should know the options by now ;-)
#This ends the map making part of this exercize, now we move onto making a scatter plot from the seismic data.
# PROJECT DATA
# Here we use the GMT code project to take all the data within a certain region and project them onto a line
awk '{print($4,$3,$5)}' $seis_data | project -C-124/47 -A90 -W-.2/.2 -L0/4 -H15 > projection.dat
# project = projects data onto a transect
# Note that the options are different for this command
# -C = defines the center of your transect
# -A = azimuth of transect (CW from N)
# -W = Width of the transect in degrees
# -L = length of transect in degrees
# -H = Header declaration
# projection.dat = new file with the original data and the projected locations
# MAKE SCATTER PLOT
#We want the scatter plot to be on the same page as the map, but just below it, so we need to redefine our
#region, projection and tick marks...
east=-120
west=-124
dmin=0
dmax=50
proj=-JX15/-5
tick=-B1:Longitude:/10:Depth:WSen
awk '{print($6,$3)}' projection.dat | psxy -R$west/$east/$dmin/$dmax $proj $tick -W1 -Sc.2 -G200 -O -K -Y-8 -P >> $out
# Columns 6 and 3 are the projected longitude and the Depth, repectively
# -Y = Shift the new plot down 8 units. You can designate if you want to shift in centimeters (c), inches (i),
# meters (m), or pixels (p). Otherwise it shifts by whatever is in your gmtdefaults.
# Last, but not least, image your map!
# Common postscript viewers: gs, gv, ggv, open, gimp
# What, you don't like postscript files? That's ok, uncomment this line:
# ps2pdf $out
open $out