#!/bin/sh  # the first line in any shell script must beingn with this

# TUTORIAL GMT SCRIPT (Not finished. There is still a lot missing.)
# Geoffrey Ely 8/7/97
# updated Joe Mount 2006-01-17 
# This script creates a earthquake and fault map of California and Nevada.
# It is the same as the script for the "cal" map from the catalog with
# exaplanations added for each line of code.
# It may be helpful to bring up the images of the "cal" map in a separate 
# browser window while you go through this script.

# Set the variable O to the name of the PostScript output file.
# Each GMT commands creates a section of Postscript that will be added to this file.
O="cal.ps"

# Set the variable D to the name of the directory where the data files reside.
D="/net/quake/gis/geoff/gmt/mapdata"

# Reset the GMT defaults to the original settings.
# This avoids using unwanted settings created by other GMT scripts.
# GMT deafaults are changed with the gmtset command in the following lines.
# Type man gmtdefaults for more info.
rm .gmtdefaults .gmtcommands

# Use inches.  Default is centimeters. 
# Set paper width.  This should always be the shortest dimension
# regardless of the orientaion of your map.
# NOTE: PAPER_WIDTH is no longer used in GMT_3.4.4 and later versions!
gmtset	MEASURE_UNIT inch		PAPER_WIDTH 8.5 (obsolete)	

# Use plain border.  Default is fancy (alternating black/white frame).
# Label longitudes are unsigned 0 to  180,  latitudes unsigned 0 to 90.
gmtset	BASEMAP_TYPE PLAIN		DEGREE_FORMAT 2

# Set font size for axis lables and contour labels (not used in this script).
# Set thickness of the basemap lines.
gmtset	LABEL_FONT_SIZE 6		FRAME_PEN 1

# Set font size for lon/lat anotations.
# Set offset between map fram and lon/lat anotations (in inches).
gmtset	ANOT_FONT_SIZE 6		ANOT_OFFSET 0.075

# Set thickness of lon/lat tickmarks.
# Set length of tickmarks (in inches).  Negative number points tickmarks inward.
gmtset	TICK_PEN 1			TICK_LENGTH -0.05

# Set the minimum the angle between the map boundary and the annotation baseline for
# which annotations are plotted.
gmtset  ANOT_MIN_ANGLE 45

# Make a Color Pallete Table (CPT).
# The cpt will be used later to color the a topography grid by elevation.
# This cpt is assignes a grey to lower elevations that get subtley lighter
# with elevation. 
# See the GMT Technical Reference for more information on CPTs.
cat					<< END			> cpt
B	200	200	200
1	200	200	200	20	210	210	210
20	210	210	210	75	220	220	220
F	220	220	220
END

# Extract the California/ Nevada region from the US topography grid and
# make a sepatate file.
# The 60-second data will be used for speed while composing the map.
# The 30-second data will be used for the final version.
# These lines get commented out after the running them once.
#grdcut -V -R-126/-113/32/43		$D/topo/usa30.grd	-G$D/topo/cal30.grd
#grdcut -V -R-126/-113/32/43		$D/topo/usa60.grd	-G$D/topo/cal60.grd
 
# Make a gradient file from the topography.  This files is used to give
# the shaded look to the topography.
#grdgradient -V -Nt -A270		$D/topo/cal30.grd	-G$D/topo/cal30i.grd
#grdgradient -V -Nt -A270		$D/topo/cal60.grd	-G$D/topo/cal60i.grd

# Preform different mathmatical functions on the gradient file
# to change the look of the of the shading.
#grdmath -V $D/topo/cal60i.grd SIGN $D/topo/cal60i.grd 25000 DIV ABS SQRT MUL = $D/topo/cal60i2.grd
#grdmath -V $D/topo/cal60i.grd .2 ADD 1.25 MUL = $D/topo/cal60i2.grd
#grdmath -V .25 $D/topo/cal60i.grd ABS .5 SUB 2 POW SUB $D/topo/cal60i.grd ADD = $D/topo/cal60i2.grd

# Set the flags that will be used in all of the GMT commands for this section.
# The two definitions for F have slightly different ranges and projections.
# 
# -R Selects the range.
# -J Selects the map projection. (M=mecator, L=lambert)
# -V Selects verbose mode, progress reports will be
#    displayed while each command is running.
# -K means more PostScript will be appended later.
#    This flag must appear in every GMT command except for the last one.
# -P Portrait mode.  Default is landscape.
F="-R-125:30/-113:30/32:00/42:45 -JM7.5 -V -K -P"
F="-R-125/32:01/-113:15/42:30r -JL-119.5/38/35/40/7.2 -V -K -P"
if [ ! -f $O.tmp ]; then
grdimage $F -O -Ccpt -I$D/topo/cal60i2.grd $D/topo/cal60.grd	> $O.tmp; fi
psclip	$F -X.65 -Y1.65 -M $D/ca_nv.xy | sed ' 9s/^/.9999 .9999 scale/' > $O
echo "2 setlinejoin"						>> $O
cat					$O.tmp			>> $O
psclip	-S -O -K						>> $O
awk '$3=="c"' $D/ca_nv.xy | psxy $F -O -W1 -M			>> $O
psxy	$F -O -L -W1 -M -G210		$D/ca_isl.xy		>> $O
awk '$3=="b"' $D/ca_nv.xy | psxy $F -O -W1t5_5_25_5_5_5:0	>> $O
psxy	$F -O -W1t5_5_25_5_5_5:0	$D/ca_nv_border.xy	>> $O
psbasemap $F -O -Lfx5.27/-.3/37/200				>> $O	
gmtset	ANOT_OFFSET -0.15
psbasemap $F -O -Bf1a2NSew					>> $O
gmtset	ANOT_OFFSET -0.20
psbasemap $F -O -Bf1a2nsEW | sed '/2156 284 M V -86.7/d'	>> $O
psxy	$F -O -M -W4/0/105/255		$D/cdmg/q		>> $O
psxy	$F -O -M -W4/205/0/205		$D/cdmg/lq		>> $O
psxy	$F -O -M -W4/255/185/0		$D/cdmg/holo		>> $O
psxy	$F -O -M -W4/255/0/0		$D/cdmg/hist		>> $O
awk '$4<=5.5' $D/eqs/cnss19* | psxy $F -O -W2/255/0/0 -Sc	>> $O
awk '$4> 5.5' $D/eqs/cnss19* | psxy $F -O -W5/255/255/0 -Sc	>> $O

F="-R0/11/0/11 -Jx1 -V -O -N"
pstext	$F -K -Y-1			<< END			>> $O	# title
3.60 9.40 16 0 29 2 EARTHQUAKES AND FAULTS IN CALIFORNIA AND NEVADA
1.00 1.25  7 0 30 1 Geoffrey P. Ely
0.25 0.70  8 0  0 1 Age of Faulting
0.60 0.55  8 0  0 5 Historical
0.60 0.40  8 0  0 5 Holocene
0.60 0.25  8 0  0 5 Late Quaternary
0.60 0.10  8 0  0 5 Quaternary
2.00 0.05  8 0  0 2 4
2.25 0.05  8 0  0 2 5
2.50 0.05  8 0  0 2 6
2.75 0.05  8 0  0 2 7
3.00 0.05  8 0  0 2 8
2.50 0.70  8 0  0 2 Earthquake Magnitude
3.70 0.29  7 0  0 1 Fault map provided by the California Division of Mines and Geology.
3.70 0.17  7 0  0 1 Earthquake locations provided by the Council of the National Seismic System.
3.70 0.05  7 0  0 1 Topographic data provided by the United States Geological Survey.
END
psxy	$F -K -W4/255/0/0		<< END			>> $O
0.25 0.55
0.50 0.55
END
psxy	$F -K -W4/255/185/0		<< END			>> $O
0.25 0.40
0.50 0.40
END
psxy	$F -K -W4/205/0/205		<< END			>> $O
0.25 0.25
0.50 0.25
END
psxy	$F -K -W4/0/105/255		<< END			>> $O
0.25 0.10
0.50 0.10
END
psxy	$F -K -W2 -Sc			<< END			>> $O
2.00 0.212 0.024
2.25 0.224 0.048
END
psxy	$F -K -W5 -Sc			<< END			>> $O
2.50 0.248 0.097
2.75 0.297 0.194
3.00 0.387 0.388
END

# To plot a custom symbol, first a Postcript definition file must be made.
# Here the definition file "ics.def" is catted into the plot file.
# It is then plotted with the psxy command.
# The psxy output is piped to a sed command which tricks psxy into plotting the
# definition "ICSBox" instead of "C4", (C4 is the definition for a circle). 
# This is the ICS logo
cat					$D/ics.def		>> $O
psxy	$F -K -Sc.7 << END	| sed 's/C4/ICSbox/'			>> $O
0.90 1.80
END

psxy	$F -X-.65 -Y-.65 -W3 -L	<< END			>> $O
0.55  0.55
0.55 10.45
7.95 10.45
7.95  0.55
END