#!/bin/bash ################################################################################ # # scripts to create the maps in the ftp directory # # please see references.?.txt for original data sources and properly refer # to the sources # # to run this script, you need a GMT installation, need to modify the gmtmedia.d # file to allow the "poster" (36'' x 72'') format, need to download the data in this directory, # and need to change the filenames to the data as listed in this script # # # $Id: world_poster.gmt,v 1.2 2014/04/07 18:20:26 becker Exp becker $ # # thorsten becker, twb@usc.edu # # # ################################################################################ shade=0 # shading on / off topo=1 # actually plot topography addage=1 # seafloor age views=${1-"0 1"} # 0: PAC 1: AFR centered view cmts=${2-"0 1 2"} # 0: engdahl seismicity color coded by depth # 1: gCMT and Endahle, color coded by depth # 2: gCMT and Endahl, < 50km, color coded by style # Location of the GMT binaries: gmtbin=$HOME/progs/src/GMT4.5.7/bin/ verbose="" # # select colormap # #cmp=topo.cpt #cmp=dave_topo.cpt # cmp=topo_topo.cpt cmp=topo_topo.c.cpt #volc_col=128/000/128 volc_col=128/000/128 # volcano color if [ -s $datadir/etopo1/ETOPO1_Ice_g_gmt4.grd ];then topo_file=$datadir/etopo1/ETOPO1_Ice_g_gmt4.grd topo_ifile=$datadir/etopo1/ETOPO1_Ice_g_gmt4.i.grd # illumination else topo_file=$datadir/etopo2/etopo2.grd topo_ifile=$datadir/etopo2/etopo2.i.grd fi #agegrid=$datadir/seafloor_age/globalage_1.6.grd agegrid=$datadir/seafloor_age/age.3.6.grd #plate_grids="$datadir/gsrm/kreemer/vx.1.grd $datadir/gsrm/kreemer/vy.1.grd" plate_grids="$datadir/plates/morvel.vx.0.25.-1.grd $datadir/plates/morvel.vy.0.25.-1.grd" # NNR-MOIRVEL56 #quake_file=$datadir/quakes/engdahl/hypocenters.xyzm quake_file=$datadir/quakes/engdahl/EHB.xyzmw echo $0: using $topo_file for topo echo $0: using $agegrid for seafloor age echo $0: using $plate_grids for plate velocities echo $0: using $quake_file for seismicity for view in $views;do for cmt in $cmts;do rm .gmt* 2> /dev/null gmtset PAPER_MEDIA dj755 fsize1=68;portrait=""; offsets='-X1.0 -Y1.0' # # scales # f1=0.01 # quakes f2=0.15 # volcanoes if [ $addage -eq 1 ];then ann="" # annotation else ann=-Bg30/g30 # annotation fi if [ $view -eq 0 ];then projection=-JN170/$fsize1 # PAC projection odir=PAC elif [ $view -eq 1 ];then projection=-JN10/$fsize1 # AFR projection odir=AFR fi region=-R-180/180/-90/90 data_region=$region if [ $cmt -eq 1 ];then # full CMT ps_filename=world_map_cmt.ps pdf_filename=world_map_cmt.pdf png_filename=world_map_cmt.png hcpt=hypo.cpt;bscl=100 elif [ $cmt -eq 2 ];then # shallow CMT ps_filename=world_map_cmt_shallow.ps pdf_filename=world_map_cmt_shallow.pdf png_filename=world_map_cmt_shallow.png hcpt=hypo.shallow.cpt;bscl=10 else # only engdahl ps_filename=world_map.ps pdf_filename=world_map.pdf png_filename=world_map.png hcpt=hypo.cpt fi #psbasemap $region $projection -Bg30 -P > $ps_filename ;exit # raster data set specific: # Plotting topo if [ $topo -eq 1 ];then # Use grdimage to create a raster map. if [ $shade -eq 0 ];then # no shading $gmtbin/grdimage $topo_file $region $projection \ -C$cmp $ann \ $verbose $portrait $offsets -K > $ps_filename else # shading $gmtbin/grdimage $topo_file \ -I$topo_ifile $region $projection -C$cmp $ann \ $verbose $portrait $offsets -K > $ps_filename fi else # only pscoast coastlines $gmtbin/pscoast $region $projection $ann \ -S200 -G100 -Dc -A70000 $verbose $portrait $offsets -K > $ps_filename fi if [ $addage -eq 1 ];then $gmtbin/grdcontour $agegrid -C10 $region $projection -O -K -W0.5 >> $ps_filename $gmtbin/grdcontour $agegrid -C30 -Af18 -G128 $region $projection -O -K -W3 >> $ps_filename fi # Add a scale beneath the plot. $gmtbin/psscale -C$cmp -N50 -D5/0.5/10/.5h -B2500:"topography":/:"m": \ $verbose -E -O -K >> $ps_filename if [ $cmt -ne 2 ];then $gmtbin/psscale -C$hcpt -N50 -D5/2.2/10/.5h -B$bscl:"hypocenter depth":/:"km": \ $verbose -O -K >> $ps_filename fi if [ $topo -eq 1 ];then # Use pscoast for possible features such as national boundaries. $gmtbin/pscoast -W0.25/128 $region $projection -Df $portrait \ -O -K $verbose >> $ps_filename fi # plate velocoties according to kreemer NNR grdvector -I7.5 $plate_grids \ $region $projection -T -S5 -Q0.15/0.2/0.15n.5 -W0.5 -G192/147/32 -O -K $verbose >> $ps_filename # Add slab contours to the plot. #$gmtbin/psxy /home/becker/progs/src/igmt_develop/allslabs_rum.gmt $verbose \ # -M $region $projection -O -K -W1/000/000/000 >> $ps_filename # Add significant USGS/NEIC 73-98 >5 #gawk '{if((substr($1,1,1)!="#")&&($1!=""))print($7,$6,$9*f)}' \ # f=$f1 /home/becker/data/quakes/usgs_neic.dat | \ # $gmtbin/psxy $verbose $region $projection -G000/000/000 -Sc -K -O >> $ps_filename if [ $cmt -eq 1 ];then # engdahl gawk '{print($1,$2,$3,$4*f)}' f=$f1 $quake_file | \ psxy $verbose $region $projection -C$hcpt -Sc -K -O >> $ps_filename # and CMT size=0.1 psmeca -Z$hcpt -D0/800 -Sd$size/-1 $region $projection \ $datadir/cmt/all.mdat -O -K >> $ps_filename elif [ $cmt -eq 2 ];then # shallow engdahl, gray gawk '{if($3<=50)print($1,$2,$4*f)}' f=$f1 $quake_file | \ psxy $verbose $region $projection -Ggray -W0.5 -Sc -K -O >> $ps_filename # and CMT, colored by style size=0.1;dr=-D0/50 psmeca $dr -Sd$size/-1 $region $projection \ $datadir/cmt/all.norm.mdat -O -K -G0/0/255 >> $ps_filename psmeca $dr -Sd$size/-1 $region $projection \ $datadir/cmt/all.ss.mdat -O -K -G0/128/0 >> $ps_filename psmeca $dr -Sd$size/-1 $region $projection \ $datadir/cmt/all.thrust.mdat -G255/0/0 -O -K >> $ps_filename else # engdahl gawk '{print($1,$2,$3,$4*f)}' f=$f1 $quake_file | \ psxy $verbose $region $projection -C$hcpt -Sc -K -O -G0.5 >> $ps_filename fi # Add volcano locations from the Smithsonian. gawk 'BEGIN{FS=","}{if(NR>2)print($9,$8)}' $datadir/volcanoes/GVP_Volcano_List.csv | \ gawk '{print($1,$2,f)}' f=$f2 | \ $gmtbin/psxy -fg $verbose $region $projection -G$volc_col -W0.5 -St -O -K >> $ps_filename # legend pstext references.$cmt.txt -R0/1/0/1 -JX$fsize1 -M -O -K -N >> $ps_filename # # logos # # psimage usc.2.ras -E100 -Xa0.2 -Ya31.5 -O -K >> $ps_filename # psimage usc.1.ras -E100 -Xa66 -Ya31.5 -O -K >> $ps_filename # Use psbasemap for basic map layout, possible title # and complete the plot $gmtbin/psbasemap $region $verbose $projection \ -B:."": -O >> $ps_filename gmtset PAPER_MEDIA letter+ # reset # convert PS to PDF and PNG eps2eps $ps_filename tmp.eps mv tmp.eps $ps_filename epstopdf $ps_filename gzip -f $ps_filename echo $ps_filename $shade $view $odir mkdir $HOME/public_html/ftp/maps/$shade/$odir/ 2> /dev/null mv *.ps.gz *.pdf $HOME/public_html/ftp/maps/$shade/$odir/ /usr/bin/convert -rotate 90 \ -background white -flatten $HOME/public_html/ftp/maps/$shade/$odir/$pdf_filename \ $HOME/public_html/ftp/maps/$shade/$odir/$png_filename echo $0: output in $HOME/public_html/ftp/maps/$shade/$odir/$png_filename done done