Takuya Miyashita

水位　プロット †

空間分布の出力 †

SDxx.tttttttt.grd というファイルに吐き出されるので，とりあえず一通り GMT で描けるようにやってみた．
アスペクト比を計算しているところで awk が必要．
後半の画像の変換には convert (Imagemagick) と ps2pdf と pdfcrop (たぶん texlive-extra-utils) が必要．

計算をした直下のディレクトリで下のスクリプトを実行すればOK．

#!/bin/bash

## check arg
if [ -z "$1" ]; then
    echo "ERROR: empty arg"
    echo "usage: $0 [number XX for (SDXX)]"
    exit 0
else
    NN=`seq -f %02g $1 $1`
fi

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

## makecpt
cpt="tmpsurf$NN.cpt"
gmt makecpt -Cpolar -T-1/1 -D > $cpt

## time interval
dt=`cat tsun.par | grep dt | sed -e 's/dt=//'`

## plot for each
for f in SD$NN.[0-9]*[0-9].grd
do

    ## filename
    outps=${f//\.grd/\.ps}

    itstep=${f//\.grd/}
    itstep=${itstep//SD$NN./}
    time_min=`echo "$itstep $dt" | awk '{printf "%d min", $1*$2/60}'`

    ## plot
    proj="X"$(gmt grdinfo $f -Cn -o0,1,2,3  | awk '{print 10"/"10*($4-$3)/($2-$1)}')
    # echo $proj  # check

    gmt grdimage $f -J$proj -Baf -BneSW+t"$time_min" -R$f -C$cpt -K > $outps
    gmt psscale -C$cpt -Bxa0.5f0.1 -By+lm -DJMR+w5.0/0.3+o1.0/0.0+e -J$proj -R$f -O >> $outps

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

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

done

rm $cpt

水位をコンター＋流速を矢印 †

流速は "-vy" のついたファイルの y 軸方向の出力を逆向きする．
これは支配方程式の緯度方向が colatitude (北極で0，南極で180)表記で，これに沿う方向で流速も出しているためだと思う．

#!/bin/bash

## check arg
if [ -z "$1" ]; then
    echo "ERROR: empty arg"
    echo "usage: $0 [number XX for (SDXX)]" 1>&2
    exit 1
else
    NN=`seq -f %02g $1 $1`
fi

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

## makecpt
cpt="tmpvel$NN.cpt"
gmt makecpt -Cpolar -T-1/1 -D -M --COLOR_NAN=gray90 > $cpt

## time interval
dt=`cat tsun.par | grep dt | sed -e 's/dt=//'`

## plot for each
for f in SD$NN.[0-9]*[0-9].grd
do

    ## filename
    outps=${f//\.grd/\-vxy.ps}
    
    itstep=${f//\.grd/}
    itstep=${itstep//SD$NN./}
    time_min=`echo "$itstep $dt" | awk '{printf "%d min", $1*$2/60}'`

    gu=${f//\.grd/-vx\.grd}
    gv=${f//\.grd/-vy\.grd}
    
    gvtmp=${gv//.grd/-neg.grd}
    gmt grdmath $gv NEG = $gvtmp=cf

    ## params
    proj="X10/`gmt grdinfo $f -Cn -o0,1,2,3  | awk '{print 10*($4-$3)/($2-$1)}'`"
    int=`gmt grdinfo $f -Cn -o6  | awk '{print 40*$1}'`
    XS=`gmt grdinfo $f -Cn -o0,1 | awk '{print $1+0.20*($2-$1)}'`
    YS=`gmt grdinfo $f -Cn -o2,3 | awk '{print $1+0.95*($2-$1)}'`

    vecscale=0.5
    vecatt="+a25+e+p1+h1+gblack+n1"
    fonts=12
    lw=1
    echo "# lon   lat     u1    u2  sig1  sig2  cor  (option) name" > scale$NN.dat
    echo "   $XS  $YS   2.00  0.00  0.00  0.00  0.0  2.0 m/s" >> scale$NN.dat

    ## plot
    gmt grdimage $f -J$proj -Baf -BneSW+t"$time_min" -R$f -C$cpt -K > $outps
    gmt psscale -C$cpt -Bxa0.5f0.1 -By+lm -DJMR+w5.0/0.3+o1.0/0.0+e -J$proj -R$f -K -O >> $outps
    gmt grdvector $gu $gvtmp -W$lw -I$int -Si$vecscale -Q$vecatt -T -J$proj -R$f -K -O >> $outps
    gmt psvelo -A$vecatt -Gblack -Se$vecscale/0/$fonts -W$lw  -J$proj -R$f -O scale$NN.dat >> $outps && rm -f scale$NN.dat # これで正しいスケールになっているかわからん
    rm $gvtmp

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

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

done

rm -f $cpt

アニメーションにまとめる †

上のアニメーションのように，生成した画像から gif を作成するスクリプトを書いてみた．
1番目と2番目の引数に数字を入れると，領域番号の開始と終了を指定できる．
引数がなければ gridfile.dat から自動で拾ってくる．
フレームレートは -r オプションで指定する．

#!/bin/bash
usage_exit() {
        echo "Usage: $0 [-r fps] [num_start] [num_end]" 1>&2
        exit 1
}

## optional arg
while getopts r:h OPT
do
    case $OPT in
        r)  fps=$OPTARG
            ;;
        h)  usage_exit
            ;;
        \?) usage_exit
            ;;
    esac
done
if [ -z "$fps" ]; then
    fps=10
fi
shift $((OPTIND - 1))


if [ $# = 0 ]; then
    SDSTART=1
    if [ -e 'gridfile.dat' ]; then
        SDEND=`cat gridfile.dat | wc -l`
    else
        SDEND=1
    fi
fi
if [ $# = 1 ]; then
    SDSTART=$1
    SDEND=$1
fi
if [ $# = 2 ]; then
    SDSTART=$1
    SDEND=$2
fi



gifdir="_gif"
if [ ! -e $gifdir ]; then
    mkdir $gifdir
fi

#skiprate=`cat tsun.par | grep itmap | sed -e 's/itmap=//'`

## animation
for N in `seq -f %02g $SDSTART $SDEND`
do
    ## directory
    figdir="_surf$N"
    ffmpeg -i $figdir/SD$N.%*.png -filter_complex "fps=$fps,scale=720:-1:flags=lanczos,split[a],palettegen,[a]paletteuse" -an $gifdir/SD$N.gif -y
done