#!/bin/bash # Written by Collin Ferebee & Scott T. Marshall 3/25/2009 # # This script plots the a map of Northern California # with a shaded DEM from the USGS NED dataset. # On top of the map is plotted earthquakes from the # Waldhauser & Schaff (2008) relocated seismicity catalog. # Earthquake data courtesy of the SCEC datacenter; processed using getEQloc.pl # # #---------------------------------------------------------------------------------------------# #--------------------------------- DEFINE SCRIPT VARIABLES ---------------------------------# #---------------------------------------------------------------------------------------------# #Should I plot the DEM? 1 for yes 0 for no. plotDEM=1 #Name the output file map_file=NorCal_3arc.eps #Define the W/E/S/N edges in degrees area=-125.5/-117.0/35.3/42.5 #The map width in inches scale=39.5i #Define the lighting angle azi=050 #How much should I move the plot over to center it on the page? xShift=1.5i #If the .eps file exists, remove it if [ -f $map_file ] then rm $map_file fi #-----------------------------------------------------------------------------------------# #--------------------------------- GMT SETUP VARIABLES ---------------------------------# #-----------------------------------------------------------------------------------------# #set the GMT variable for eps paper size gmtset PAPER_MEDIA Custom_42ix44.5i #set the GMT variable for dots per inch (resolution) gmtset DOTS_PR_INCH 300 #set the GMT variable for font size gmtset ANNOT_FONT_SIZE_PRIMARY 20 #set the GMT variable for map scale height gmtset MAP_SCALE_HEIGHT 0.25i #----------------------------------------------------------------------------------------# #--------------------------------- PROCESS GRID FILES ---------------------------------# #----------------------------------------------------------------------------------------# #convert from the binary .bil file to a GMT .grd file. GMT will only plot .grd files #Processing piece 1 #xyz2grd bil/NED_25060208.bil -G1arc/1.grd -R-126.0/-123.75/39.0/42.5 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/1.grd #Processing piece 2 #xyz2grd bil/NED_72510048.bil -G1arc/2.grd -R-123.75/-121.5/39.0/42.5 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/2.grd #Processing piece 3 #xyz2grd bil/NED_81525458.bil -G1arc/3.grd -R-121.5/-119.25/39.0/42.5 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/3.grd #Processing piece 4 #xyz2grd bil/NED_90029135.bil -G1arc/4.grd -R-119.25/-117.0/39.0/42.5 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/4.grd #Processing piece 6 #xyz2grd bil/NED_25152779.bil -G1arc/6.grd -R-123.75/-121.5/35.0/39.0 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/6.grd #Processing piece 7 #xyz2grd bil/NED_44713122.bil -G1arc/7.grd -R-121.5/-119.25/35.0/39.0 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/7.grd #Processing piece 8 #xyz2grd bil/NED_88882090.bil -G1arc/8.grd -R-119.25/-117.0/35.0/39.0 -I1c -ZTLh -V -F #give info for the newly created grid #grdinfo 1arc/8.grd #resample .grd files down to 3 arc seconds #grdsample 1arc/1.grd -G3arc/1.grd -I3c #grdsample 1arc/2.grd -G3arc/2.grd -I3c #grdsample 1arc/3.grd -G3arc/3.grd -I3c #grdsample 1arc/4.grd -G3arc/4.grd -I3c #grdsample 1arc/6.grd -G3arc/6.grd -I3c #grdsample 1arc/7.grd -G3arc/7.grd -I3c #grdsample 1arc/8.grd -G3arc/8.grd -I3c #make intensity files for to figure out the range of sigma values for all grd files #grdgradient 3arc/1.grd -A$azi -G3arc/1_int.grd -Nt -V #grdgradient 3arc/2.grd -A$azi -G3arc/2_int.grd -Nt -V #grdgradient 3arc/3.grd -A$azi -G3arc/3_int.grd -Nt -V #grdgradient 3arc/4.grd -A$azi -G3arc/4_int.grd -Nt -V #grdgradient 3arc/6.grd -A$azi -G3arc/6_int.grd -Nt -V #grdgradient 3arc/7.grd -A$azi -G3arc/7_int.grd -Nt -V #grdgradient 3arc/8.grd -A$azi -G3arc/8_int.grd -Nt -V #Now that I know the sigmas from these grids, apply the avg sigma value using the -Nt[amp/sigma/offset] option. #This should give uniform shading for all grid tiles. #grdgradient 3arc/1.grd -A$azi -G3arc/1_int.grd -Nt1/14000/0 #grdgradient 3arc/2.grd -A$azi -G3arc/2_int.grd -Nt1/14000/0 #grdgradient 3arc/3.grd -A$azi -G3arc/3_int.grd -Nt1/14000/0 #grdgradient 3arc/4.grd -A$azi -G3arc/4_int.grd -Nt1/14000/0 #grdgradient 3arc/6.grd -A$azi -G3arc/6_int.grd -Nt1/14000/0 #grdgradient 3arc/7.grd -A$azi -G3arc/7_int.grd -Nt1/14000/0 #grdgradient 3arc/8.grd -A$azi -G3arc/8_int.grd -Nt1/14000/0 #--------------------------------------------------------------------------------------------# #--------------------------------- START PLOTTING THE MAP ---------------------------------# #--------------------------------------------------------------------------------------------# if [ $plotDEM -eq 1 ] then #plot topography and print a progress report echo printing 3arc/1.grd... grdimage 3arc/1.grd -X$xShift -R$area -JM$scale -CCaliColorMap.cpt -I3arc/1_int.grd -P -K >> $map_file echo printing 3arc/2.grd... grdimage 3arc/2.grd -R -JM -CCaliColorMap.cpt -I3arc/2_int.grd -O -K >> $map_file echo printing 3arc/3.grd... grdimage 3arc/3.grd -R -JM -CCaliColorMap.cpt -I3arc/3_int.grd -O -K >> $map_file echo printing 3arc/4.grd... grdimage 3arc/4.grd -R -JM -CCaliColorMap.cpt -I3arc/4_int.grd -O -K >> $map_file echo printing 3arc/6.grd... grdimage 3arc/6.grd -R -JM -CCaliColorMap.cpt -I3arc/6_int.grd -O -K >> $map_file echo printing 3arc/7.grd... grdimage 3arc/7.grd -R -JM -CCaliColorMap.cpt -I3arc/7_int.grd -O -K >> $map_file echo printing 3arc/8.grd... grdimage 3arc/8.grd -R -JM -CCaliColorMap.cpt -I3arc/8_int.grd -O -K >> $map_file #plot basemap mark every degree psbasemap -R -JM -Bf60ma60m/f60ma60m:horizontal:WSne -O -K --ANNOT_FONT_SIZE_PRIMARY=28 --FRAME_WIDTH=0.20i --ANNOT_OFFSET_PRIMARY=0.25i >> $map_file else #plot basemap mark every degree psbasemap -X$xShift -R$area -JM$scale -Bf60ma60m/f60ma60m:horizontal:WSne -P -K --ANNOT_FONT_SIZE_PRIMARY=28 --FRAME_WIDTH=0.20i --ANNOT_OFFSET_PRIMARY=0.25i >> $map_file fi #Plot the coast and a fancy scale. Make the water blue and the land white pscoast -R -JM -Na/1 -W4 -S102/153/204 -Df -O -K >> $map_file #plot seismicity as red circles echo printing earthquake locations... psxy NorCalEQ.mag -R -JM -Sc0.20c -W1 -H1 -G190/30/45 -P -O -K >> $map_file #plot fault traces #psxy ../../SoCal/Lin_et_al/fault.lin -R -JM -H0 -m -W8/0/0/0 -P -O -K >> $map_file #plot beachballs psmeca BerkeleyMomentTensors.txt -R -JM -Sa1/20/100 -W1p -H1 -P -O -K >> $map_file #Plot significant historical earthquakes as stars psxy historicalEQ.txt -R -JM -H1 -W4 -G210/172/043 -Sa0.5i -O -K >> $map_file pstext historicalEQ.txt -R -JM -H1 -D0/0.35i -G210/172/043 -S4 -O -K >> $map_file #plot a legend pslegend -R -JM -Dx0.5i/0.5i/13.5i/8.5i/BL -F -G255/255/255 -P -O -K <> $map_file G 0.25i H 36 1 Northern California Earthquakes > M@-W@-0 (1984-2003) G 0.15i H 28 1 301,888 Total Events in 20 Years; 15,094 Events/Average Year G 0.25i H 24 1 Compiled by Collin L. Ferebee and Scott T. Marshall G 0.15i H 24 1 Department of Geology, Appalachian State University G 0.35i S 1.0i c 0.2c 190/30/45 0.5 1.5i @:20:Earthquake Epicenters from Waldhauser and Schaff (2008)@:: G 0.20i S 1.0i c 0.0001i 255/255/255 0.0 1.5i @:20:Focal Mechanisms for > M@-W@- 5.0 from the Berkeley @:: S 1.0i c 0.0001i 255/255/255 0.0 1.5i @:20:Seismological Lab http://www.ncedc.org/ncedc/mt.html@:: G 0.2i S 1.0i a 0.5i 210/172/43 4.0 1.5i @:20:Large Historical Epicenters from Stover and Coffman (1993)@:: G 0.50i #S 1.0i y 0.4i 0/0/0 8.0 1.5i @:20:Fault Traces from Jennings (1992)@:: #G 0.25i M 242.5 34.5 100k+u f G 0.25i M 242.5 34.5 100m+u f G 0.5i B CaliColorMap.cpt 0.5i 0.5i -B:"Elevation Above Sea Level (meters)": --LABEL_FONT_SIZE=20 --LABEL_FONT=0 END #plot the north arrow and directional rose pscoast -R -JM -W0 -Df -Tfx12.25i/5.25i/2i/2 -O -K --HEADER_FONT_SIZE=26 --HEADER_OFFSET=0.2c >> $map_file #plot a beachball for the legend psmeca -R -JM -Sa1/20/100 -W1p -O -K <> $map_file -125.165 36.227 24.0 320 90 180 5.0 0 0 END #This is what you need to finish off an .eps file correctly echo showpage >> $map_file #------------------------------------------------------------------------------------------------------------------# #--------------------------------- DISPLAY THE MAP AND CONVERT TO OTHER FORMATS ---------------------------------# #------------------------------------------------------------------------------------------------------------------# #display $map_file in a program that can read .eps files #gv $map_file & #acrobat $map_file & #illustrator $map_file & #make a pdf, jpg, and transparent png from the eps files echo converting to pdf ps2raster $map_file -Tf #Now open the pdf with acrobat acrobat NorCal_3arc.pdf & echo converting to jpg ps2raster $map_file -Tj echo converting to png ps2raster $map_file -TG #print a done message and ring the system bell echo -e "Done! \a" exit #--------------------------------------------------------# #--------------------------------------------------------#