格子化データの作成のため、格子点番号(i,j)を割り振る。この例では一度メッシュとする。
[2015年 5月 25日 月曜日 10:15:07 JST]
[~/Sato.Manda.Tachibana.FOG/ICOADS.Map]
[am@aofd165]
$ set.grid.number.sh
Shell script, set.grid.number.sh starts.
Input: /work4/data/ICOADS/2.5/ECS_SCS_JS/IMMA.MANDA107587.1748.08-1910.08.txt
Output: output.set.grid.number/IMMA.MANDA107587.1748.08-1910.08.grid.txt
-rw-rw-r-- 1 am am 2.5M 5月 25 10:15 output.set.grid.number/IMMA.MANDA107587.1748.08-1910.08.grid.txt
.....
.....
.....
#で始まるヘッダー行を除いて出力させてみる。
$ awk '{if ($1!="#") print $0}' output.set.grid.number/IMMA.MANDA107587.2014.08-2015.02.grid.txt |head
YR MO DY HR LAT LON VI VV AT WBT DPT SST B1 idx jdx
2014 8 1 0.00 38.80 132.00 -9 98 28.0 -999.9 -999.9 -999.9 82 18 17
2014 8 1 0.00 35.10 121.10 -9 95 28.0 26.0 25.3 28.0 51 7 14
2014 8 1 0.00 30.30 123.00 -9 95 26.3 -999.9 25.8 27.0 3 9 9
2014 8 1 0.00 27.30 123.60 -9 96 28.0 26.0 25.3 -999.9 73 9 6
2014 8 1 0.00 24.40 118.00 -9 97 28.0 26.0 25.2 29.9 48 4 3
2014 8 1 0.00 23.60 122.10 -9 99 28.0 26.0 25.0 29.0 32 8 2
2014 8 1 0.00 23.30 117.90 -9 97 25.0 22.0 20.5 27.1 37 3 2
2014 8 1 0.00 22.70 116.70 -9 98 29.2 -999.9 24.1 29.0 26 2 1
2014 8 1 2.00 29.10 123.20 -9 96 28.6 26.6 25.8 28.0 93 9 8
2014 8 1 3.00 29.50 131.60 -9 97 30.0 27.0 26.0 28.0 91 17 8
..........
1 YR 5 year UTC 1600 2024 year
2 MO 3 month UTC 1 12 month
3 DY 3 day UTC 1 31 day
4 HR 6 hour UTC 0 23.99 hour
5 LAT 7 latitude -90.0 90.0 deg N
6 LON 8 longitude 0. 359.0 deg E
7 VI 2 VV indic. 0 2 *
8 VV 3 visibility 90 99 *
9 AT 6 air temp. -99.9 99.9 deg C
10 WBT 6 web-bulb temp. -99.9 99.9 deg C
11 DPT 6 dew-point temp. -99.9 99.9 deg C
12 SST 6 sea surfc temp. -99.9 99.9 deg C
13 B1 3 1 deg. box num. 0 99
14 idx Grid number in longitudinal direction
15 jdx Grid number in meridional direction
set.grid.number.shにFortranプログラム, set.grid.number.f90が埋め込まれており、自動でコンパイル、実行を行う。
wlon, elon, slat, nlatの範囲内のデータのみ抜き出し、南西角の座標 (wlon, slat)を原点として、格子間隔dlon, dlatで格子を切った場合に、各データがどの格子に存在するか判別して、その結果をidx,jdxに代入し出力ファイルに書き出す。
$ cat set.grid.number.sh
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.fish.nagasaki-u.ac.jp
# Directory: /work1/am/Sato.Manda.Tachibana.FOG/ICOADS.Map
#
# Revision history:
# This file is created by /usr/local/mybin/nbscr.sh at 17:44 on 03-15-2015.
# Reference
# R2.5-imma_short.pdf
echo "Shell script, $(basename $0) starts."
echo
exe=$(basename $0 .sh)
indir=/work4/data/ICOADS/2.5/ECS_SCS_JS
prefix=IMMA.MANDA107587
inlist=$(ls $indir/${prefix}*.txt )
outdir="output.$(basename $0 .sh)"; mkdir -vp $outdir
wlon=115
elon=132
slat=22
nlat=45
dlon=1
dlat=1
#
# Create program source file
#
src=${exe}.f90
cat <<EOF>$src
program $(echo $exe |sed -e "s/\./_/g")
character(len=500)::infle,ofle
integer lentrm
integer vyr,vmo,vdy,vvi,vvv,vb1
real vhr,vlat,vlon,vat,vwbt,vdpt,vsst
character(len=5000)::strm
real wlon, elon, slat, nlat, dlon, dlat
integer,parameter::imiss=-999
namelist /para/infle, ofle, wlon, elon, slat, nlat, dlon, dlat
read(*,nml=para)
open(10, file=infle, action='read')
print '(A,A)','Input: ',trim(infle)
print '(A,A)','Output: ',trim(ofle)
open(20,file=ofle)
write(20,'(A,A)')'# Input = ',trim(infle)
write(20,'(A,A)')'# '
write(20,'(A)')'# Idx Abbr FW ElementDesc Min. Max. Units* MissingVal IMMA Comp./Sec.'
write(20,'(A)')'# --- ---- -- ---------------- ----- ------ ------ ---------- ---------------'
write(20,'(A)')'# 1 YR 5 year UTC 1600 2024 year -9999 core/location '
write(20,'(A)')'# 2 MO 3 month UTC 1 12 month -99 core/location '
write(20,'(A)')'# 3 DY 3 day UTC 1 31 day -99 core/location '
write(20,'(A)')'# 4 HR 6 hour UTC 0 23.99 hour -99.99 core/location '
write(20,'(A)')'# 5 LAT 7 latitude -90.0 90.0 deg N -999.99 core/location '
write(20,'(A)')'# 6 LON 8 longitude 0. 359.0 deg E -9999.99 core/location '
write(20,'(A)')'# 7 VI 2 VV indic. 0 2 * -9 core/regular '
write(20,'(A)')'# 8 VV 3 visibility 90 99 * -99 core/regular '
write(20,'(A)')'# 9 AT 6 air temp. -99.9 99.9 deg C -999.9 core/regular '
write(20,'(A)')'# 10 WBT 6 web-bulb temp. -99.9 99.9 deg C -999.9 core/regular '
write(20,'(A)')'# 11 DPT 6 dew-point temp. -99.9 99.9 deg C -999.9 core/regular '
write(20,'(A)')'# 12 SST 6 sea surfc temp. -99.9 99.9 deg C -999.9 core/regular '
write(20,'(A)')'# 13 B1 3 1 deg. box num. 0 99 -99 ICOADS/Box '
write(20,'(A)')'# 14 idx'
write(20,'(A)')'# 15 jdx'
write(20,'(A,f10.5)')'# wlon = ',wlon
write(20,'(A,f10.5)')'# elon = ',elon
write(20,'(A,f10.5)')'# dlon = ',dlon
write(20,'(A,f10.5)')'# slat = ',slat
write(20,'(A,f10.5)')'# nlat = ',nlat
write(20,'(A,f10.5)')'# dlat = ',dlat
write(20,'(A,f10.5)')'# '
write(20,'(A,f10.5)')'# VV (horizontal visibility at the surface in kilometers) according to WMO Code 4377 from'
write(20,'(A,f10.5)')'# which, in reporting visibility at sea, WMO (2009a; Reg. 12.2.1.3.2) states that the decile'
write(20,'(A,f10.5)')'# 90-99 shall be used (moreover Reg. 12.2.1.3.1: when the horizontal visibility is not the'
write(20,'(A,f10.5)')'# same in different directions, the shortest distance shall be given for VV):'
write(20,'(A,f10.5)')'# 90 - less than 0.05 kilometer'
write(20,'(A,f10.5)')'# 91 - 0.05'
write(20,'(A,f10.5)')'# 92 - 0.2'
write(20,'(A,f10.5)')'# 93 - 0.5'
write(20,'(A,f10.5)')'# 94 - 1'
write(20,'(A,f10.5)')'# 95 - 2'
write(20,'(A,f10.5)')'# 96 - 4'
write(20,'(A,f10.5)')'# 97 - 10'
write(20,'(A,f10.5)')'# 98 - 20'
write(20,'(A,f10.5)')'# 99 - 50 or more'
write(20,'(A,f10.5)')'# VI shows whether VV was:'
write(20,'(A,f10.5)')'# 0 : estimated (or unknown method of observation)'
write(20,'(A,f10.5)')'# 1 : measured'
write(20,'(A,f10.5)')'# 2 : fog present (obsolete)'
write(20,'(A,f10.5)')'# Background: The “Cloud height and visibility measuring indicator” from IMMT-4 is'
write(20,'(A,f10.5)')'# separated into independent indicators in IMMA format, HI (see field 39) and VI.'
write(20,'(A,f10.5)')'# [Note: When VI=2, and VV=93, it meant that fog was present and visibility was'
write(20,'(A,f10.5)')'# not reported (NCDC 1968). This “fog present” combination of VI=2 with VV=93'
write(20,'(A,f10.5)')'# appears to originate from “overpunch” procedures that took effect in the IMMPC'
write(20,'(A,f10.5)')'# format around 1966 (see Table B1) as translated into the TDF-11 format.]'
do
read(10, '(A)', iostat=ios) strm
if(ios<0)exit
read(strm,*)vyr,vmo,vdy,vhr,vlat,vlon,vvi,vvv,vat,&
& vwbt,vdpt,vsst,vb1
idx=(vlon - wlon)/dlon + 1
jdx=(vlat - slat)/dlat + 1
if(vlon > elon .or. vlon < wlon )idx=imiss
if(vlat > nlat .or. vlat < slat )jdx=imiss
if(idx==imiss .or. jdx==imiss)then
idx=imiss
jdx=imiss
endif
if(idx/=imiss .and. vvv>0.0 )then !.and. vvi==1)then
write(20, '(i5,x,i3,x,i3,x,f6.2,x,f7.2,x,f8.2,x,i2,x,i3,x,f6.1,&
x,f6.1,x,f6.1,x,f6.1,x,i3,x2i5)') &
& vyr,vmo,vdy,vhr,vlat,vlon,vvi,vvv,vat,vwbt,vdpt,vsst,vb1,&
& idx,jdx
endif
enddo
end program
function lentrm(str)
character str*(*)
do lentrm=len(str),1,-1
if(str(lentrm:lentrm) .ne. ' ') return
enddo
end function
EOF
#
# Compile the program
#
fort=ifort
opt="-traceback -check"
${fort} ${opt} ${src} -o ${exe}
if [ $? -ne 0 ]; then
echo
echo Compile Error.
exit 1
fi
namelist=${exe}.namelist.txt
#
# Run the program
#
for in in $inlist; do
out=$(basename $in .txt).grid.txt
cat <<EOF>$namelist
¶
infle="$in"
ofle="${outdir}/$out"
wlon=$wlon
elon=$elon
slat=$slat
nlat=$nlat
dlon=$dlon
dlat=$dlat
&end
EOF
${exe} < $namelist
ls -lh ${outdir}/$out
done
echo "Done $0"
echo