vortex_9696: 7月 2011

2011年7月31日日曜日

番号一覧　線の種類

線種

線種番号	線種
0	(なし)
1	実線
2	長破線
3	短破線
4	長破線+短破線
5	点線
6	一点鎖線
7	二点鎖線

番号一覧マーク、色

番号一覧

色

番号	名前	色
0	background	■■■
1	foreground	■■■
2	red	■■■
3	green	■■■
4	dark blue	■■■
5	light blue	■■■
6	magenta	■■■
7	yellow	■■■
8	orange	■■■
9	purple	■■■
10	yellow/green	■■■
11	medium blue	■■■
12	dark yellow	■■■
13	aqua	■■■
14	dark purple	■■■
15	gray	■■■

マークの種類

番号	マーク
0	(なし)
1	＋
2	○
3	●
4	□
5	■
6	×
7	◇
8	△
9	▲
10	○
11	●

twin cycloneをプロット

~/wkyoshimura2_3/draw4.gs
draw4.gs

*tmp1=0
*tmp=read(data.dat)
*while(tmp1=0)
tmp=read(data.dat)
tmp1=sublin(tmp,1)
tmp2=sublin(tmp,2)
nlon=subwrd(tmp2,1)
nlat=subwrd(tmp2,2)
slon=subwrd(tmp2,3)
slat=subwrd(tmp2,4)
say tmp2
sym = 3
siz = 0.15
'q w2xy 'nlon' 'nlat''
x1 = subwrd(result,3)
y1 = subwrd(result,6)
'draw mark 'sym' 'x1' 'y1' 'siz
'q w2xy 'slon' 'slat''
x2 = subwrd(result,3)
y2 = subwrd(result,6)
'draw mark 'sym' 'x2' 'y2' 'siz
'draw line 'x1' 'y1' 'x2' 'y2''

2011年7月27日水曜日

gradsのラベル位置、文字位置、大きさ　gradsの文字非表示

gradsのラベル位置
ga>run cbarn.gs

ラベル文字位置、大きさ

set xlopts 色番号 [太さ [大きさ]]
 set ylopts 色番号 [太さ [大きさ]]

ラベルに用いるフォントの色、太さ、大きさを設定します。デフォルト値は色番号=1、太さ=4、大きさ=0.12です。xloptsはX軸(横軸)、yloptsはY軸(縦軸)のラベルを設定します。

(例)x軸のラベル：色1、太さ2、大きさ0.2、を設定する
```
ga-> set xlopts 1 2 0.2
```
(例)y軸のラベル：色1、太さ1、大きさ0.2、を設定する
```
ga-> set ylopts 1 1 0.2
```

x軸のラベルを描く場所を指定する

set xlpos 位置 側

「位置」にはx軸からのずれ、「側」にはt(上側)、b(下側)のどちらかを指定します。

(例)x軸のラベルを、下側の通常より0.5下に描画する
```
ga-> set xlpos -0.5 b
```

y軸のラベルを描く場所を指定する

set ylpos 位置 側

「位置」にはy軸からのずれ、「側」にはr(右側)、l(左側)のどちらかを指定します。

(例)y軸のラベルを、左側の通常より0.5左に描画する
```
ga-> set ylpos -0.5 l
```

gradsの文字非表示

GrADSのロゴ（左下）と日付（右下）を表示する／しない

set grads on/off

日付（右下）を表示する／しない

set timelab on/off

2011年7月26日火曜日

NHMの計算結果からpminとvmaxそして、それぞれの位置を求めるシェルスクリプト

NHMの計算結果からpminとvmaxそして、それぞれの位置を求めるシェルスクリプト

find_n_v_pmin.sh

#! /bin/tcsh -f

#####################################
filename=find_n_v_pmin

#type yyyymmddhhhh start and end time
: ${start:=199201020000}
: ${end:=199201081000}
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`
hour=`echo ${start} |cut -c 9-10`
h=00

echo start time ${year} ${month} ${day} ${hour} ${h}
echo end   time ${end}

: ${time:=${start}}

#timestep
timestep=6

time_step=0

i=0
while [ ${time} -le ${end} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}

##decide m
if [ ${month} -eq 01 ]; then
    m=JAN
elif [ ${month} -eq 02 ]; then
    m=FEB
elif [ ${month} -eq 03 ]; then
    m=MAR
elif [ ${month} -eq 04 ]; then
    m=APR
elif [ ${month} -eq 05 ]; then
    m=MAY
elif [ ${month} -eq 06 ]; then
    m=JUN
elif [ ${month} -eq 07 ]; then
    m=JUL
elif [ ${month} -eq 08 ]; then
    m=AUG
elif [ ${month} -eq 09 ]; then
    m=SEP
elif [ ${month} -eq 10 ]; then
    m=OCT
elif [ ${month} -eq 11 ]; then
    m=NOV
elif [ ${month} -eq 12 ]; then
    m=DEC
fi

##decide m

#echo $hour
time=${year}${month}${day}${hour}${h}
#echo yoshiyoshi
echo time $time
echo month $m

    t=`expr ${day} \* 4 - 3`
    t=`expr ${hour} / 6 + ${t}`

: ${prefix:=anl_p}    #want data

    cat > data.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time}.ctl'

'set gxout fwrite'
'set undef dfile'
'set fwrite ${filename}.grd'
*'set t ${t}'
'set x 1 400'
'set y 1 300'
'd psea'
'disable fwrite'
'c'
'set gxout fwrite'
'set undef dfile'
'set fwrite v_${filename}.grd'
*'set t ${t}'
'set x 1 400'
'set y 1 300'
'd sqrt(u*u+v*v)'
'disable fwrite'
'c'

'quit'

#157
EOF
    grads -blc data.gs
########################################
####program read_data###################
########################################

if [ ${time} -eq ${start} ]; then
##start file content
cat <<EOF > start.f90
program start
implicit none
open (unit=7, file='${filename}.txt', form='FORMATTED')
write(7,*) '${filename} NH'
write(7,*) 'time_step mini_psea xgrid ygrid lon lat vmax vlon vlat'
write(*,*) 'time_step mini_psea xgrid ygrid lon lat vmax vlon vlat'
close(7)
end program start
EOF
##output data.txt by f90
f90 -o start.out start.f90
./start.out
##start
fi

####find min psea
cat <<EOF > ${filename}.f90
program find_pmin
implicit none
integer, parameter :: xgrid=400, ygrid=300, rgrid=10
integer :: n, i, j, kx, ky, is, ie, js, je
real(kind=4), dimension(xgrid,ygrid) :: p_d
real(kind=4), dimension(xgrid,ygrid) :: v_d
real(kind=4) :: p, pp, lon, lat, v, vv, vlon ,vlat
open (unit=17,file='${filename}.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=17,rec=1) p_d
open (unit=18,file='v_${filename}.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=18,rec=1) v_d

open (unit=7, file='${filename}.txt', position='append', form='FORMATTED')

p=2000
do i = 139, 277 !155~180
   do j = 150, ygrid !150 @NH
    pp= p_d(i-1,j) + p_d(i+1,j) + p_d(i,j-1) + p_d(i,j-1)
    pp= pp / 4
      if (p_d(i,j) < pp) then
   if (p_d(i,j) < p) then
        kx=i
        ky=j
        p=p_d(i,j)
!         write(*,*) n, ${time_step}, i, j, p
!         write(7,*) n, ${time_step}, i, j, p
   endif
      endif
   n=n+1
   enddo
enddo

lon=kx*0.18+130
lat=ky*0.18-30

is=kx-10
js=ky-10
ie=kx+10
je=ky+10

do i = is, ie
    do j = js, je
    vv= v_d(i-1,j) + v_d(i+1,j) + v_d(i,j-1) + v_d(i,j-1)
    vv= vv / 4
      if (v_d(i,j) > vv) then
   if (v_d(i,j) > v) then
        kx=i
        ky=j
        v=v_d(i,j)
   endif
      endif
    enddo
enddo

vlon=kx*0.18+130
vlat=ky*0.18-30

write(*,*) ${time_step}, p, kx, ky, lon, lat, v, vlon, vlat
write(7,*) ${time_step}, p, kx, ky, lon, lat, v, vlon, vlat

close(7)
close(17)
close(18)

end program find_pmin
EOF

##output data.txt by f90
f90 -o ${filename}.out ${filename}.f90
./${filename}.out

#echo ${day}
#echo "yoshiyoshi"
i=`expr $i + 1`
#echo $i

#hour=`expr $hour + 1`
hour=`expr $hour + $timestep`
hour=`echo ${hour} | awk '{printf("%02d",$1)}'`

if [ $hour -eq 24 ]; then
    hour=00
#    day=`expr $day + 1`
       day=`echo ${day} | awk '{printf("%02d",$1+1)}'`

   if [ $day -ge 32 ];then
        day=01
#        month=`expr $month + 1`
       month=`echo ${month} | awk '{printf("%02d",$1+1)}'`

       if [ $month -ge 13 ];then
            month=01
            year=`expr $year + 1`

       fi
   fi
fi
time=${year}${month}${day}${hour}${h}

done

NHM計算結果のtrackとJTWCのtrackデータを描く

ＪＴＷＣのbest_trackをgradsで表示する

台風再現実験において、best_trackデータをgradsで再現する　gsスクリプトを作成
~/yoshi/best_trackで作成する

準備１
ＪＴＷＣのサイトから、描きたいbest_trackデータを持ってくる
datファイルとしてファイルを作成する。
以下に中身を示す
bwp011992.datの中身

WP, 01, 1992010306, , BEST, 0, 33N, 1757E, 20
・・・・

その次に、
・draw_track3_1.gsでreadするデータのファイル名を変える（今、作ったファイルを参照させる）
・whileループの繰り返し数を参照するデータの数に合わせる

ここまでは前回と同じ。
後述するfind_n_ｖ_pmin.shによって作られた
find_n_ｖ_pmin.txtに手を加える

１．txtファイルからdatファイルに変える（ファイルの語尾をかえればOK ）
２．ファイル名が長すぎるので、find_n.datなどにする
３．読み込みたいデータを先頭にする（余計なものは消す）
30 1008.134338 236 197 171.9400024 3.840001345 11.22286606 172.4799957 5.460001469
36 1007.647217 214 177 170.3200073 1.320001245 6.509649754 168.5200043 1.860001326
42 1007.682495 234 199 170.6800079 4.020001411 10.66093826 172.1199951 5.820001602
48 1007.458923 226 195 169.2400055 3.300001383 10.83714390 170.6800079 5.100001335
54 1007.141785 216 190 170.1399994 2.400001287 8.402082443 168.8800049 4.200001240

・・・・・・・・・・・・・・
その次に、
・draw_track_mhm1_1.gsでreadするデータのファイル名を変える（今、作ったファイルを参照させる）
・whileループの繰り返し数を参照するデータの数に合わせる

準備完了
まずgradsを立ち上げる
$grads
ga>open gmap.ctl
ga>blankmap.gs
これで白地図が完成
（白地図が作れないときは、適当に地形データを作成、開いたものに記述する）

ga>draw_track3_1.gs
台風の強度別（Vmax）に色を変えてプロット

つぎにNHMの計算結果をプロット
ga>draw_track_nhm1_1.gs

2011年7月25日月曜日

fortran ファイルへの書き足し

open (unit=7, file='psea.txt', position='append', form='FORMATTED')

position='append'を書けばOK
ふつうは、status='UNKNOWN'
でよいはずなのだが・・・・

それか、linux系だから

a.out>>output
として、outputに書かせる方法もあり！

#! /bin/tcsh -f

##make pmin.txt
#type yyyymmddhhhh start and end time
: ${start:=199201020000}
: ${end:=199201050000}
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`
hour=`echo ${start} |cut -c 9-10`
h=00

echo start time ${year} ${month} ${day} ${hour} ${h}
echo end   time ${end}

: ${time:=${start}}

#timestep
timestep=12

time_step=0

i=0
while [ ${time} -le ${end} ]
do
time_step=`expr ${timestep} \* ${i}`
echo time_step ${time_step}

##decide m
if [ ${month} -eq 01 ]; then
    m=JAN
elif [ ${month} -eq 02 ]; then
    m=FEB
elif [ ${month} -eq 03 ]; then
    m=MAR
elif [ ${month} -eq 04 ]; then
    m=APR
elif [ ${month} -eq 05 ]; then
    m=MAY
elif [ ${month} -eq 06 ]; then
    m=JUN
elif [ ${month} -eq 07 ]; then
    m=JUL
elif [ ${month} -eq 08 ]; then
    m=AUG
elif [ ${month} -eq 09 ]; then
    m=SEP
elif [ ${month} -eq 10 ]; then
    m=OCT
elif [ ${month} -eq 11 ]; then
    m=NOV
elif [ ${month} -eq 12 ]; then
    m=DEC
fi

##decide m

#echo $hour
time=${year}${month}${day}${hour}${h}
#echo yoshiyoshi
echo time $time
echo month $m

    t=`expr ${day} \* 4 - 3`
    t=`expr ${hour} / 6 + ${t}`

: ${prefix:=anl_p}    #want data

    cat > data.gs <<EOF
'reinit'
'open _NHMMRPPFCSVSTD1_${time}.ctl'

'set gxout fwrite'
'set undef dfile'
'set fwrite p1.grd'
*'set t ${t}'
'set x 1 400'
'set y 1 300'
'd psea'
'disable fwrite'
'c'
*'set gxout fwrite'
*'set undef dfile'
*'set fwrite u1.grd'
**'set t ${t}'
*'set x 1 288'
*'set y 1 145'
*'d ugrdsfc'
*'disable fwrite'
*'c'
*'set gxout fwrite'
*'set undef dfile'
*'set fwrite v1.grd'
*'set t ${t}'
**'set x 1 288'
*'set y 1 145'
*'d vgrdsfc'
*'disable fwrite'
*'c'
*'set gxout fwrite'
*'set undef dfile'
*'set fwrite rv.grd'
*'set t ${t}'
*'set x 1 288'
*'set y 1 145'
*'a=hcurl(ugrdprs,vgrdprs)'
*'d a'
*'disable fwrite'
*'c'

'quit'

#157
EOF
    grads -blc data.gs

####program read_data###################

if [ ${time} -eq ${start} ]; then
##start file content
cat <<EOF > start.f90
program start
implicit none
open (unit=7, file='psea.txt', form='FORMATTED')
write(7,*) 'time_step mini_psea xgrid ygrid lon lat'
write(*,*) 'time_step mini_psea xgrid ygrid lon lat'
close(7)
end program start
EOF
##output data.txt by f90
f90 -o start.out start.f90
./start.out
##start
fi

####find min psea
cat <<EOF > f_find_pmin.f90
program find_pmin
implicit none
integer, parameter :: xgrid=400, ygrid=300
integer :: n, i, j, kx, ky
real(kind=4), dimension(xgrid,ygrid) :: p_d
real(kind=4) :: p, lon, lat
open (unit=17,file='p1.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=17,rec=1) p_d
open (unit=7, file='psea.txt', position='append', form='FORMATTED')

p=2000
do i = 1, xgrid
   do j = 1, ygrid
      if (p_d(i,j) < p) then
         kx=i
         ky=j
        p=p_d(i,j)
         write(*,*) n, ${time_step}, i, j, p
!         write(7,*) n, ${time_step}, i, j, p
      endif
   n=n+1
   enddo
enddo

lon=kx*20*0.18+130
lat=ky*20*0.18-30

write(*,*) ${time_step}, p, kx, ky, lon, lat
write(7,*) ${time_step}, p, kx, ky, lon, lat

close(7)
close(17)

end program find_pmin
EOF

##output data.txt by f90
f90 -o output.out f_find_pmin.f90
./output.out

#echo ${day}
#echo "yoshiyoshi"
i=`expr $i + 1`
#echo $i

#hour=`expr $hour + 1`
hour=`expr $hour + $timestep`
hour=`echo ${hour} | awk '{printf("%02d",$1)}'`

if [ $hour -eq 24 ]; then
    hour=00
#    day=`expr $day + 1`
       day=`echo ${day} | awk '{printf("%02d",$1+1)}'`

   if [ $day -ge 32 ];then
        day=01
#        month=`expr $month + 1`
       month=`echo ${month} | awk '{printf("%02d",$1+1)}'`

       if [ $month -ge 13 ];then
            month=01
            year=`expr $year + 1`

       fi
   fi
fi
time=${year}${month}${day}${hour}${h}

done

ファイルの行数をカウントするプログラム例

program count_lines
  implicit none
  character(1000) filename
  integer counter
  
  print *, 'Enter file name to count lines: '
  read '(A)', filename
  open(11,file=filename,status='old')
  counter = 0
  do
    read(11,'()',end=100)
    counter = counter + 1
  end do
100 close(11)
  print *,'Number of lines is',counter
end program count_lines

fortran,gnuplt 参考ページ

全般

2009年度　地球惑星物理学演習
http://www-geoph.eps.s.u-tokyo.ac.jp/ta/resume2009/

fortran
東大高木先生
http://www-aos.eps.s.u-tokyo.ac.jp/~takagi/lecture/f90-enshu/index.html

gnuplot
http://t16web.lanl.gov/Kawano/gnuplot/

2011年7月24日日曜日

gradsからデータを読み込み、pminを探すf90

gradsからデータを読み込み、pminを探すf90

だたし、gxoutを使ってgradsにおいてp1.grdを作ってから走らせること。
@yoshi/analysis/f_find_pmin.f90

program find_pmin
implicit none
integer, parameter :: xgrid=400, ygrid=300
integer :: n, i, j, kx, ky
real(kind=4), dimension(xgrid,ygrid) :: p_d
real(kind=4) :: p
open (unit=17,file='p1.grd',form='unformatted',access='direct',recl=xgrid*ygrid*4)
read (unit=17,rec=1) p_d

open (unit=7, file='psea.txt', form='FORMATTED')
write(7,*) 'write mini psea xgrid ygrid lon lat'

p=2000
do i = 1, xgrid
   do j = 1, ygrid
      if (p_d(i,j) < p) then
         kx=i
         ky=j
        p=p_d(i,j)
         write(*,*) i, j, n
         write(7,*) n, 'min psea', p, 'kx', kx, 'ky', ky
      endif
   n=n+1
   enddo
enddo

write(7,*) n, ' out put   min psea', p, 'kx', kx, 'ky', ky
write(7,*) 'yoshiyoshi'

close(7)
close(17)

end program find_pmin

2011年7月5日火曜日

JRA25 SST

JRA25 から　SST
mk_sst_jra25.sh @yoshi/tool/jra25
#! /bin/tcsh -f

#type yyyymmddhhhh start and end time
: ${start:=199201010000}
: ${end:=199201160000}
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`
hour=`echo ${start} |cut -c 9-10`
h=00
: ${prefix:=fcst_phy2m25}
: ${suffix:=gr}

echo start time ${year} ${month} ${day} ${hour} ${h}
echo end   time ${end}

: ${time:=${start}}

#timestep
timestep=12
cp ~/tool/color.gs .

i=1
while [ ${time} -le ${end} ]
do

##decide m
if [ ${month} -eq 01 ]; then
    m=JAN
elif [ ${month} -eq 02 ]; then
    m=FEB
elif [ ${month} -eq 03 ]; then
    m=MAR
elif [ ${month} -eq 04 ]; then
    m=APR
elif [ ${month} -eq 05 ]; then
    m=MAY
elif [ ${month} -eq 06 ]; then
    m=JUN
elif [ ${month} -eq 07 ]; then
    m=JUL
elif [ ${month} -eq 08 ]; then
    m=AUG
elif [ ${month} -eq 09 ]; then
    m=SEP
elif [ ${month} -eq 10 ]; then
    m=OCT
elif [ ${month} -eq 11 ]; then
    m=NOV
elif [ ${month} -eq 12 ]; then
    m=DEC
fi

time=${year}${month}${day}${hour}
tim=${hour}z${day}${m}${year}
echo time $time

##change 3 time
    t=`expr ${day} \* 4 - 3`
    t=`expr ${hour} / 6 + ${t}`

#make press wind rv image
    cat > sst1.gs <<EOF
'reinit'
'open ${prefix}.ctl'

'set t ${t}'

#
'set lon 100 200'
'set lat -25 25'

'set gxout shaded'
'set display color white'
'c'
'color 0 33 3 -kind white->yellow->orange->red'
'd wtmpsfc-273.15'
'cbar'
'set gxout contour'
'd wtmpsfc-273.15'

'q dims'
temp=sublin(result,5)
tim=subwrd(temp,6)

'draw title SST[deg C] time='tim''

'enable print t='tim'.gmf'
'print'
'disable print'
'!gxps -c -i t='tim'.gmf -o sst_t='tim'.ps'
'!rm t='tim'.gmf'
'!ps2pdf sst_t='tim'.ps ${i}_sst_t='tim'.pdf'
'!rm sst_t='tim'.ps'

'quit'
EOF
    grads -blc sst1.gs

#echo "yoshiyoshi"
i=`expr $i + 1`
#echo $i

#### next time step
#hour=`expr $hour + 1`
hour=`expr $hour + $timestep`
hour=`echo ${hour} | awk '{printf("%02d",$1)}'`

if [ $hour -eq 24 ]; then
    hour=00
#    day=`expr $day + 1`
       day=`echo ${day} | awk '{printf("%02d",$1+1)}'`

   if [ $day -ge 32 ];then
        day=01
#        month=`expr $month + 1`
       month=`echo ${month} | awk '{printf("%02d",$1+1)}'`

       if [ $month -ge 13 ];then
            month=01
            year=`expr $year + 1`

       fi
   fi
fi

time=${year}${month}${day}${hour}${h}

done

2011年7月4日月曜日

JAR25からデータを取ってくるシェルスクリプト

JAR25からデータを取ってくるシェルスクリプト

getby6by.sh

#!/bin/sh
######################################################################
# initialization
######################################################################
set -eu

: ${user:=jra01479}
: ${pass:=*******}
: ${option:="-N -r -nd -np -w 3 --level=0"}
: ${URL:=http://ds.data.jma.go.jp/gmd/jra/download/data/GribFinal} #when by 6hour
#: ${URL:=http://ds.data.jma.go.jp/gmd/jra/download/data/MonthFinal} #when by day or month
: ${prefix:=anl_p}    #want data
: ${suffix:=gr}
: ${start:=20110101} #start time
: ${end:=20110229}    #finish time

: ${WGET:=wget}

######################################################################
# get JRA-25 data
######################################################################
year=`echo ${start} | cut -c 1-4`
month=`echo ${start} | cut -c 5-6`
day=`echo ${start} | cut -c 7-8`   #when by 6hour
hour=00                            #when by 6hour

#get control file
mkdir ${start}${prefix}
cd ${start}${prefix}

${WGET} \
    ${option} \
    ${user:+--http-user=${user}} \
    ${pass:+--http-passwd=${pass}} \
    ${URL}/${prefix}/${year}${month}/${prefix}.ctl
${WGET} \
    ${option} \
    ${user:+--http-user=${user}} \
    ${pass:+--http-passwd=${pass}} \
    ${URL}/${prefix}/${year}${month}/${prefix}.idx #when by 6hour

#download
while [ ${year}${month}${day} -le ${end} ]
do
    ${WGET} \
    ${option} \
    ${user:+--http-user=${user}} \
    ${pass:+--http-passwd=${pass}} \
    ${URL}/${prefix}/${year}${month}/${prefix}.${year}${month}${day}${hour}
        #${URL}/${prefix}/${prefix}.${year}${month}.${suffix} #when by day or month

#tomorrow
    if [ ${hour} -eq 18 ]; then
    hour=00
    if [ ${month} -eq 2 ]; then
        if [ `expr ${year} % 4` -eq 0 -a \
        `expr ${year} % 1000` -ne 0]; then
        if [ ${day} -eq 30 ]; then
            day=01
            month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
            rm -f ${prefix}.ctl
            rm -f ${prefix}.idx
            ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.ctl
            ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.idx
        else
            day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
        fi
        else
        if [ ${day} -eq 29 ]; then
            day=01
            month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
            rm -f ${prefix}.ctl
            rm -f ${prefix}.idx
            ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.ctl
            ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.idx
        else
            day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
        fi
        fi
    elif [ ${month} -eq 04 -o ${month} -eq 06 -o \
        ${month} -eq 09 -o ${month} -eq 11 ]; then
        if [ ${day} -eq 31 ]; then
        day=01
        month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
        rm -f ${prefix}.ctl
        rm -f ${prefix}.idx
        ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.ctl
        ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.idx
        else
        day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
        fi
    else
        if [ ${day} -eq 32 ]; then
        day=01
        if [ ${month} -eq 12 ]; then
            month=01
            year=`expr ${year} + 1`
        else
            month=`echo ${month} | awk '{printf("%02d",$1+1)}'`
        fi
        rm -f ${prefix}.ctl
        rm -f ${prefix}.idx
        ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.ctl
        ${WGET} \
            ${option} \
            ${user:+--http-user=${user}} \
            ${pass:+--http-passwd=${pass}} \
            ${URL}/${prefix}/${year}${month}/${prefix}.idx
        else
        day=`echo ${day} | awk '{printf("%02d",$1+1)}'`
        fi
    fi
    else
    hour=`echo ${hour} | awk '{printf("%02d",$1+6)}'`
    fi
done

exit 0

2011年7月3日日曜日

gradsの図をpdfにすぐ変換するgsスクリプト

gradsの図をpdfにすぐ変換するgsスクリプト
mkpdf.gs @~/yoshi/best_track

*'make .ps'
'enable print pspdf.gmf'
'print'
'disable print'
'!gxps -c -i pspdf.gmf -o pspdf.ps'
'!rm pspdf.gmf'

'!ps2pdf pspdf.ps pspdf.pdf'

ＪＴＷＣのbest_trackをgradsで表示する

ＪＴＷＣのbest_trackをgradsで表示する

台風再現実験において、best_trackデータをgradsで再現する　gsスクリプトを作成
~/yoshi/best_trackで作成する

準備１
ＪＴＷＣのサイトから、描きたいbest_trackデータを持ってくる
datファイルとしてファイルを作成する。
以下に中身を示す
bwp011992.datの中身

WP, 01, 1992010306,   , BEST,   0, 33N, 1757E, 20
WP, 01, 1992010312,   , BEST,   0, 32N, 1769E, 20
WP, 01, 1992010318,   , BEST,   0, 35N, 1778E, 20
WP, 01, 1992010400,   , BEST,   0, 40N, 1787E, 20
WP, 01, 1992010406,   , BEST,   0, 46N, 1792E, 20
WP, 01, 1992010412,   , BEST,   0, 50N, 1791E, 25
WP, 01, 1992010418,   , BEST,   0, 53N, 1788E, 25
WP, 01, 1992010500,   , BEST,   0, 56N, 1784E, 25
WP, 01, 1992010506,   , BEST,   0, 59N, 1777E, 25
WP, 01, 1992010512,   , BEST,   0, 60N, 1768E, 30
WP, 01, 1992010518,   , BEST,   0, 61N, 1760E, 30
WP, 01, 1992010600,   , BEST,   0, 62N, 1752E, 35
WP, 01, 1992010606,   , BEST,   0, 62N, 1743E, 40
WP, 01, 1992010612,   , BEST,   0, 62N, 1735E, 45
WP, 01, 1992010618,   , BEST,   0, 60N, 1727E, 50
WP, 01, 1992010700,   , BEST,   0, 59N, 1719E, 50
WP, 01, 1992010706,   , BEST,   0, 58N, 1710E, 55
WP, 01, 1992010712,   , BEST,   0, 58N, 1701E, 60
WP, 01, 1992010718,   , BEST,   0, 58N, 1690E, 65
WP, 01, 1992010800,   , BEST,   0, 59N, 1678E, 70
WP, 01, 1992010806,   , BEST,   0, 60N, 1665E, 70
WP, 01, 1992010812,   , BEST,   0, 60N, 1651E, 70
WP, 01, 1992010818,   , BEST,   0, 60N, 1636E, 70
WP, 01, 1992010900,   , BEST,   0, 62N, 1621E, 70
WP, 01, 1992010906,   , BEST,   0, 65N, 1606E, 65
WP, 01, 1992010912,   , BEST,   0, 69N, 1593E, 65
WP, 01, 1992010918,   , BEST,   0, 73N, 1580E, 60
WP, 01, 1992011000,   , BEST,   0, 77N, 1567E, 60
WP, 01, 1992011006,   , BEST,   0, 83N, 1556E, 55
WP, 01, 1992011012,   , BEST,   0, 89N, 1545E, 50
WP, 01, 1992011018,   , BEST,   0, 93N, 1534E, 45
WP, 01, 1992011100,   , BEST,   0, 95N, 1523E, 45
WP, 01, 1992011106,   , BEST,   0, 96N, 1512E, 40
WP, 01, 1992011112,   , BEST,   0, 97N, 1501E, 40
WP, 01, 1992011118,   , BEST,   0, 98N, 1490E, 35
WP, 01, 1992011200,   , BEST,   0, 99N, 1479E, 35
WP, 01, 1992011206,   , BEST,   0, 101N, 1469E, 35
WP, 01, 1992011212,   , BEST,   0, 104N, 1459E, 35
WP, 01, 1992011218,   , BEST,   0, 109N, 1450E, 35
WP, 01, 1992011300,   , BEST,   0, 116N, 1442E, 30
WP, 01, 1992011306,   , BEST,   0, 122N, 1434E, 30
WP, 01, 1992011312,   , BEST,   0, 128N, 1426E, 30
WP, 01, 1992011318,   , BEST,   0, 135N, 1420E, 30
WP, 01, 1992011400,   , BEST,   0, 143N, 1413E, 30
WP, 01, 1992011406,   , BEST,   0, 152N, 1409E, 30
WP, 01, 1992011412,   , BEST,   0, 164N, 1408E, 30
WP, 01, 1992011418,   , BEST,   0, 179N, 1412E, 30
WP, 01, 1992011500,   , BEST,   0, 202N, 1428E, 30
WP, 01, 1992011506,   , BEST,   0, 230N, 1470E, 25

その次に、
・draw_track3_1.gsでreadするデータのファイル名を変える（今、作ったファイルを参照させる）
・whileループの繰り返し数を参照するデータの数に合わせる

準備完了
まずgradsを立ち上げる
$grads
ga>open gmap.ctl
ga>blankmap.gs
これで白地図が完成
（白地図が作れないときは、適当に地形データを作成、開いたものに記述する）

ga>draw_track3_1.gs
台風の強度別（Vmax）に色を変えてプロット
*vmax
vmax=subwrd(tmp2,9)
say vmax_ vmax
if (vmax>64)
say T vmax red
'set line 2'
endif
if (vmax<64 & vmax>48)
say STS vmax orange
'set line 8'
endif
if (vmax<48 & vmax>34)
say TS vmax green
'set line 3'
endif
if (vmax<34)
say TD vmax blue
'set line 4 '
endif

色を変えないバージョンはdraw_track2_1.gs
ちなみに、学会で使ったtwin TCsの表示させるプログラムはdraw4.gs

以下にdraw_track3_1.gsのプログラムを示す

'set display color white'
i=1
while (i<=49)
*set lon and lat
say
say n_hemisphere
tmp=read(bwp011992.dat)
tmp1=sublin(tmp,1)
tmp2=sublin(tmp,2)
nlon=subwrd(tmp2,8)
nlat=subwrd(tmp2,7)
say tmp2
say nlon
say nlat
nnlon=substr(nlon,1,4)
nnlon=nnlon*0.1
say nnlon_is_ nnlon
str_len=strlen(nlat)
say lat_of_length_ str_len
if (str_len<5)
nnlat=substr(nlat,1,2)
nnlat=nnlat*0.1
say nnlat_is_ nnlat
else
nnlat=substr(nlat,1,3)
nnlat=nnlat*0.1
say nnlat_is_ nnlat
endif
say
say lon_and_lat_is
say nnlon__ nnlon
say nnlat__ nnlat
*vmax
vmax=subwrd(tmp2,9)
say vmax_ vmax
if (vmax>64)
say T vmax red
'set line 2'
endif
if (vmax<64 & vmax>48)
say STS vmax orange
'set line 8'
endif
if (vmax<48 & vmax>34)
say TS vmax green
'set line 3'
endif
if (vmax<34)
say TD vmax blue
'set line 4 '
endif
sym = 3
siz = 0.15
'q w2xy 'nnlon' 'nnlat''
x = subwrd(result,3)
y = subwrd(result,6)
'draw mark 'sym' 'x' 'y' 'siz
i=i+1
endwhile

i=1
while (i<=44)
say
say s_hemisphere
tmp=read(bsh111992.dat)
tmp1=sublin(tmp,1)
tmp2=sublin(tmp,2)
slon=subwrd(tmp2,8)
slat=subwrd(tmp2,7)
say tmp2
say slon
say slat
sslon=substr(slon,1,4)
sslon=sslon*0.1
say sslon_is_ sslon
str_len=strlen(slat)
say lat_of_length_ str_len
if (str_len<5)
sslat=substr(slat,1,2)
sslat=-sslat*0.1
say sslat_is_ sslat
else
sslat=substr(slat,1,3)
sslat=-sslat*0.1
say sslat_is_ sslat
endif
say
say lon_and_lat_is
say sslon__ sslon
say sslat__ sslat
*vmax
vmax=subwrd(tmp2,9)
say vmax_ vmax
if (vmax>64)
say T vmax red
'set line 2'
endif
if (vmax<64 & vmax>48)
say STS vmax orange
'set line 8'
endif
if (vmax<48 & vmax>34)
say TS vmax green
'set line 3'
endif
if (vmax<34)
say TD vmax blue
'set line 4 '
endif
sym = 3
siz = 0.15
'q w2xy 'sslon' 'sslat''
x = subwrd(result,3)
y = subwrd(result,6)
'draw mark 'sym' 'x' 'y' 'siz
i=i+1
endwhile

draw4.gsは以下に示す
whileがこのときは使えなかったので、繰り返しになってるｗ

*tmp1=0
*tmp=read(data.dat)
*while(tmp1=0)
tmp=read(data.dat)
tmp1=sublin(tmp,1)
tmp2=sublin(tmp,2)
nlon=subwrd(tmp2,1)
nlat=subwrd(tmp2,2)
slon=subwrd(tmp2,3)
slat=subwrd(tmp2,4)
say tmp2
sym = 3
siz = 0.15
'q w2xy 'nlon' 'nlat''
x1 = subwrd(result,3)
y1 = subwrd(result,6)
'draw mark 'sym' 'x1' 'y1' 'siz
'q w2xy 'slon' 'slat''
x2 = subwrd(result,3)
y2 = subwrd(result,6)
'draw mark 'sym' 'x2' 'y2' 'siz
'draw line 'x1' 'y1' 'x2' 'y2''

gmap.ctl中身
DSET    ^gmap.dat
TITLE   MAP DRAW
UNDEF   1.E36
XDEF          2 LINEAR 100.000   80.00000
YDEF          2 LINEAR -30.000   60.00000
ZDEF          1 LINEAR     .000    1.00000
TDEF    1 LINEAR 00JAN1997 60MN
VARS 1
PT 1 99 PT
ENDVARS

blankmap.gsの中身
'c'
'set vpage 0 8.5 0 5.5'
'set font 1'

'set xlint 10'
'set ylint 10'
'set xlopts 1 5 0.13'
'set ylopts 1 5 0.13'
'set grads off'
'set mpdset mres'
'set map 1 1 5'

'set cmin 1'
'd pt'

'set vpage off'

2011年7月1日金曜日

もじもじ君　ｇｓスクリプト　ｓｔｒ　演算子

ｇｓスクリプト　ｓｔｒ　　演算子

文字の扱い、数字の扱い

文字列を取り出す

ret = substr(string, start, length)

string	文字列
start	整数値（開始位置）
length	整数値（長さ）
ret	指定した文字列

stringが短すぎる場合、NULL（""）が返ります。

行を指定して文字列を取り出す

ret = sublin(string, n)

string	文字列
n	整数値
ret	n行目の文字列

stringの行数がnより小さい場合、NULL（""）が返ります。

空白で区切られた文字列から、n番目の文字列を取り出す

ret = subwrd(string, n)

string	文字列
n	整数値
ret	n番目の文字列

stringが短すぎる場合、NULL（""）が返ります。

空白で区切られた文字列から、n番目の文字列の位置を取得する

ret = wrdpos(string, int)

string	文字列
int	整数値
ret	n番目の文字列の開始位置

テキストファイルの読み込み

ret = read(filename)

演算子

一覧表

優先順位	演算子	意味
1	-	マイナス符号(単項演算子)
1	!	否定(単項演算子)
2	/	除算
2	*	乗算
3	+	加算
3	-	減算
4	%	連結
5	=	等しい
	!=	等しくない
	>	大なり
	>=	大なり又は等しい
	<	小なり
	<=	小なり又は等しい
6	&	論理AND
7	\|	論理OR

(例)a=1またはb=1のときYesと表示する

if(a=1 | b=1)
    say 'Yes'
 else
    say 'No'
 endif

参照　東北大ＨＰ
http://wind.geophys.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGrADS%2FGrADS%A5%B9%A5%AF%A5%EA%A5%D7%A5%C8%A4%CETips

gradsでbar バーの位置を変える

gradsでbar バーの位置を変える

set gxout shaded(p.60)
図にハッチをかける。
入力例:ga>set gxout shaded
!!:ハッチの色や間隔は何も指定しないとgradsが勝手に決めてしまいます。gradsは何段階にも分けてハッチをかけますが、白黒の場合、実際に認識できるのはせいぜい３～４種類です。それでハッチの配色や間隔を指定するコマンドがあります。これはset ccols とset clevs です。
!:gxout shadedを使うと等値線にラベルが付きません。ラベルを付けたいときはgxout contourで上書きすることになります。この逆はできません。
!:カラーバーはshadedは図を描いてからga>run cbar.gs (or cbarn.gs or cbarc.gs)とするとでます。 cbarn.gsは指定しないとやたらとでかいカラーバー
が表示されます。cbarn.gsは場所、大きさを指定することができます。
書式: run cbarn.gs <scale > <num > <xmid > <ymid >
scale.......カラーバーの大きさを決めます。defaultは1.0です。半分の大きさにしたければ0.5にします。
num.........カラーバーを縦にするのか、横にするのかを決めます。縦にする場合は1、横にする場合は0とする。
xmid,ymid...カラーバーの位置の指定をする。数字はカラーバーの中心がvirtual pageのどこに位置するのかを指定している。
入力例:ga>run cbarn.gs 0.8 0 4.0 0.8
!: カラーバーを描くものとして、cbar.gs, cbarn.gs, cbarc.gsと言うのがあります。一番使い勝手のいいのはcbarn.gsだとおもいます

注意　デフォルトで.gsスクリプトがgradsに入っている模様　cbarn.gsが使える

from http://mausam.hyarc.nagoya-u.ac.jp/~hatsuki/grads/node8.html

登録: 投稿 (Atom)