高山試験地AWSの月報処理のために

岐阜大学流域圏科学研究センターの高山試験地には、いくつもの観測サイトがあります。
本研究室では、高山常緑針葉樹林サイト(TKC)と試験地庁舎AWSの 2つを管理しています。
ここでは、そのうちAWSの方の月毎の統計データ処理を試みます。

AWS のページ
日報処理のソフトやや、やる気なく＆忙しい時に書いたので汚い。

まずは、何をするのか。。

月報とは何か、というと、月毎に、各日のデータ(平均値、最大値、最小値、あるいは降水量の積算値)を出力するものである。
そうすると、まず、いつからいつまでのデータをどの間隔で処理するのか、などという基本的な情報を処理せねばならない。
日付データの処理は、日時の処理の例に基本的なものを用意してあるが、月の長さというようなものを出力するのは作ってなかったので、
これから作成する。
date_functions.f95 に月の長さサブルーチンを追加したもの
本質的にはこれだけ↓


subroutine days_of_month(year, month, days)
! number of days in year/month
integer :: year, month, days
  logical, external :: leepyear

 integer :: m(12)
  m(1) = 31; m(2) = 28; m(3) = 31; m(4) = 30
  m(5) = 31; m(6) = 30; m(7) = 31; m(8) = 31
  m(9) = 30; m(10) = 31; m(11) = 30; m(12) = 31

 if(leepyear(year)) then
  m(2) = 29
 end if

 days = m(month)


end


logical function leepyear( year )
 integer :: year

 if( mod(year, 4) > 0 ) then
   leepyear  = .FALSE.
 else if( mod(year, 100) > 0 ) then
   leepyear = .TRUE.
 else if( mod(year, 400) > 0 ) then
   leepyear = .FALSE.
 else
   leepyear = .TRUE.
 end if

return
end

元データが、前10分値なので、分を基本単位に処理を考える。
年と月が、与えられると、開始年月日時刻と、処理期間(分)を返すようにするには、 test1.f90 のようにすれば良い。

次に、読み込むデータファイルは２つあって、

http://sateco-archive.green.gifu-u.ac.jp/AWS1/Raw/00_All_AWS1_10Min.csv
http://sateco-archive.green.gifu-u.ac.jp/AWS1/Raw/00_All_AWS2_10Min.csv

"TOA5","AWS1","CR1000","7964","CR1000.Std.14","CPU:AWS1_112.CR1","34813","TBL110"
"TIMESTAMP","RECORD","Wind_spd_WVc(1)","Wind_spd_WVc(2)","Air_Temp_Avg","Humidity_Avg","Sun_Wm2_Avg","Rain_Tot","Press_Avg","Issui_Tot"
"TS","RN","","","","","","","",""
"","","WVc","WVc","Avg","Avg","Avg","Tot","Avg","Tot"
"2007-07-20 12:20:00",1,0.155,70.95,17.37,93.2,107.7723,0,862.1425
"2007-07-20 12:30:00",2,0.079,69.33,17.65,91.1,51.95923,0,861.0716
"2007-07-20 12:40:00",3,0.011,265.2,18.02,97.9,80.16187,5,861.0716

のようなフォーマットであることを考える。(のちに項目が増える)

データは、"," で区切られており、ASCII文字なので、
1行を文字列で読んで、"," を見て分解( "" の処理くらいは付ける) して、文字列からのREADで数値にするだけであるが、
開始時期や欠測の位置の違う２つのファイルを併せて処理するので、ちょっと工夫が必要である。

10分毎のデータであることを積極的に使うと、
必要期間の10分毎のデータを格納する配列を作っておいて、そこにデータをコピーし、後でいろいろ統計処理をする手がある。

あるいは、前から読んで処理するのに拘って

読んだデータを、１日毎に統計計算してその結果のみを残す。

という処理も可能である。
後者の方なら、測定間隔が30分でも1時間でもソフトの変更なしに使えそうで惹かれるが、
複雑な統計処理を、データの読みだしと同時に行うことになり、やや込み入ったソフトになる。
ここで、降水量を考えると、月毎には、２時間や２４時間降水量の最大値なども欲しくなるし、
もちろん、気温など他のデータも日単位ばかりではなく、月単位での最大値や最小値なども欲しいので、結構込み入って来る。
どうしようか。。。。
日報はこのパターンで作ったんで、実はこいつをちょこっといじってすますこともできなくはなさそうなんだけれど。。。。
ちなみに、風速の平均は、ベクトルにしてから取るべきだけど、最大風速は風速そのものの最大値だし、降水量は積算だし。。。。

もし、10分毎にデータを残すとすると、
6個/h * 24 * 31 = 4464個/月のデータを各測定項目(および計算項目)毎に用意することになる。
これくらいなら、１年の統計でもたいしたことはないから、今回はこっちで行こうと思う。

ということで、まずは２つのデータファイルから該当部分を読み出して、
１つの配列にいれるプログラムがこれ monthly_report_tmp.f90 。
ちょっと解説を入れると、メインプログラム以外に、以下のサブルーチン、関数を作ってある。

subroutine testout( data, total_num, ch1, shift1, ch2, shift2 ) テスト出力
subroutine initdata( data, total_num ) データ配列初期化 (NODATAを入れる)
subroutine read_copy( file1, ch1, shift1 。。。データをファイルから読んで、必要部分を配列 data(*) にコピー
subroutine devidecsv( buf, eachcol, maxch, num ) カンマ区切りのアスキー文字列を、部分文字列に分解、" " の処理は付けたはず。。
subroutine setparam( year, month, start, duration ) 年月から、処理開始日時および計算期間(分)をもとめる。
subroutine days_of_month(year, month, days) １ヵ月の日数を計算
subroutine datedif( year1, month1, day1, year2, month2, day2, diffday ) 日時の差 (日) を計算
integer function calc_diffmin( now, start ) ２つの日時の差 (分) を計算
logical function leepyear( year ) 閏年かどうかを判定
integer function dayofyear(year, month, day) 年の通算日を計算

また、データや日時の表現には、構造体を使ってまとめてある。
他に、引数の処理には、メインプログラムの最初の方にある
command_argument_cout, get_command_argument を使った( Fortran 2003 標準らしい。既に gfortran44 で使える)
時間の処理を分単位にして、整数で扱っているので、計算誤差がなく、
== の成立や、区切りを微小に越えたりすることを深く考えないで良いようになっている。

コンパイルして、データファイルを同じディレクトリに置いておいて、
./a.out 2010 1 とかして、実行すると、その月のデータを横に並べた表にしてくれる。
プログラム中の extent_min を１日とかにすれば、前後に１日分の余裕をとって、データを配列にコピーしてくれる機能もつけた。(２４時間降水量等対応用)
あとは、統計部分を作って出力を整えるだけだ(2010.3.24)

2010.3.25, 今日は、あまりやってない。
統計処理の前に、風向風速から、UV成分、気温、相対湿度から、比湿、水蒸気圧への変換
そういうのを書くための、ファイル情報 (ファイル名、チャンネル番号など )の module への分離、24時間積算降水量を移動平均的に出す、などという機能を
作ったところで、止まってしまった。
今日の成果品
積算処理や、風速等の変換処理でもいちいち NODATA の確認をしているあたりがちょっとした注意点。(当然の処理だけど)

あとは、１日毎の統計、月での統計 (積算、平均、最大最小) をとって出力すれば、普通の月報にはなる。
でも、ここまできたら、 windrose とか、スカラー量の頻度分布とか取っておきたい。
それに、日射同士、降水量同士で統計もみたいなぁ。
誰のためにするんだか、って感じだけど。

ということで、まずは、ある時間間隔で、最大、最小、平均、合計くらいの統計処理を行う機能をつけた
monthly_report_tmp3.f90 また、コンパイルが通る程度に、修正してある。

それから、頻度分布を作るサンプルサブルーチンを作っておいたこっち (2010.4.1)
もうちょっと、作ると monthly_report_tmp4.f90 となって、後は頻度分布と風配図、それと表示の整理 (2010.4.1)

結局、出力用に年/月 (たとえば 2010/03) のディレクトリを作成し、その下に
月の統計、日の統計、ヒストグラム、風配図などを csv ファイルで出力する仕様にした。
monthly_report_tmp6.f90 (2010.4.6)

で、さらに、gnuplot で作図するスクリプトを吐いて、call system で実行するとか、
HTML ファイルも出力するようにした。gnuplot と convert (ImageMagick) を call system する。
monthly_report_tmp7.f90
あと、windorose のレーダチャートを追加して、HTML の体裁を整え、 cron 実行に備えて、gnuplot 等をフルパスに直せば、多分出来上がる。
なかなか面倒だ。。。。(2010.4.6)

もういい加減おわりにしたい。
と言う訳で、ちょっと頑張って作った。できたのは、 monthly_report_tmp8.f90 。
日データの時系列、および、頻度分布や風配図をプロットして、HTMLにまとめるとともに、 csv ファイルへのリンクを貼ってある。
計算結果のサンプル

その後、2011年1月に積雪計の処理を追加したバージョンを作った。
monthly_report_tmp10.f90 。しかし、ファイル名の tmp はついたまま。。。
表を吐いたり、図を描くスクリプトをだしたりするところがとても汚い
作者としては、整理したいが、その気力はない。。。

なのに、さらに2011年11月に、地温、土壌水分のセンサーが加わった。
とりあえず作ったのが、 monthly_report_tmp11.f90 。
時系列図の部分に追加が必要だが。。。。(2011.12.09)
時系列図を追加した (2011.12.14) monthly_report_tmp12.f90 。

ちなみに、tailコマンドのオプションが、Linux の　GNU版と FreeBSD の　BSD版で違うようで、tail +3 みたいな書き方は、Linux では tail -n +3 だったかな？こんな感じで違う。注意。

AWS の月報ページ

サンプルプログラムのページへ, 玉川研究室のページへ, 玉川のページへ