3-10-大地电磁视电阻率曲线绘制

20 Jul 2018

#!/bin/bash

DAT=$1
SITE=`basename $DAT .asc`
PS=${SITE}.ps




#gmt set MAP_GRID_PEN_PRIMARY 0.1p,green,-
#gmt set MAP_FRAME_TYPE plain
gmt set MAP_FRAME_WIDTH 2p
gmt set FONT_LABEL 12p

gmt pstext -JX12c/1c -R10/20/10/20 -Y23c -P -K <<EOF >$PS
15 15 SITE:$SITE
EOF

# find max and min value, there is two ways
#Y_MIN=`gawk '{print $3"\n"$4}' $DAT|gawk 'NR==1{min=$1;next}{min=min<$1?min:$1}END{print min*0.1}'`
#Y_MAX=`gawk '{print $3"\n"$4}' $DAT|gawk 'NR==1{max=$1;next}{max=max>$1?max:$1}END{print max*10}'`
gawk '{print $3"\n"$4}' $DAT|gmt info -C >min_max.txt
Y_MIN=`gawk '{print $1*0.1}' min_max.txt`
Y_MAX=`gawk '{print $2*10}' min_max.txt`
#
JR=-JX12cl/10cl
JR_text=-JX12c/10c
RR=-R1/1000000/$Y_MIN/$Y_MAX

JP=-JX12cl/5c
RP=-R1/1000000/-9.9/99.9

#Tx_MIN1=`cat $DAT|gawk 'NR==1{min=$30;next}{min=min<$30?min:$30}END{print min}'`
#Tx_MAX1=`cat $DAT|gawk 'NR==1{max=$30;next}{max=max>$30?max:$30}END{print max}'`
#Tx_MIN=`cat $DAT|gawk 'NR==1{min=$Tx_MIN1;next}{min=min<$31?min:$31}END{print min*0.1}'`
#Tx_MAX=`cat $DAT|gawk 'NR==1{max=$Tx_MAX1;next}{max=max>$31?max:$31}END{print max*10}'`
JTx=-JX12cl/3c
JTx_text=-JX12c/3c
RTx=-R1/1000000/-0.6/0.6

gawk '{print 1/$1,$3,2*$11}' $DAT | gmt psxy -Y-10c -Bwsne -Ba1f1g3p $JR $RR -Ey4p -Ss0.1c -W1p,blue \
	--MAP_GRID_PEN_PRIMARY=0.1p,lightgray,- --MAP_FRAME_TYPE=inside  -K -O >> $PS

gawk '{print 1/$1,$4,2*$12}' $DAT | gmt psxy -BWsne -Ba1g1p -By+l"Apparent Resisvity(Ohm-m)" $JR $RR -Ey4p -Sc0.1c -W1p,red \
	--MAP_GRID_PEN_PRIMARY=0.1p,green,- --MAP_TICK_LENGTH=0p --MAP_FRAME_TYPE=plain -K -O >> $PS
gmt pslegend -DjTR+w2c/1c+o0.2c/0.2c  $JR_text $RR  -K -O << EOF >> $PS
S 0.4c s 0.1c  - 1p,blue 0.8c @~\162@~@-xy@-
S 0.4c c 0.1c  - 1p,red 0.8c @~\162@~@-yx@-
EOF


gawk '{print 1/$1,$7,2*$15}' $DAT | gmt psxy -Y-5.1c -Bwsne -Bxa1f1g3p -Bya20g20 $JP $RP -Ey4p/2p,yellow -Ss0.1c -W1p,blue \
	--MAP_GRID_PEN_PRIMARY=0.1p,lightgray,- --MAP_TICK_LENGTH=3p --MAP_FRAME_TYPE=inside -K -O >> $PS

gawk '{print 1/$1,$8+180,2*$16}' $DAT | gmt psxy -BWsne -Bxa1g1p -Bya20g20+l"Phase(Degrees)" $JP $RP -Ey4p -Sc0.1c -W1p,red \
	--MAP_GRID_PEN_PRIMARY=0.1p,green,- --MAP_TICK_LENGTH=0p --MAP_FRAME_TYPE=plain -K -O >> $PS

gawk '{print 1/$1,$30,2*$27}' $DAT | gmt psxy -Y-3.1c -Bwsne -Bxa1f1g3p -Bya0.5g0.5 $JTx $RTx -Ey4p -SS0.1c -W1p,blue \
	--MAP_GRID_PEN_PRIMARY=0.1p,lightgray,- --MAP_TICK_LENGTH=3p --MAP_FRAME_TYPE=inside -K -O >> $PS

gawk '{print 1/$1,$31,2*$27}' $DAT | gmt psxy -BWsne -Bxa1g1p -Bya0.5g0.5+l"Tzx" $JTx $RTx -Ey4p -Sc0.1c -W1p,red \
	--MAP_GRID_PEN_PRIMARY=0.1p,green,- --MAP_TICK_LENGTH=0p --MAP_FRAME_TYPE=plain -K -O >> $PS
gmt pslegend -DjTR+w2c/1c  $JTx_text $RTx  -K -O << EOF >> $PS
S 0.4c s 0.1c -  1p,blue 0.8c real
S 0.4c c 0.1c  - 1p,red 0.8c img
EOF

gawk '{print 1/$1,$32,2*$28}' $DAT | gmt psxy -Y-3.1c -Bwsne -Bxa1f1g3p -Bya0.5g0.5 $JTx $RTx -Ey4p -Ss0.1c -W1p,blue \
	--MAP_GRID_PEN_PRIMARY=0.1p,lightgray,- --MAP_TICK_LENGTH=3p --MAP_FRAME_TYPE=inside -K -O >> $PS

gawk '{print 1/$1,$33,2*$28}' $DAT | gmt psxy -BWSne -Bxa1g1p+l"Period(secs)" -Bya0.5g0.5+l"Tzy" $JTx $RTx -Ey4p -Sc0.1c -W1p,red \
	--MAP_GRID_PEN_PRIMARY=0.1p,green,- --MAP_TICK_LENGTH=0p --MAP_FRAME_TYPE=plain -O >> $PS

更多资料
GMT目录