Skip to content

Commit 5432831

Browse files
danielpeterdanielpeter
committed
updates plotting script
1 parent 1b4646e commit 5432831

1 file changed

Lines changed: 54 additions & 34 deletions

File tree

‎utils/Visualization/GMT/plot_wavefield.pl‎

Lines changed: 54 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -236,7 +236,7 @@
236236
print CSH "echo\n";
237237

238238
if($gmt5==1){
239-
print CSH "gmtset MAP_FRAME_TYPE plain PS_MEDIA letter PROJ_LENGTH_UNIT inch FORMAT_GEO_MAP D MAP_TICK_LENGTH_PRIMARY $tick PS_CHAR_ENCODING Standard+ COLOR_NAN $cgray\n";
239+
print CSH "gmt gmtset MAP_FRAME_TYPE plain PS_MEDIA letter PROJ_LENGTH_UNIT inch FORMAT_GEO_MAP D MAP_TICK_LENGTH_PRIMARY $tick PS_CHAR_ENCODING Standard+ COLOR_NAN $cgray\n";
240240
}else{
241241
# older gmt
242242
print CSH "gmtset BASEMAP_TYPE plain PAPER_MEDIA letter MEASURE_UNIT inch PLOT_DEGREE_FORMAT D TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 CHAR_ENCODING Standard+ COLOR_NAN $cgray\n";
@@ -260,7 +260,7 @@
260260
$Ts[$k] = sprintf("-T%3.3f/%3.3f/%3.3f",-$ss[$k]*1.05,$ss[$k]*1.05,$ds[$k]);
261261
print "color range Ts = $Ts[$k]\n";
262262

263-
print CSH "makecpt -C$colorbar $Ts[$k] -D > color_${k}.cpt\n";
263+
print CSH "gmt makecpt -C$colorbar $Ts[$k] -D > color_${k}.cpt\n";
264264
}
265265

266266
# colorbar
@@ -291,13 +291,14 @@
291291

292292
$name = "wavefield_${tlab}_${ttag}";
293293
$psfile = "$name.ps";
294+
$pdffile = "$name.pdf";
294295
$jpgfile = "$name.jpg";
295296

296297
print CSH "echo\n";
297298
print CSH "echo \"postscript file = $psfile\";echo\n";
298299
print CSH "echo\n";
299300
# clean
300-
print CSH "rm -f $psfile\n\n";
301+
print CSH "rm -f $psfile $pdffile $jpgfile\n\n";
301302

302303
#$numf = 1;
303304

@@ -395,28 +396,28 @@
395396
#} # temp
396397

397398
if ($i == $imin && $k==1) {
398-
print CSH "psbasemap $J $R $B $BG -K -V $orient $origin > $psfile\n";
399+
print CSH "gmt psbasemap $J $R $B $BG -K -V $orient $origin > $psfile\n";
399400
} # START
400401
else {
401-
print CSH "psbasemap $J $R $B $BG -K -O -V $shift >> $psfile\n";
402+
print CSH "gmt psbasemap $J $R $B $BG -K -O -V $shift >> $psfile\n";
402403
}
403404

404405
# PLOT THE FORWARD WAVEFIELD
405406
if ($icolor==1) {
406407
if ($igrd==1) {
407408
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f > dfile\n";
408409
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f | nearneighbor -G$grdfile $R $interp\n";
409-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f | xyz2grd -G$grdfile $R $interp\n";
410-
print CSH "grdimage $grdfile -C$cfile $J -K -O -V -Q >> $psfile\n";
410+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f | gmt xyz2grd -G$grdfile $R $interp\n";
411+
print CSH "gmt grdimage $grdfile -C$cfile $J -K -O -V -Q >> $psfile\n";
411412
} else {
412-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f | psxy $R $J $minfo -C$cfile -K -O -V >> $psfile\n";
413+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $nm}' $snapshot_f | gmt psxy $R $J $minfo -C$cfile -K -O -V >> $psfile\n";
413414
}
414415

415416
# mask points above topography
416417
if ($imask==1) {
417418
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' $snapshot_f > maskpts\n";
418419
printf CSH "grep NaN $snapshot_f > maskpts\n";
419-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' maskpts | psxy $R $J $minfo -C$cfile -K -O -V >> $psfile\n";
420+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' maskpts | gmt psxy $R $J $minfo -C$cfile -K -O -V >> $psfile\n";
420421
}
421422
}
422423

@@ -429,7 +430,7 @@
429430
die("check if bfile $bfile exist or not\n");
430431
}
431432
#print CSH "awk '{print \$1/1000,\$2/1000}' $bfile | psxy $J $R -W1.5p -m -K -O -V >> $psfile\n";
432-
print CSH "awk '{print \$1,\$2}' $bfile | psxy $J $R -W1p -m -K -O -V >> $psfile\n";
433+
print CSH "awk '{print \$1,\$2}' $bfile | gmt psxy $J $R -W1p -m -K -O -V >> $psfile\n";
433434
}
434435

435436
# plot seismicity
@@ -438,11 +439,11 @@
438439
if (not -f $bfile) {
439440
die("check if bfile $bfile exist or not\n");
440441
}
441-
print CSH "awk '{print \$1,\$2}' $bfile | psxy $J $R -Sc3p -G0 -K -O -V >> $psfile\n";
442+
print CSH "awk '{print \$1,\$2}' $bfile | gmt psxy $J $R -Sc3p -G0 -K -O -V >> $psfile\n";
442443
}
443444

444-
if ($i == $imin) {print CSH "psscale -C$cfile $Dscale $BscaleS -K -O -V >> $psfile \n";}
445-
print CSH "psbasemap $J $R $B -K -O -V >> $psfile\n";
445+
if ($i == $imin) {print CSH "gmt psscale -C$cfile $Dscale $BscaleS -K -O -V >> $psfile \n";}
446+
print CSH "gmt psbasemap $J $R $B -K -O -V >> $psfile\n";
446447

447448
# label the figures
448449
if (0==1) {
@@ -453,19 +454,19 @@
453454
print CSH "pstext $J $R -K -O -V >>$psfile<<EOF\n -1 0.75 14 0 $fontno CM ice\nEOF\n";
454455
print CSH "pstext $J $R -K -O -V >>$psfile<<EOF\n 1.0 0.75 14 0 $fontno CM water\nEOF\n";
455456
}
456-
print CSH "psxy -W1p $J $R -K -O -V >>$psfile<<EOF\n1.0 0.6\n 0.7 -0.1\nEOF\n";
457-
print CSH "psxy -W1p $J $R -K -O -V >>$psfile<<EOF\n-1.0 0.6\n -0.7 -0.05\nEOF\n";
457+
print CSH "gmt psxy -W1p $J $R -K -O -V >>$psfile<<EOF\n1.0 0.6\n 0.7 -0.1\nEOF\n";
458+
print CSH "gmt psxy -W1p $J $R -K -O -V >>$psfile<<EOF\n-1.0 0.6\n -0.7 -0.05\nEOF\n";
458459
}
459460
}
460461

461462
# plot stations and labels
462463
if($iplotrec==1) {
463464
#print CSH "awk '{print \$3/1000,\$4/1000}' $recfile2 | psxy $J $R -K -O -V $rec2 >> $psfile\n";
464-
print CSH "awk '{print \$3/1000,\$4/1000}' $recfile | psxy $J $R -K -O -V $rec >> $psfile\n";
465+
print CSH "awk '{print \$3/1000,\$4/1000}' $recfile | gmt psxy $J $R -K -O -V $rec >> $psfile\n";
465466
#if($i==$imin && $iplotreclab==1) {print CSH "awk '{print \$3/1000,\$4/1000,8,0,$fontno,\"CM\",\"S\"\$7}' $recfile | pstext $textrec $J $R -K -O -V >> $psfile\n";}
466-
if($i==$imin && $iplotreclab==1 && $k==1) {print CSH "awk '{print \$3/1000,\$4/1000,8,0,$fontno,\"CM\",\$1}' $recfile | pstext $textrec $J $R -K -O -V >> $psfile\n";}
467+
if($i==$imin && $iplotreclab==1 && $k==1) {print CSH "awk '{print \$3/1000,\$4/1000,8,0,$fontno,\"CM\",\$1}' $recfile | gmt pstext $textrec $J $R -K -O -V >> $psfile\n";}
467468
}
468-
print CSH "psxy $J $R -K -O -V $src >> $psfile<<EOF\n $srcx $srcz\nEOF\n";
469+
print CSH "gmt psxy $J $R -K -O -V $src >> $psfile<<EOF\n $srcx $srcz\nEOF\n";
469470

470471
# plot the time of the snapshot (for some reason, it won't work inside the B command)
471472
#$xtext = $xmin-0.3*$xran;
@@ -532,7 +533,7 @@
532533
if($ipsv==1) {$tlab0 = "PSV"} else {$tlab0 = "SH"}
533534
$title = "$tlab0 wavefield -- $tlab";
534535
if($gmt5==1){
535-
print CSH "pstext -N $J $R $Utag -K -O -V $shift -F+f+a+j >>$psfile<<EOF\n $xmin $zmax $fsize0,$fontno 0 LM $title\nEOF\n";
536+
print CSH "gmt pstext -N $J $R $Utag -K -O -V $shift -F+f+a+j >>$psfile<<EOF\n $xmin $zmax $fsize0,$fontno 0 LM $title\nEOF\n";
536537
}else{
537538
print CSH "pstext -N $J $R $Utag -K -O -V $shift >>$psfile<<EOF\n $xmin $zmax $fsize0 0 $fontno LM $title\nEOF\n"; # LM or CM
538539
}
@@ -570,14 +571,15 @@
570571

571572
$name = "kernel_${tlab}_${ttag}";
572573
$psfile = "$name.ps";
574+
$pdffile = "$name.pdf";
573575
$jpgfile = "$name.jpg";
574576

575577
print CSH "echo\n";
576578
print CSH "echo \"kernel file = $kfile\"\n";
577579
print CSH "echo \"postscript file = $psfile\";echo\n";
578580
print CSH "echo\n";
579581
# clean
580-
print CSH "rm -f $psfile\n\n";
582+
print CSH "rm -f $psfile $pdffile $jpgfile\n\n";
581583

582584
if($gmt5==1){
583585
$BscaleS = sprintf("-B%2.2e+l\"%s (x, z, t), 10\@+%2.2i\@+ m\@+-2\@+\" -BS",$bs[0],"K",$pwr[0]);
@@ -619,9 +621,9 @@
619621

620622
# START
621623
if ($k == $kmin) {
622-
print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n";
624+
print CSH "gmt psbasemap $J $R $B -K -V $orient $origin > $psfile\n";
623625
} else {
624-
print CSH "psbasemap $J $R $B -K -O -V $dY >> $psfile\n";
626+
print CSH "gmt psbasemap $J $R $B -K -O -V $dY >> $psfile\n";
625627
}
626628

627629
print CSH "echo;echo \"norm = $norm[$k]\"\n";
@@ -632,20 +634,20 @@
632634
print CSH "echo;echo \"plotting by grdimage\";echo;\n";
633635
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile > dfile\n";
634636
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile | nearneighbor -G$grdfile $R $interp\n";
635-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile | xyz2grd -G$grdfile $R $interp\n";
636-
print CSH "grdinfo $grdfile\n";
637-
print CSH "grdimage $grdfile -Ccolor_${k}.cpt $J -K -O -V -Q >> $psfile\n";
637+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile | gmt xyz2grd -G$grdfile $R $interp\n";
638+
print CSH "gmt grdinfo $grdfile\n";
639+
print CSH "gmt grdimage $grdfile -Ccolor_${k}.cpt $J -K -O -V -Q >> $psfile\n";
638640
} else {
639641
print CSH "echo;echo \"plotting by psxy\";\n";
640-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile | psxy $R $J $minfo -Ccolor_${k}.cpt -K -O -V >> $psfile\n";
642+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2) / $norm[$k]}' $kfile | gmt psxy $R $J $minfo -Ccolor_${k}.cpt -K -O -V >> $psfile\n";
641643
}
642644

643645
# mask points above topography
644646
if($imask==1) {
645647
print CSH "echo;echo \"using mask\";\n";
646648
#print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' $kfile > maskpts\n";
647649
printf CSH "grep NaN $kfile > maskpts\n";
648-
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' maskpts | psxy $R $J $minfo -Ccolor_${k}.cpt -K -O -V >> $psfile\n";
650+
print CSH "awk '{print \$1/1000,\$2/1000,\$($comp+2)}' maskpts | gmt psxy $R $J $minfo -Ccolor_${k}.cpt -K -O -V >> $psfile\n";
649651
}
650652
}
651653

@@ -659,13 +661,13 @@
659661
die("check if bfile $bfile exist or not\n");
660662
}
661663
#print CSH "awk '{print \$1/1000,\$2/1000}' $bfile | psxy $J $R -W1.5p -m -K -O -V >> $psfile\n";
662-
print CSH "awk '{print \$1,\$2}' $bfile | psxy $J $R -W1p -m -K -O -V >> $psfile\n";
664+
print CSH "awk '{print \$1,\$2}' $bfile | gmt psxy $J $R -W1p -m -K -O -V >> $psfile\n";
663665
}
664666

665667
#print CSH "pscoast $J $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
666668
#print CSH "awk '{print \$2,\$1}' INPUT/oms_shelf |psxy $J $R $Wshelf -K -O -V >> $psfile\n";
667669
if ($k == $kmin) {
668-
print CSH "psscale -Ccolor_${k}.cpt $Dscale $BscaleS -K -O -V >> $psfile \n";
670+
print CSH "gmt psscale -Ccolor_${k}.cpt $Dscale $BscaleS -K -O -V >> $psfile \n";
669671
}
670672
#print CSH "awk '{print \$3/1000,\$4/1000}' $recfile | psxy $J $R -K -O -V $rec >> $psfile\n";
671673
#print CSH "psxy -N $J $R -K -O -V $src >> $psfile<<EOF\n $srcx $srcz\nEOF\n";
@@ -692,7 +694,7 @@
692694
$textinfo = "-G255"; #"-G0 -W255 -C2p -N";
693695
print "\ntitle plotting at $xtx, $ztx info: $textinfo $titles[$k]\n";
694696
if($gmt5==1){
695-
print CSH "pstext $textinfo $J $R -K -O -V -F+f+a+j >>$psfile<<EOF\n $xtx $ztx 16,$fontno 0 RT $titles[$k]\nEOF\n";
697+
print CSH "gmt pstext $textinfo $J $R -K -O -V -F+f+a+j >>$psfile<<EOF\n $xtx $ztx 16,$fontno 0 RT $titles[$k]\nEOF\n";
696698
}else{
697699
# older gmt
698700
print CSH "pstext $textinfo $J $R -K -O -V >>$psfile<<EOF\n $xtx $ztx 16 0 $fontno RT $titles[$k]\nEOF\n";
@@ -712,7 +714,7 @@
712714
if ($ipsv==1) {$tlab0 = "PSV";} else {$tlab0 = "SH";}
713715
$title = "$tlab0 kernel -- $tlab";
714716
if($gmt5==1){
715-
print CSH "pstext -N $J $R $Utag -K -O -V $shift -F+f+a+j >>$psfile<<EOF\n $xmin $zmax $fsize0,$fontno 0 LM $title\nEOF\n";
717+
print CSH "gmt pstext -N $J $R $Utag -K -O -V $shift -F+f+a+j >>$psfile<<EOF\n $xmin $zmax $fsize0,$fontno 0 LM $title\nEOF\n";
716718
}else{
717719
# older gmt
718720
print CSH "pstext -N $J $R $Utag -K -O -V $shift >>$psfile<<EOF\n $xmin $zmax $fsize0 0 $fontno LM $title\nEOF\n";
@@ -725,14 +727,32 @@
725727
#-------------------------
726728
# FINISH
727729
#print CSH "pstext $J -R0/1/0/1 -O -V >>$psfile<<EOF\n 10 10 $fsize0 0 $fontno CM junk \nEOF\n";
728-
print CSH "psxy $J -R0/1/0/1 -O -V >>$psfile<<EOF\nEOF\n\n";
730+
print CSH "gmt psxy $J -R0/1/0/1 -O -V >>$psfile<<EOF\nEOF\n\n";
729731
# cleanup
730732
print CSH "rm -f temp.grd color_0.cpt color_1.cpt color_2.cpt\n";
731733
# user output
732734
print CSH "\necho\necho \"see output psfile: $psfile\"\necho\n";
733735

734-
# converts (using ImageMagick command)
735-
#print CSH "convert $psfile $jpgfile\n";
736+
# convert ps to pdf
737+
# check if ps2pdf command is available
738+
my $ps2pdf = `which ps2pdf 2>/dev/null`;
739+
chomp($ps2pdf);
740+
if (!defined $ps2pdf || $ps2pdf eq '') {
741+
print "ps2pdf command not found. Please install Ghostscript to convert ps-file to pdf-file.\n";
742+
}else{
743+
print CSH "echo \"converting to pdf...\"\n";
744+
print CSH "ps2pdf $psfile $pdffile\n";
745+
}
746+
# converts ps to jpg (using ImageMagick command)
747+
my $magick = `which magick 2>/dev/null`;
748+
chomp($magick);
749+
if (!defined $magick || $magick eq '') {
750+
print "magick command (ImageMagick) not found. Please install ImageMagick to convert ps-file to jpg-file.\n";
751+
}else{
752+
print CSH "echo \"converting to jpg...\"\n";
753+
print CSH "magick -density 300 -quality 90 $psfile $jpgfile\n";
754+
}
755+
736756
close (CSH);
737757

738758
# executes script

0 commit comments

Comments
 (0)