Takuya Miyashita

定点 (tgs) のプロット メモ †

tgs ファイルの成型 †

awk を叩いて直接 tgs ファイルをプロットしても良いけど，ムダが多いので一旦出力データ間引いて出力する．
全ステップを描く必要もなく， 5ステップに1回分に間引いて新たにファイルを吐くとすれば

#!/bin/bash

if [ -z "$1" ]; then
    simdir="."
else
    ## simulation
    simdir="$1"
    if [ ! -e "$simdir" ]; then
        echo "ERROR: not found $simdir"
        exit 1
    fi
fi

## output dir
datdir="$simdir/_dattgs"
if [ ! -e "$datdir" ]; then
    mkdir $datdir
fi

## number of files
ntgs=`ls -1 $simdir/tgs* | wc -l`
echo "found $ntgs files..."

## processing
for f in $simdir/tgs*
do
    fname=$(basename $f)
    outdat="${fname//tgs/station}.dat"
    echo "$fname -> $outdat ..."

    dep=`awk 'NR==1 {print $6}' $f | sed -e 's/depth=//'`

    seaflag=`echo "$dep < 0.0" | bc`
    ## output
    if [ "$seaflag" = 1 ]; then
        awk -v dep=$dep \
	'{if (NR>1 && NR%5==1) {\
	 if ($5+dep > 0) \
              {printf "%12.2f %12.5f %12.5f %12.5f\n", $3, $5+dep, $7, -1*$9}  \
         else {printf "%12.2f %12.5f %12.5f %12.5f\n", $3, 0, $7, -1*$9} \
	} }' $f  > $outdat
    else
        awk '{ if (NR>1 && NR%5==1) {printf "%12.2f %12.5f %12.5f %12.5f\n", $3, $5, $7, -1*$9} }' $f > $outdat
    fi

    mv $outdat $datdir/

done

地点が陸地にある場合は浸水深を出力するように，1行目の水深を読み取って，水位（=平均海水面から高さ）から標高（=負の水深）分を差し引くようにしている．
地形変化を反映させるように設定している場合はこれにその地形の鉛直上昇 or 沈下分を含めないと浸水深にならないので注意．

水位・流速の時系列プロット †

10 地点ずつに分けて出力するスクリプトを書いてみた．
地点数が多く，手早く出力して概要を把握したいときに使う．

#!/bin/bash

var="eta"
#var="vel"

## number of tgs in a single figure
ntgs=10

## plot
if [ "$var" = "eta" ]; then
    prefix="waveform"
    proj="X14/6"
    region="0/365/-1/4"
    xframe=xa30f15g30+l"time (min)"
    yframe=ya1.0f0.5g1.0+l"height (m)"
    frame=neSW
else
    prefix="velocity"
    proj="X14/6"
    region="0/365/-0.2/3.0"
    xframe=xa30f15g30+l"time (min)"
    yframe=ya1.0f0.5g1.0+l"velocity (m/s)"
    frame=neSW
fi

## legend setting
legendfile="tmp_legend.txt"
dx1="0.8"
dx2="1.8"
sym="-"
ms="1.0"
lw="0.5p"


## makecpt
cpt="tmp_tgslc.cpt"
gmt makecpt -Ccategorical -T1/256/1 > $cpt

## directory
figdir="_tgs"
if [ ! -e $figdir ]; then
    mkdir $figdir
fi


## plot for each
cnt=0
for f in tgs*
do
    if [ $cnt = 0 ]; then
        outps="${prefix}_${f}_$ntgs.ps"
        echo $outps
        echo "A $cpt" > $legendfile

	## base plot
	gmt psbasemap -J$proj -R$region -B"$xframe" -B"$yframe" -B"$frame" -K --MAP_GRID_PEN_PRIMARY=thinnest,gray > $outps
    fi

    y=`awk 'NR==1 {print $4}'  $f | sed -e 's/lat=//'`
    x=`awk 'NR==1 {print $5}'  $f | sed -e 's/lon=//'`
    dep=`awk 'NR==1 {print $6}' $f | sed -e 's/depth=//'`


    sID=`awk 'NR==1 {print $3}' $f | sed -e 's/No.=//'`
    lc=`awk -v row=$sID 'NR==row {print $2}' $cpt`

    f_aligned="_dattgs/${f//tgs/station}.dat"

    ## legend file
    #     Spec  dx1   symbol  size  fill  pen     dx2   text"
    echo "S     $dx1  $sym    $ms   -     1p,$lc  $dx2  Station $sID" >> $legendfile

    if [ "$var" = "eta" ]; then
        awk '{printf "%12.5f %12.5f\n", $1/60, $2 }' $f_aligned | \
        gmt psxy -J$proj -R$region -W$lw -C$cpt -Z$sID -O -K >> $outps
    else
        awk '{printf "%12.5f %12.5f\n", $1/60, sqrt($3^2+$4^2) }' $f_aligned | \
        gmt psxy -J$proj -R$region -W$lw -C$cpt -Z$sID -O -K >> $outps
    fi


    cnt=`expr $cnt + 1`
    if [ $cnt = $ntgs ]; then
        ## legend
	gmt pslegend -J -R -DJMR+o0.8/0.0 -F+p0.5+gwhite -O $legendfile >> $outps
	rm -f $legendfile

   　    ## convert
        gmt psconvert -A -Tf $outps # PS -> PDF
　　　　    gmt psconvert -A -TG $outps # PS -> PNG

        ## move
        \mv ${outps//.ps/.pdf} ${outps//.ps/.jpg} $figdir/

	## reset counter
        cnt=0
    fi

done

if [ $cnt != 0 ]; then
        ## legend
	gmt pslegend -J -R -DJMR+o0.8/0.0 -F+p0.5+gwhite -O $legendfile >> $outps
	rm -f $legendfile

   　    ## convert
        gmt psconvert -A -Tf $outps # PS -> PDF
　　　　    gmt psconvert -A -TG $outps # PS -> PNG

        ## move
        \mv ${outps//.ps/.pdf} ${outps//.ps/.jpg} $figdir/
fi

rm -f $cpt