Fortranでアメダスデータを処理して平均気温を求める
みなさんこんにちは。おこめまんです。
前回はFortranで大間アメダス10分間値データを読み込むことに挑戦しました。
今回は、Fortranで読み込んだ1か月分の気温データから平均気温を求めてみたいと思います。
今回使用するデータ
今回使用するデータは以前スクレイピングによってゲットした「アメダス大間10分間値データ」です。
アメダス10分間値データをゲットする方法は以下をご覧ください。
<ファイル名>
「amd10_31001_201801.csv」
このcsvファイルの中身はこんな感じです。
各列はそれぞれ、
「年」「月」「日」「時」「分」「秒」「降水量」「気温」「平均風速」「風向(平均風速)」「最大瞬間風速」「風向(最大瞬間風速)」「日照時間」
となっています。
よって、最初の1行目には
2018年1月1日0時10分0秒の降水量、気温、平均風速、風向(平均風速)、最大瞬間風速、風向(最大瞬間風速)、日照時間が書かれています。こんな感じ1月31日まで1か月分のデータがあります。
使用した言語
Fortran95
Fortranプログラム
以下を書き込んだプログラム「open_amedas_data.f95」というプログラムを作成します。
program open_amedas_data
implicit none
integer :: year, month, day, hour, minute, second
integer ::counter_temp
real :: prc, temp, ws, wd, gws, gwd, sun
real :: sum_temp, ave_temp
character :: ifile*100
ifile="amd10_31001_201801.csv"
open(30, file=trim(ifile), status='old')
sum_temp = 0
counter_temp = 0
do
read (30, *, end = 99, err=98) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
if (temp > -900.0) then
sum_temp = sum_temp + temp
counter_temp = counter_temp + 1
end if
!write (6,*) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
end do
98 write (6, *) "ERROR"
stop
99 write (6, *) "Finish"
close (30)
write(6,*) "counter=", counter_temp
ave_temp = sum_temp/counter_temp
write(6,*) "ave_temp= ", ave_temp
stop
end program open_amedas_data
プログラムの説明
プログラムの説明は以前ご紹介した「Fortranでcsvファイルを開いて1行ずつ読み込む」から追加した内容を扱います。付け足した部分はプログラムをオレンジ色で着色しました。
--------------------------------------------------------
implicit none
integer :: year, month, day, hour, minute, second
integer ::counter_temp
real :: prc, temp, ws, wd, gws, gwd, sun
real :: sum_temp, ave_temp
character :: ifile*100
⇒ 計算処理する際に必要な変数を指定します。counter_tempはデータの個数(10分間値データ1か月分なので、欠測がなければ4464個)を表します。sum_tempは1か月分の気温データを合計したものに相当します。ave_tempは今回求めたい「平均気温」に相当します。
--------------------------------------------------------
sum_temp = 0
counter_temp = 0
⇒ 単純に気温に足せばいいじゃんと思う人もいることでしょう。しかし、はじめのsum_tempとcounterの値を指定していなければうまく計算することができません。
あなたも、「+4は?」とだけ言われても困ってしまうでしょう。「何」に4を足すのかをはじめに指定しなければなりません。Fortanでは、ループを回す前に最初の値は「0」だと教えてあげましょう。
--------------------------------------------------------
if (temp > -900.0) then
sum_temp = sum_temp + temp
counter_temp = counter_temp + 1
end if
⇒ readで読み込んだ要素のうち、気温に注目して積算していきます。if (temp > -900.0) thenとしているのは、アメダスデータの欠測値を「-999.9」としているので、このときのデータをはじくためです。「もし気温が-900.0℃よりも大きければ、sum_temp = sum_temp + tempとcounter_temp = counter_temp + 1を実行してください」という意味です。正しく観測されていれば、気温は-900.0℃よりは大きい値のはずですよね。counter_tempでは、欠測値を除いたデータ数をカウントしています。欠測がなければ4464になるはずです。
--------------------------------------------------------
write(6,*) "counter=", counter_temp
⇒ さて、doループが終わりました。このとき、sum_tempには1か月分の気温データを積算した値が、counter_tempには欠測値を除いたデータ数が相当していることになります。ここでは、欠測を除いたデータ数を画面に表示させています。
--------------------------------------------------------
ave_temp = sum_temp/counter_temp
write(6,*) "ave_temp= ", ave_temp
⇒ ここで、ついに1か月の平均気温を求めます。平均気温は、積算した気温をデータ数で割れば求まりますよね。計算した平均気温はave_tempに入ります。最後に、write文を書いて、画面に平均気温を表示させましょう。
プログラムを実行して、
counter= 4464
ave_temp= 0.185748115
と画面に出れば正解です。
最後に
いかがでしたか。平均気温は、データ処理のはじめによく扱われます。時間があれば、平均気温のほかに、平均風速も求めてみてください。
それでは、また次回お会いしましょう。おこめまんでした。