戻る
気象庁CD(月報)からのアメダス読みだし
(Fortran90)
read_amedas.f90 プログラムのダウンロード
Subject: アメダスデータ読みだしソフト
Date: Fri, 22 May 2009 19:29:13 +0900
研究室の皆さん、特に薛君
玉川です。薛君の学会発表にむけて、今アメダスデータを読んでいると思いますが、
読みだしのソフトを書いてみました。
急いで書いているので、何かバグがあるかも知れませんが、
一応、動いているようなので、送ります。
添付のソフト read_amedas.f90 を
どこかにおいて、そこに、データCDの amedas フォルダを
コピーしてください。
そして、read_amedas.f90 をコンパイルして
./a.out > datafile
と動かすと、若干のメッセージの後に
year, month, day, hour
と聞かれますから、必要な 年月日をいれてください。
(もちろん 2007 8 まではこれしかないですが)
そうすると datafile に
write(*,*) find, pnum(i), lat(i), lon(i), windsp, winddir, temp, prec, sun
でメッセージが書かれます。
find : 1 しかないはず
pnum(i) : 地点番号
lat(i) : latitude
lon(i): longitude
windsp : 風速 (m/s)
winddir : 風向 (0-360度) 0 が北，90度が東
temp : 気温
prec : 雨量
sun : 日照
です。
若干注意して使ってください。
この出力を加工すれば、GMTで図が書けると思います。
例えば
awk '{print $4, $3, $5, $6}' < datafile > datafile2
で、lon lat windsp winddir のデータが並ぶ出力が datafile2 に
得られると思います。
P.S. バグってたら直してください。
--------------------------------------------------------------------------------
〒 501-1193 岐阜市柳戸1-1
岐阜大学 流域圏科学研究センター
玉川一郎 (e-mail: tama@green.gifu-u.ac.jp)
phone/fax 兼用: 058-293-2430
WWW: http://tama.green.gifu-u.ac.jp/~tama
--------------------------------------------------------------------------------
! maybe ....
! bt I. Tamagawa 2009 May
implicit none
integer, parameter :: maxpoint = 20000
integer :: pnum(maxpoint), windtemp(maxpoint), pointnum
real :: lat(maxpoint), lon(maxpoint)
character(len=1024) :: amdmaster, datafile
integer :: i, find
integer :: year, month, day, hour
real :: windsp, winddir, temp, prec, sun
amdmaster = 'amedas/amdmaster.index'
write(0,*) trim(amdmaster)
! read madmaster.index
call readindex( amdmaster, pnum, lat, lon, windtemp, maxpoint, pointnum)
write(0,*) 'read ', trim(amdmaster), pointnum
call select_newest( pnum, lat, lon, windtemp, maxpoint, pointnum)
write(0,*) 'select last record for each point', pointnum
write(0,*) 'year, month, day, hour'
read(*,*) year, month, day, hour
do i=1, pointnum
write(datafile, "('amedas/hourly/',i0.4, i0.2,'/area',i0.2,'/ah',i0.5,'_', i0.4, i0.2,'.csv')") &
year, month, &
int(pnum(i)/1000), pnum(i), year, month
! write(*,*) datafile
if( windtemp(i) == 1 ) then
call readdata( datafile, year, month, day, hour, windsp, winddir, temp, prec, sun, find )
if(find == 1 ) then
write(*,*) find, pnum(i), lat(i), lon(i), windsp, winddir, temp, prec, sun
end if
end if
end do
end
subroutine readindex( amdmaster, pnum, lat, lon, windtemp, maxpoint, pointnum)
character(len=*) :: amdmaster
integer :: maxpoint, pointnum
integer :: pnum(maxpoint), windtemp(maxpoint)
real :: lat(maxpoint), lon(maxpoint)
integer :: lat1, lon1
real :: lat2, lon2
integer :: ios=0
integer :: count = -2 ! skip first 2 line
character(len=1024) :: buf
integer :: prec, windspeed, temp, sunshine, snow_depth
open( file=trim(amdmaster), &
status='old', unit = 10, err = 900 )
write(0,*) 'open ', trim(amdmaster)
do while(ios == 0 )
read(10,'(a)', iostat=ios ) buf
count = count +1
! write(0,*) count
if( count > maxpoint) then
write(0,*) ' Too many data, change maxpoint'
stop
end if
! write(0,*) count, 'QQ: ', buf
if( count > 0 ) then
! write(0,*) count, buf(1:5), ':', buf(75:90)
read(buf(1:5), '(i5)') pnum(count)
read(buf(75:90), * ) lat1, lat2, lon1, lon2
lat(count) = lat1 + lat2/60.0
lon(count) = lon1 + lon2/60.0
read(buf(103:112), *) prec, windspeed, temp, sunshine, snow_depth
if( windspeed == 1 .AND. temp == 1 ) then
windtemp(count) = 1
else
windtemp(count) = 0
endif
end if
!write(0,*) count, buf(103:112)
! write(0,*) count, pnum(count), lat(count), lon(count), &
! prec, windspeed, temp, sunshine, snow_depth
end do
pointnum = count
close(10)
! write(0,*) trim(buf)
return
900 write(0,*) 'Can not open ', amdmaster
stop
end
subroutine select_newest( pnum, lat, lon, windtemp, maxpoint, pointnum)
integer :: pointnum
integer :: pnum(maxpoint), windtemp(maxpoint)
real :: lat(maxpoint), lon(maxpoint)
integer :: i, j
j=1
do i=2, pointnum
if( pnum(i) /= pnum(i-1) ) then
j=j+1
end if
! write(0,*) 'i,j', i, j
pnum(j) = pnum(i)
lat(j) = lat(i)
lon(j) = lon(i)
windtemp(j) = windtemp(i)
end do
write(0,*) pointnum, '==>', j
pointnum = j
end
subroutine readdata( datafile, year, month, day, hour, windsp, winddir, temp, prec, sun, find )
implicit none
character(len=*) :: datafile
integer :: year, month, day, hour
real :: windsp, winddir, temp, prec, sun
character(len=1024) :: buf
integer :: i, ios=0
integer :: dday, dhour, dpfl, dwdfl, dwsfl
real :: dprec, dwinddir, dwindsp, dtemp, dsun
integer :: dtfl, dsfl
integer :: find
find = 0
open( file = datafile, unit = 20, status='old', err = 999 )
do i=1,4 ! skip 4 lines
read(20,'(a)') buf
end do
do while( ios== 0 )
read(20, *, iostat = ios, end=910) dday, dhour, dprec, dpfl, dwinddir, dwdfl, dwindsp, dwsfl, &
dtemp, dtfl, dsun, dsfl
write(0,*) dday, dhour, dprec, dwinddir, dwindsp, dtemp, dsun
if( dday == day .AND. dhour == hour ) then
! write(*,*) 'FIND!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
windsp = dwindsp
! 16: N,
if(dwinddir /= 16 ) then
winddir = (dwinddir)*(360.0/16.0)
else
winddir = 0
end if
temp = dtemp
prec = dprec
sun = dsun
find = 1
close(20)
return
goto 910
end if
end do
910 close(20)
if(find == 0 ) then
write(0,*) 'Can not find data'
windsp = -9999 ; winddir = -9999 ; temp = -9999 ; prec = -9999 ; sun = -9999
stop
endif
return
999 write(0,*) 'Can not open ', datafile
windsp = -9999 ; winddir = -9999 ; temp = -9999 ; prec = -9999 ; sun = -9999
end