はじめての気象データ処理

このブログでは、気象を学ぶ方でも超初心者を対象として気象データ解析のための基礎基本を解説します。主にFORTRAN、Python、GMT、GrADSのTipsをご紹介します。

Fortranでアメダスデータを処理して平均気温を求める

みなさんこんにちは。おこめまんです。

前回はFortranで大間アメダス10分間値データを読み込むことに挑戦しました。

今回は、Fortranで読み込んだ1か月分の気温データから平均気温を求めてみたいと思います。

今回使用するデータ

今回使用するデータは以前スクレイピングによってゲットした「アメダス大間10分間値データ」です。

アメダス10分間値データをゲットする方法は以下をご覧ください。

気象庁アメダス10分間値をダウンロードする方法 

 

<ファイル名>

「amd10_31001_201801.csv

このcsvファイルの中身はこんな感じです。

f:id:met_learner:20191217003801p:plain

amd10_31001_201801.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 

と画面に出れば正解です。

 

最後に

いかがでしたか。平均気温は、データ処理のはじめによく扱われます。時間があれば、平均気温のほかに、平均風速も求めてみてください。

それでは、また次回お会いしましょう。おこめまんでした。