GMTの図
古気候データのグラフ作成
#!/bin/sh
awk -f ngeo.awk ngeo1186-s1.prn > ngeo1186-s1.dat
proj=x-0.4/11l
region=0/66/80/3000
gmtset ANNOT_FONT_SIZE_PRIMARY 20
gmtset ANNOT_FONT_PRIMARY 0
gmtset TICK_LENGTH -0.2c
f1=ngeo1186-s1.dat
#proxyによって色わけ
awk '{if($1=="Boron") print $2,$5,0.5}' $f1 >tmp
awk '{if($1=="Phytoplankton") print $2,$5,1.5}' $f1 >>tmp
awk '{if($1=="Stomata") print $2,$5,2.5}' $f1 >>tmp
awk '{if($1=="Paleosols") print $2,$5,3.5}' $f1 >>tmp
awk '{if($1=="B/Ca") print $2,$5,4.5}' $f1 >>tmp
awk '{if($1=="Liverworts") print $2,$5,5.5}' $f1 >>tmp
awk '{if($1=="Nahcolite") print $2,$5,6.5}' $f1 >>tmp
echo 0 50 50 255 1 50 50 255 >tmp.cpt
echo 1 255 0 0 2 255 0 0 >>tmp.cpt
echo 2 0 128 0 3 0 128 0 >>tmp.cpt
echo 3 129 129 129 4 129 129 129 >>tmp.cpt
echo 4 239 143 15 5 239 143 15 >>tmp.cpt
echo 5 243 213 16 6 243 213 16 >>tmp.cpt
echo 6 0 0 0 7 0 0 0 >>tmp.cpt
#エラーバーを取得
awk -f getuncco2line.awk $f1 > tmp2
awk -f getunctimeline.awk $f1 >> tmp2
awk '{if($1=="Nahcolite" && $5<-1.0) print $2,$6,6.5}' $f1 >>tmp
awk '{if($1=="Stomata" && $5<-1.0) print $2,$6,2.5}' $f1 >>tmp
awk -f getunctimelowerco2line.awk $f1 >> tmp2
psxy tmp2 -J$proj -R$region -K -M -W1,black -X2 -Y1.5>test5.eps
cs=c0.35
psxy tmp -J$proj -R$region -K -S$cs -Ctmp.cpt -W1,black -BSWnea10/a3 -O>>test5.eps
# 凡例をかきこんでいく
echo ">" > tmp2
echo 20 3000 >>tmp2
echo 20 1000 >>tmp2
echo ">" >> tmp2
echo 0 1000 >>tmp2
echo 20 1000 >>tmp2
psxy tmp2 -J$proj -R$region -K -M -W2,black -O>>test5.eps
echo 10 10 0.5 >tmp
echo 10 11 1.5 >>tmp
echo 10 12 2.5 >>tmp
echo 10 13 3.5 >>tmp
echo 10 14 4.5 >>tmp
echo 10 15 5.5 >>tmp
echo 10 16 6.5 >>tmp
psxy tmp -Jx2/0.67 -R0/20/0/20 -K -O -Sc0.4 -W1,black -Ctmp.cpt -Y6 -X-0.5>>test5.eps
fo1=16
fo2=0
echo 10 10 $fo1 0 $fo2 "ML" Boron >tmp
echo 10 11 $fo1 0 $fo2 "ML" Phytoplankton >>tmp
echo 10 12 $fo1 0 $fo2 "ML" Stomata >>tmp
echo 10 13 $fo1 0 $fo2 "ML" Paleosols >>tmp
echo 10 14 $fo1 0 $fo2 "ML" B/Ca >>tmp
echo 10 15 $fo1 0 $fo2 "ML" Liverworts >>tmp
echo 10 16 $fo1 0 $fo2 "ML" Nahcolite >>tmp
pstext -J -R -K -O -N tmp -D0.5/0>>test5.eps
ngeo.awkの中身は
{
printf("%s\t", substr($0, 1,13)) #proxy
printf("%s\t", substr($0, 53, 4)) #LON
if(substr($0, 69,1) ~ /^[0-9]/) {
printf("%s\t", substr($0, 66, 4))}
else {printf("%s\t", -999)}#
if(substr($0, 83,1) ~ /^[0-9]/) {
printf("%s\t", substr($0, 80, 4))}
else {printf("%s\t", -999)}
if(substr($0, 92,1) ~ /^[0-9]/) {
printf("%s\t", substr($0, 89, 4))}
else {printf("%s\t", -999)}
if(substr($0,107,1) ~ /^[0-9]/) {
printf("%s\t", substr($0, 104, 4))}
else {printf("%s\t", -999)}
if(substr($0,122,1) ~ /^[0-9]/) {
printf("%s\t", substr($0,119, 4))}
else {printf("%s\t", -999)}
print ""
}
getuncco2line.awk
{
if($5>80.0 && $6>-1.0) {
print ">"
print $2,$6
print $2,$7
}
}
getunctimeline.awk
{
if($5>80.0 && $3>-1.0) {
print ">"
print $3,$5
print $4,$5
}
}
getunctimelowerco2line.awk(lower limitのみ示したデータ)
{
if($5<-1.0) {
print ">"
print $2,$6
print $2,3000
print ">"
print $3,$6
print $4,$6
}
}
TSダイアグラム
pw=999.842594 + 6.793952*1e-2*T(i,j,k) - 9.095290*1e-3*T(i,j,k)**2.0 &
+ 1.001685*1e-4*T(i,j,k)**3.0 - 1.120083*1e-6*T(i,j,k)**4.0 + 6.536332*1e-9*T(i,j,k)**5.0
ro(i,j,k)=pw+ S(i,j,k)*( 0.824493 + (-4.0899*1e-3)*T(i,j,k) + (7.6438*1e-5)*T(i,j,k)**2.0 &
+ (-8.2467*1e-7)*T(i,j,k)**3.0 + (5.3875*1e-9)*T(i,j,k)**4.0) &
+ S(i,j,k)**1.5* (-5.72466*1e-3 + (1.0227*1e-4)*T(i,j,k) + (-1.6546*1e-6)*T(i,j,k)**2.0) &
+ S(i,j,k)**2.0*(4.8314*1e-4)
でてくるのはnewtral density anomalyです。
温度塩分をとりあえず0.05きざみでかえた(塩分,温度,密度)のテーブルdensity.datを作って、
xyz2grd -R28/38/-2/32 -I0.05/0.05 -Gdensity.grd density.dat
とすればテーブルができるので、
grdcontour -J -R -C0.5 -A1 density.grd
で描画できます。
もうつくってある。→density.grd
上の図のスクリプトは、
#!/bin/sh
proj=x2.8/0.65
region=28/39/-2/36
wid=8
sc=0.4
psxy -J$proj -R$region tmp -K -W$wid,blue>test.eps
awk '{if($4==11 || $4==23 || $4==39 || $4==55 || $4==72 || $4==97 || $4==121 || $4==138 || $4==154 || $4==170 || $4==182) print $1,$2}' tmp > tmp0
psxy -J$proj -R$region tmp0 -K -Sc$sc -W$wid,blue -Gwhite -O -BSWnea1/a4>>test.eps
grdcontour -J -R -K -O -C0.5 -A1 density.grd>>test.eps
南北熱輸送
#!/bin/sh
proj=x0.2/1.7
region=-90/90/-7/7
a1=ATHTTP_ANNUAL
ao=AOHTTP_ANNUAL
o2=OHTTP_ANNUAL
a5=WVTP_ANNUAL
out=test3.eps
rm -f $out
gmtset TICK_LENGTH -0.1c
gmtset ANNOT_FONT 0
gmtset ANNOT_FONT_SIZE_PRIMARY 28
gmtset GRID_PEN_PRIMARY 0.01c
gmtset ANNOT_OFFSET_PRIMARY 0.2c
wid=15
wid2=10
awk '{print $1,$2*1e-15} ' $a1 | psxy -J$proj -R$region -K -BSWneg90a15/g10a1 -W$wid,red -X3.5 -Y3>$out
awk '{print $1,$2*1e-15} ' $ao | psxy -J$proj -R$region -K -W$wid,black -O>>$out
awk '{print $1,$2*1e-15} ' $a5 | psxy -J$proj -R$region -K -W$wid,0/64/0 -O>>$out
awk '{print $1,$2*1e-15} ' $o2 | psxy -J$proj -R$region -K -W$wid,blue -O>>$out
awk '{print $1,$2*1e-15} ' $a1.ctl | psxy -J$proj -R$region -K -BSWneg90a15/g10a1 -W$wid2,red,.. -O>>$out
awk '{print $1,$2*1e-15} ' $ao.ctl | psxy -J$proj -R$region -K -W$wid2,black,.. -O>>$out
awk '{print $1,$2*1e-15} ' $a5.ctl | psxy -J$proj -R$region -K -W$wid2,0/64/0,.. -O>>$out
awk '{print $1,$2*1e-15} ' $o2.ctl | psxy -J$proj -R$region -K -W$wid2,blue,.. -O>>$out
echo -102 0 28 90 1 "MC" "Meridional Heat Transport [PW]" | pstext -J -R -K -O -N >>$out
echo 0 -8 28 0 1 "MC" "Latitude" | pstext -J -R -K -O -N >>$out
echo -86 6.5 24 0 1 "ML" ALL heat transport | pstext -J -R -K -O -N -Gblack -Y0.3>>$out
echo -86 6.0 24 0 1 "ML" Ocean heat transport | pstext -J -R -K -O -N -Gblue>>$out
echo -86 5.5 24 0 1 "ML" Atmosphere heat transport | pstext -J -R -K -O -N -Gred>>$out
echo -86 5.0 24 0 1 "ML" Latent heat transport | pstext -J -R -K -O -N -G0/128/0>>$out
echo -80 4.5 >tmp
echo -65 4.5 >>tmp
psxy tmp -J -R -K -O -N -W$wid2,black,.. >>$out
echo -62 4.5 24 0 1 "ML" "285 ppm with ice (CTL)" | pstext -J -R -K -O -N -Gblack>>$out
echo -80 4 >tmp
echo -65 4 >>tmp
psxy tmp -J -R -K -O -N -W$wid2,black>>$out
echo -62 4 24 0 1 "ML" 1142 ppm with ice | pstext -J -R -K -O -N -Gblack>>$out
海底地形地図の作成
#!/bin/sh
proj=x0.05
region=0/360/-90/90
cpt=tmp.cpt
gray=210
out=test.eps
rm -f $out
cont=1
gmtset ANNOT_FONT_SIZE_PRIMARY 20
echo 0 95 110 50 50 95 110 50 >$cpt
echo 50 133 113 55 100 133 113 55 >>$cpt
echo 100 172 152 70 200 172 152 70 >>$cpt
echo 200 224 204 90 300 224 204 90>>$cpt
echo 300 250 232 113 500 250 232 113 >>$cpt
echo 500 251 236 139 700 251 236 139 >>$cpt
echo 700 254 252 235 1000 254 252 235 >>$cpt
echo 1000 191 248 189 2000 191 248 189 >>$cpt
echo 2000 129 221 176 3000 129 221 176 >>$cpt
echo 3000 43 173 192 4000 43 173 192 >>$cpt
echo F 0 115 162 >>$cpt
f1=bathymetry
f2=coastline
f3=land
echo 180 0 | psxy -J$proj -R$region -K -X2 -Y10>$out
psxy -J$proj -R$regiond $f1 -M -L -K -O -A -C$cpt>> $out
psxy -J$proj -R$region $f3 -M -L -K -A -G$gray/$gray/$gray -O>>$out
psxy -J$proj -R$region $f2 -M -L -K -A -W3 -O -BSWnea30/a15>>$out
psscale -C$cpt -D19/4/8/0.5 -K -O -Ef0.8 -L>>$out
cp test.eps test2.eps
awk '{if($1=="%%BoundingBox:" || $1=="%%HiResBoundingBox:"){print $1, 295,$3,595,685}else{print $0}}' < test2.eps > test.eps
地図にデータを重ね合わせ
#!/bin/sh
region=672000/672800/159000/159800
proj=x0.02175/0.0214
line=SWneg100a100/g100a100
psxy -J$proj -R$region boundary20000824.dat -W8,red -K -B$line -P -X0 -Y0>test.eps
psxy -J$proj -R$region boundary2008.dat -W10,brown -K -B$line -O -P>>test.eps
echo 672700 159450 3.9 1.9 | psxy -J -R -Sr -Gwhite -K -O>>test.eps
echo 672650 159470 0.7 0.07 | psxy -J -R -Sr -Gred -K -O>>test.eps
echo 672650 159430 0.7 0.07 | psxy -J -R -Sr -Gbrown -K -O>>test.eps
echo 672730 159470 18 0 0 "MC" 2000 | pstext -J -R -Gred -K -O>>test.eps
echo 672730 159430 18 0 0 "MC" 2008 | pstext -J -R -Gbrown -K -O>>test.eps
awk '{print $2,$3}' sone.dat >tmp
psxy -J$proj -R$region tmp -Sc0.5 -Gorange -K -O -P>>test.eps
awk '{print $2,$3,14,0,0,"MC",$1}' sone.dat >tmp
pstext -J$proj -R$region tmp -K -O -P>>test.eps
size=0.25
awk '{if($4==1) print $2,$3}' 2013.dat >tmp
psxy -J$proj -R$region tmp -K -O -P -Sc0.22 -Gblack>>test.eps
awk '{if($4==2) print $2,$3}' 2013.dat >tmp
psxy -J$proj -R$region tmp -K -O -P -Sd$size -Gblack>>test.eps
awk '{if($4==3) print $2,$3}' 2013.dat >tmp
psxy -J$proj -R$region tmp -K -O -P -Si$size -Gblack>>test.eps
awk '{if($4==4) print $2,$3}' 2013.dat >tmp
psxy -J$proj -R$region tmp -K -O -P -St$size -Gblack>>test.eps
echo 672370 159300 135 2.5 1.5 | psxy -Se -J -R -K -O -W5,yellow>>test.eps
echo 672474.59 159232 | psxy -Sc0.8 -J -R -K -O -W5,yellow>>test.eps
echo 672515 159224 | psxy -Sc0.7 -J -R -K -O -W5,yellow>>test.eps
echo 672560.92 159209.00 | psxy -Sc0.6 -J -R -K -O -W5,yellow>>test.eps
echo 672370 159300 20 0 5 "MC" A | pstext -J -R -K -O -D1/1>>test.eps
echo 672474.59 159232 20 0 5 "MC" B | pstext -J -R -K -O -D0/1>>test.eps
echo 672515 159224 20 0 5 "MC" C | pstext -J -R -K -O -D0/1>>test.eps
echo 672560.92 159209.00 20 0 5 "MC" D | pstext -J -R -K -O -D0/1>>test.eps
awk '{if($1=="%%BoundingBox:" || $1=="%%HiResBoundingBox:"){print $1, $2,$3,492,484}else{print $0}}' < test.eps > test2.eps
convert test2.eps test2.png
convert rhone672000672800159000159800.png test2.png -composite tmp.png
地図にデータを重ね合わせその2
氷床モデル
#!/bin/sh
#gfortran read.f90 -L/usr/local/netcdf-4.3.0/lib -lnetcdf -L/usr/lib -lnetcdff
out=test.eps
cpt=tmp.cpt
cpt2=tmp2.cpt
region=270/32/90/47r
#region=0/360/60/90
proj=S315/90/19
echo 0 95 110 50 100 95 110 50 >$cpt
echo 100 133 113 55 200 133 113 55 >>$cpt
echo 200 172 152 70 300 172 152 70 >>$cpt
echo 300 224 204 90 400 224 204 90>>$cpt
echo 400 250 232 113 500 250 232 113 >>$cpt
echo 500 254 252 235 600 254 252 235 >>$cpt
echo 600 191 248 189 700 191 248 189 >>$cpt
echo 700 129 221 176 800 129 221 176 >>$cpt
echo 800 43 173 192 900 43 173 192 >>$cpt
echo 900 21 123 176 1000 21 123 176 >>$cpt
echo 1000 11 63 165 1500 11 63 165 >>$cpt
echo 1500 5 23 155 2000 5 23 165 >>$cpt
echo 2000 0 0 115 2500 0 0 115 >>$cpt
echo 2500 0 0 65 3000 0 0 65 >>$cpt
echo F 0 0 31 >>$cpt
echo B 255 220 220 >>$cpt
makecpt -Cglobe >$cpt
makecpt -Chot -T0/4000/200 >$cpt2
sc=0.2
i=101
gmtset GRID_PEN_PRIMARY 0.25p,,0.25_2:0
while [ $i -lt 902 ]
do
j=`expr $i - 100`
k=`expr 801 - $j`
l=`echo "scale=7 ; "$k" * 0.5" | bc`
echo $l
echo $j > time
./a.out
psxy -J$proj -R$region -K -M -L -A -C$cpt tmp2 -P -X1.2 -Y1.2> $out
psxy -J -R -K -O -M -L -A tmp1 -C$cpt2 -BSWNEg10/g5>>$out
i=`expr $i + 1`
echo 190 50 5 2 | psxy -J -R -K -O -Sr -Gwhite >>$out
echo 190 50 32 0 27 "MC" $l"ka" | pstext -J -R -K -O >>$out
cp test.eps test2.eps
awk '{if($1=="%%BoundingBox:" || $1=="%%HiResBoundingBox:"){print $1, $2,$3,595,595}else{print $0}}' < test2.eps > test.eps
convert test.eps ice$i.png
done
人文地理的図（国勢調査）
#!/bin/sh
proj=m1:125000 #5/7
region=139.875/140.11/35.79/35.97 #33=>23.5, 25=>18
gmtset GRID_PEN_PRIMARY 0.25p,,0.25_1:0
gmtset BASEMAP_TYPE plain
gmtset ANNOT_FONT_SIZE_PRIMARY 12
gmtset FRAME_PEN 1.0p
gmtset ANNOT_FONT_PRIMARY 1
gmtset GRID_PEN_PRIMARY 0.25p,,0.25_1:0
gmtset BASEMAP_TYPE plain
gmtset ANNOT_FONT_SIZE_PRIMARY 15
gmtset PLOT_DEGREE_FORMAT +D
gmtset TICK_LENGTH -0.05c
gmtset ANNOT_OFFSET_PRIMARY 0.15c
gmtset GRID_PEN_PRIMARY 0.25p,,0.25_3:0
echo 1 1 | psxy -Jx1 -R0/2/0/2 -Ss90 -G219/219/219 -K -N -Y1 -X2>test.eps
echo 1 0 0 121 100 0 0 121 > tmp.cpt
makecpt -Cpanoply -Q -T2/4.2/0.16 >>tmp.cpt
#makecpt -Cpanoply -T1/20000/1000 >tmp.cpt
echo B 255 255 255>>tmp.cpt
echo N 0 99 0 >>tmp.cpt
#echo B 0 0 99 >>tmp.cpt
#echo F 99 0 0 >>tmp.cpt
#pscoast -J$proj -R$region -C200/200/255 -K -O -Df -Ia>>test.eps
awk -f jinko.awk kashiwajinko.gmt > tmp1
awk -f jinko.awk nagareyamajinko.gmt >> tmp1
awk -f jinko.awk abikojinko.gmt >> tmp1
psxy -J$proj -R$region tmp1 -Ctmp.cpt -M -L -K -A -O -W1,black>>test.eps
psxy -J$proj -R$region roadchiba.gmt -M -K -W1,99/99/99 -O >>test.eps
psxy -J -R stations1995.gmt -M -L -A -W6,0/0/0 -O -K -Sc0.2 >>test.eps
psxy -J -R tx.dat -W6,0/0/0 -O -K -Sc0.2 >>test.eps
psxy -J -R railways19552.gmt -M -L -A -W6,49/49/49,... -O -K>>test.eps
psxy -J$proj -R$region lakes22.gmt -M -K -G169/169/169 -O -W1,black>>test.eps
gmtset ANNOT_FONT_SIZE_PRIMARY 18
gmtset ANNOT_OFFSET_PRIMARY 0.28c
psscale -D21.4/9.7/17/1 -Ctmp.cpt -K -O -L -E0.9>>test.eps
gmtset ANNOT_FONT_PRIMARY 0
gmtset GRID_PEN_PRIMARY 0.25p,,0.25_3:0
gmtset ANNOT_FONT_SIZE_PRIMARY 16
echo 140.195 35.76 15 0 1 "MC" "[MUs/km^2]" | pstext -J -R -K -O -N -BSWneg0.02a0.04/g0.02a0.02>>test.eps
cp test.eps test2.eps
awk '{if($1=="%%BoundingBox:" || $1=="%%HiResBoundingBox:"){print $1, $2,$3,595,750}else{print $0}}' < test2.eps > test.eps
convert test.eps test.png
convert -rotate 90 test.png test.png
BEDMAP2を用いた氷床断面図
#!/bin/sh
out=test.eps
p=x1.26/2
d=0/20/-2/4
gmtset FRAME_PEN 2p
gmtset PAPER_MEDIA A3
gmtset ANNOT_FONT_SIZE_PRIMARY 36
gmtset ANNOT_FONT_PRIMARY 1
gmtset GRID_PEN_PRIMARY 0.5p,,0.5_3:0
gmtset ANNOT_MIN_SPACING 0
gmtset ANNOT_OFFSET_PRIMARY 0.3c
gmtset ANNOT_OFFSET_SECONDARY 0.2c
gmtset TICK_LENGTH 0.02c
gmtset BASEMAP_FRAME_RGB +0/0/0
echo 0 0 > tmp
echo 2500 0 >>tmp
wl=8
wd=8
cbed=199/199/199
psxy tmp -J -R -K -O -W6,black,- >>$out
awk '{print $1*0.01,$2*0.001}' beda.dat | psxy -J$p -R$d -K -W$wl,black -P -G$cbed -X3.4 -Y3>$out
awk '{print $1*0.01,$2*0.001}' sura.dat | psxy -J$p -R$d -K -W$wl,black -O>>$out
awk '{print $1*0.01,$2*0.001}' basea.dat | psxy -J$p -R$d -K -W$wl,black -O -BSWnea5/a1>>$out
psxy tmp -J -R -K -O -W$wd,black,- >>$out
echo 10 -2.9 32 0 0 "MC" "Distance ( x 100 km)" | pstext -J -R -K -O -N>>$out
echo -2 7 32 90 0 "MC" "Altitude ( x 1000 m)" | pstext -J -R -K -O -N>>$out
awk '{print $1*0.01,$2*0.001}' bedb.dat | psxy -J$p -R$d -K -W$wl,black -Y13 -O -G$cbed>>$out
awk '{print $1*0.01,$2*0.001}' surb.dat | psxy -J$p -R$d -K -W$wl,black -O >>$out
awk '{print $1*0.01,$2*0.001}' baseb.dat | psxy -J$p -R$d -K -W$wl,black -O -BSWne/a1>>$out
psxy tmp -J -R -K -O -W$wd,black,- >>$out
awk '{print $1*0.01,$2*0.001}' bedc.dat | psxy -J$p -R$d -K -W$wl,black -Y13 -O -G$cbed>>$out
awk '{print $1*0.01,$2*0.001}' surc.dat | psxy -J$p -R$d -K -W$wl,black -O >>$out
awk '{print $1*0.01,$2*0.001}' basec.dat | psxy -J$p -R$d -K -W$wl,black -O -BSWne/a1>>$out
psxy tmp -J -R -K -O -W$wd,black,- >>$out
bedb.datなどのデータセットは,grdtrackをつかってつくっています。
#!/bin/sh
gmtset ANNOT_FONT_SIZE_PRIMARY 20
#gfortran shelfbase.f90
#./a.out
bedmap='bedelev.10km'
surface='surface.10km'
base='shelfbase.10km'
conts=7000
out='test.eps'
rm -rf $outfile
startx=370
starty=390
endx=415
endy=30
a=$(($endx-$startx))
echo $a
b=$(($endy-$starty))
echo $b
danmenmap=$startx/$endx/-2500/4500
line=g100a100/g1000a1000
proj=x0.017/0.0037
surface -I1 -R1/561/1/561 -Gtmp1.grd $bedmap
surface -I1 -R1/561/1/561 -Gtmp2.grd $surface
surface -I1 -R1/561/1/561 -Gtmp3.grd $base
xyz2grd -ZTLa -I1 -R1/561/1/561 -G561x561-9999.grd 561x561.-9999
grdmath tmp1.grd 561x561-9999.grd NAN = $bedmap.grd
grdmath tmp2.grd 561x561-9999.grd NAN = $surface.grd
grdmath tmp3.grd 561x561-9999.grd NAN = $base.grd
int=nearneighbor
II=1/0.1
echo 000
gr=200
awk 'BEGIN{for(x=1;x<='$gr';x+=1) print '$startx+$a*x/$gr','$b*x/$gr+$starty'}' > nxny
grdtrack nxny -G$bedmap.grd >bed.dat
grdtrack nxny -G$surface.grd >sur.dat
grdtrack nxny -G$base.grd >base.dat
awk '{print 10*(('$startx'-$1)*('$startx'-$1)+('$starty'-$2)*('$starty'-$2))^0.5,$3}' bed.dat >bed2.dat
awk '{if($3 != "NaN") print 10*(('$startx'-$1)*('$startx'-$1)+('$starty'-$2)*('$starty'-$2))^0.5,$3}' sur.dat >sur2.dat
awk '{if($3 == "NaN") print 10*(('$startx'-$1)*('$startx'-$1)+('$starty'-$2)*('$starty'-$2))^0.5,0.0}' sur.dat >sur3.dat
awk '{if($3*$3>0.01) print 10*(('$startx'-$1)*('$startx'-$1)+('$starty'-$2)*('$starty'-$2))^0.5,$3}' base.dat >base2.dat
head -n 1 sur3.dat >>sur2.dat
d=0/2200/-2500/4500
psxy bed2.dat -J$proj -R$d -K -W8,black >$out
psxy sur2.dat -J$proj -R$d -K -W8,black -O >>$out
psxy base2.dat -J$proj -R$d -K -W8,black -O -BSWnea200/a1000>>$out
grdcontour -Jx0.06 -R1/561/1/561 -C200 -K $surface.grd -X1 -Y1>test2.eps
echo $startx $starty>tmpo
echo $endx $endy>>tmpo
echo 0 0 > tmp
echo 2500 0 >>tmp
psxy -J -R -K -O -W5,red tmpo>>test2.eps
psxy tmp -J -R -K -O -W6,black,- >>$out
メルカトル地図で等距離線を描画
#!/bin/sh
echo 1 1 | psxy -Jx1 -R0/2/0/2 -Ss1000 -N -G244/244/244 -K -X1.5 -Y1.5>test.eps
pscoast -Jq150/0/1:100000000 -R-10/240/-30/86 -Ggray -Ba30g30 -K -O -W1 -N1>> test.eps
lon=140
lat=36
W=5,red
echo $lon $lat 0 2000 2000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 4000 4000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 6000 6000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 8000 8000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 10000 10000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 12000 12000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 14000 14000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 16000 16000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 18000 18000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo $lon $lat 0 20000 20000 | psxy -J -R -SE -K -O -W$W>>test.eps
echo ">" > tmp
echo "140 36" >>tmp
echo "20 62">> tmp
echo "-157 71">> tmp
echo "140 36">> tmp
psxy -J -R -K -O -W3,black -M tmp -G0/160/0>>test.eps
back