[気象庁の観測データを使う (序論)]のページで、 気象庁ウェブサイトの「過去の気象データ ダウンロード」の機能をつかって、 例として、館林の2年ぶんの日ごとの値をダウンロードした。 (ダウンロードしたものを、つぎのところにも置いておく。)
ここでは、そのデータを R で読んで、図化や基本的な統計処理をこころみる。 【説明は Windows 上の R 4.0.3 での実行結果によるが、他の OS でも同様にできるはずである。 ただし Linux 版の R では R処理系とその外のエディタとを行き来する処理になるだろう。】
作業用ディレクトリをつくっておく。 気象庁からダウンロードしたCSVファイルを、 ファイル名を「daily_tatebayashi.csv」として、作業用ディレクトリに置く。 Rを起動し、(Windows上のばあい)「ファイル」の中の「ディレクトリの変更」を実行して、 作業用ディレクトリをえらんでおく。 (「ディレクトリの変更」は、Rを起動するたびに実行する必要があるようだ。)
Rには「データフレーム」というデータ構造がある。 おおざっぱには、縦横の表と思えばよい。 ただし、行はおたがいに同質のデータだが、列はそれぞれにちがった意味をもつことがある。
データフレームには名まえが必要で、名まえは自由につけられるが、 ここでは、プログラム内で1つだけのデータフレームをつかうので、 データフレーム名を(データフレームであることをわすれないように)「df」としてみる。
Rには、CSVファイルをデータフレームに読みこむ「read.csv」という関数がある。 (関数名がピリオド(「.」)をふくんでいることに注意。) これをつかって、Rのコンソールからつぎのようなコマンドを打ちこめば、 CSVファイル「daily_tatebayashi.csv」の内容をデータフレーム「df」に読みこめる。
df <- read.csv("daily_tatebayashi.csv", skip=6, header=F, encoding="Shift_JIS")
この作業をするだけの R プログラムをつぎのところに用意した。 これをダウンロードして作業ディレクトリに置き、 (Windows上のRのばあい) 「ファイル」の「スクリプトを開く」で これを読みこみ、 「編集」から「全て実行」で実行してもよい。
【Rの変数名は ASCII コードであらわせるアルファベット・数字・ 下線(「_」、アンダースコア)で構成する。はじめの字はアルファベットにする。 大文字と小文字は R では区別されるが、人がまちがえやすいので、 わたしの教材では、自動的に大文字になるもの以外は、小文字だけをつかうことにする。】
CSVファイルを表と見たときの各列の内容は、 「df$V1」...「df$V18」という形で参照できる (データフレーム名、「$」(ドル記号)、大文字の「V」、列番号の数字が、空白をはさまずにつなげられている)。 Rのコンソール画面から、単に「df$V4」などと打ちこむと、その列の数値のリストが表示される。
[気象庁の観測データを使う (1)]のページでつくった館林の日ごとの値のばあい、 列番号、Rでの列名、表計算での列記号、気象データとしての項目は、つぎのように対応している。 (つけなおした列名は、つぎであつかう。)
| 列番号 | Rでの列名 | 表計算での列記号 | 項目 | つけなおしたRでの列名 |
|---|---|---|---|---|
| 1 | df$V1 | A | 年 | df$year |
| 2 | df$V2 | B | 月 | df$month |
| 3 | df$V3 | C | 日 | df$day |
| 4 | df$V4 | D | 平均気温 | df$temp |
| 5 | df$V5 | E | 平均気温の品質情報 | df$q_temp |
| 6 | df$V6 | F | 平均気温の均質番号 | df$h_temp |
| 7 | df$V7 | G | 最高気温 | df$tempmax |
| 8 | df$V8 | H | 最高気温の品質情報 | df$q_tempmax |
| 9 | df$V9 | I | 最高気温の均質番号 | df$h_tempmax |
| 10 | df$V10 | J | 最低気温 | df$tempmin |
| 11 | df$V11 | K | 最低気温の品質情報 | df$q_tempmin |
| 12 | df$V12 | L | 最低気温の均質番号 | df$h_tempmin |
| 13 | df$V13 | M | 降水量 | df$precip |
| 14 | df$V14 | N | 降水量の品質情報 | df$q_precip |
| 15 | df$V15 | O | 降水量の均質番号 | df$h_precip |
| 16 | df$V16 | P | 日照時間 | df$sundur |
| 17 | df$V17 | Q | 日照時間の品質情報 | df$q_sundur |
| 18 | df$V18 | R | 日照時間の均質番号 | df$h_sundur |
列名が番号だとおぼえきれないので、実際の気象要素に対応する名まえにつけなおすことにする。 日本語ローマ字でもかまわないのだが、ここでは英語の略語にした。
列名の変更は、R の colnames という関数によってできる。 18列の名まえを一度に変更するプログラムをつぎのところに用意した。 ただし、1行に書いてしまったので、長い行になっている。 上記の readf_{sjis,utf8}.R を実行したあとに、これを読みこんで実行すれば、 以後は、つけかえられた列名をつかうことができる。
Rの作図機能をつかって、データをグラフにすることができる。 ここでは R の base にふくまれるグラフィックス機能をつかう。 R の plot は、2つの変数列を指定して散布図をつくるのが基本的な機能である。 しかし、1つの変数列だけを指定すれば、行を順にかぞえた番号を横軸にとったグラフをつくることができる。
たとえば、つぎのようにすれば、日照時間の時系列グラフをつくることができる。 (気温や降水量についても変数名がちがうだけで同様である。)
plot(df$sundur)
それにグリッドをかさねてみる。
grid()
図はつぎのようなものになる。
ここで、横軸は「Index」となっているが、これはデータの行を順にかぞえた番号で、 daily_tatebayashi.csv では、2018年12月1日を1とした とおしの日付にあたる。 これをわかりやすく表現しなおすこともやりたいのだが、 わたしがまだ R での日付のあつかいに熟練していないので、 ひとまずこのままにしておく。 (日付を日付としてあつかう方法をわたしが理解してプログラム例をつくれるようになったら、補足説明を書くことにしたい。)
ひとつひとつの観測値が、まる(〇)であらわされている。 記号(ここでは「まる」)の大きさを変えるには、たとえばつぎのようにする。
plot(df$sundur, cex=0.2)
記号を変えるには、たとえばつぎのようにする。
plot(df$sundur, pch=4)
日照時間の品質情報 df$q_sundur の値を見ると、ほとんど「8」なのだが、ときどき「4」がまざっている。 その日には、(おそらく機器点検があったので) 24時間の観測値がそろっていないのだ。 df$sundur にはその日の数値がはいっているが、それは統計処理からはずしたほうがよいだろう。 そこで、 df$q_sundur の値が 8 でない日の df$sundur の値を、欠損 (Rの表現で「NA」) に変えておくことにする。
この手つづきは、手順だて型のプログラム言語に共通なループの機能をつかうならば、 つぎのようなプログラムでできる。 「!=」は「ひとしくない」という、2つの数値から論理値を得る演算子、 「<-」は代入演算 (左の変数に右の計算結果を代入する)である。
なお、dfの内容の行数(各列のデータ個数)は、 たとえば、length(df$sundur) で得られる。 (他の列でも行数は同じだが、実際につかう列をあげておくことにする。 length(df) では 列数になってしまう。 ) 実際の数は、うるう年を1つふくむ2年間の日数なので、365×2+1 である。】
for(kt in 1:length(df$sundur)) {
if(df$q_sundur[kt] != 8) {
df$sundur[kt] <- NA
}
}
しかし、R のデータフレーム演算をつかうと、つぎのような短いプログラムで実行できる。
df$sundur <- ifelse(df$q_sundur != 8, NA, df$sundur)
【これは、左辺の各要素に、もし「df$q_sundurの値が8でない」という条件がなりたっていればNAを、 なりたっていなければ df$sundurの値を代入する、という意味である。 左辺はdfとは別の変数列(Rの用語でいう「ベクトル」)でもよい。】
気温の品質情報には、ときどき「5」がまざっている。 ここでは、気温や降水量についても、品質情報が「8」でない観測値は欠損にしておこう。
df$temp <- ifelse(df$q_temp != 8, NA, df$temp)
df$tempmax <- ifelse(df$q_tempmax != 8, NA, df$tempmax)
df$tempmin <- ifelse(df$q_tempmin != 8, NA, df$tempmin)
df$precip <- ifelse(df$q_precip != 8, NA, df$precip)
(もし上記のforループによるプログラム断片をつかうならば、「if」からの3行をコピーして 変数名を書きかえたものを追加すればよい。)
散布図 (英語では scatterplot) の機能をためしてみよう。 つぎのようにして、日照時間と最高気温の散布図をつくることができる。
plot(df$sundur, df$tempmax)
つぎのようにして軸のラベルを変え、表題をつけることができる。
plot(df$sundur, df$tempmax, xlab="日照時間 (時)", ylab="最高気温 (℃)", main="館林, 2018年12月～2020年11月")
R には、3つ以上のデータ列について、総あたりの散布図をまとめてつくる機能がある。 縦横とも列の数だけ正方形の枠がならび、対角線のところに列名が書かれ、 そのほかの枠に散布図がかかれた形になる。 このようなグラフ表現を、英語では scatterplot matrix という。 日本語表現としては「散布図方陣」がよいと思う。
3列以上のデータフレーム df があれば、単に plot(df) とすれば、散布図方陣がつくられる。 しかし、ここでの df のうち、年月日、品質情報、均質番号の列は、総あたりからはずしたほうがよい。 そこで、subset機能をつかって df から観測値のはいった列だけをとりだして plot にわたすことにする。
plot(subset(df, select=c(temp, tempmax, tempmin, precip, sundur)))
結果の図はつぎのようになる。(軸の目盛りなどの調整はまだしていない。)
データの頻度分布(度数分布)をみるのには、ヒストグラム(histogram)が基本的な方法である。 観測値のヒストグラムは、観測値を階級わけして集計し、 階級ごとの頻度を長方形の面積であらわす柱をならべた形のグラフにしたものである。
R には ヒストグラムをつくる関数「hist」がある。ひとまずこれをすなおに実行してみる。
hist(df$sundur)
縦軸のラベルが「Frequency」となっているが、これは「頻度」である (ここでは「振動数」や「周波数」ではない)。
hist にはさまざまなオプションがあって、データ値の階級のわけかたや、 グラフの縦軸・横軸の尺度を変えることができる。
つぎのように hist に2つめの引き数を指定すると、階級わけの区間数が概算でその値になる。
hist(df$sundur, 20)
【hist 関数をつかわずに、集計と作図を別々の手つづきによって、ヒストグラムをつくることもできる。 実例をまだ用意できていないが、できたら追加したい。】
ヒストグラムの情報を要約した表現として 「箱ひげ図」(box and whisker plot)がある。 たとえば、つぎのようにして、日平均気温、最高気温、最低気温の箱ひげ図をならべてつくることができる。
boxplot(df$temp, df$tempmax, df$tempmin)
箱の上下の端はデータの四分位値である。データ値の半分は箱の中にふくまれることになる。 中央値のところに横線がひかれている。 「ひげ」の上下の端はかならずしもデータの最大値・最小値ではなく、 「はずれ値」を「ひげ」からはずして個別に表示することがある。 ひとことで書けないので、Rのドキュメントの boxplot のところを参照してほしい。
【なお、中央値や四分位値の数値を数字で見たいときは、 R の 基本統計機能のうちにある summary 関数を たとえば summary(df$temp) のような形でつかうのが便利である。 最大値・最小値・平均値もあわせて表示される。】
1年を季節にわけることは、それ自体、研究課題になりうる問題でもある。 しかし、多数の地点の観測値をいっしょにあつかうばあいは、 各地点ごとに季節区分を考えるわけにいかず、 同じ季節区分をいっせいにあてはめたい。
大気上端にはいってくる太陽放射は、中緯度では、夏至のころに最大、冬至のころに最小になる。 地上の日射量や日照時間も、もし雲の影響が一定ならば、夏至のころに最大、冬至のころに最小になる。 したがって、夏至・冬至などを基準とした季節のわけかたも、有意義だ。
しかし、気象観測値は、グレゴリオ暦の月ごとに整理する慣習ができている。 そこで、気象学で世界の気候・天候をあつかうときは、 1年を、12・1・2月、3・4・5月、6・7・8月、9・10・11月にわけて「季節」(英語では season)とみなすことが多い。 北半球の中高緯度の感覚で、それぞれ、冬、春、夏、秋とよぶこともある。 しかし、季節名は南半球では逆になるし、熱帯では不適切なこともあるので、 英語の月の名の頭文字をとって、DJF, MAM, JJA, SON とすることも多い。
ここでは、CSVファイルをつくる際に日付を年、月、日にわけておいたから、 その「月」のところを読みとって、年は問わずに、「季節」ごとにまとめてヒストグラムをつくってみよう。
Rでは、subset という関数をつかって、データフレームの部分集合をつくることができる。 「subset(df, 条件)」で、条件にあった行だけが抽出される。 「subset(...)$列名」の形で、その列を参照できる。 条件の指定には、「%in%」(このように「%」の字をふくめる)という演算子をつかって、 「集合にふくまれる」という指定をすることができる。 これによって、たとえば、つぎのように、日照時間のデータを「季節」ごとにわけることができる。 つぎの各行の文字列をそれぞれ hist( ) の中にいれれば、季節別のヒストグラムをかくことができる。
subset(df, month %in% c(12,1,2) )$sundur
subset(df, month %in% c(3,4,5) )$sundur
subset(df, month %in% c(6,7,8) )$sundur
subset(df, month %in% c(9,10,11) )$sundur
ひとつのグラフィックス画面に複数のグラフをならべることは、 あらかじめつぎのように par で mfrow を指定しておけばできる。 つぎの例文のように mfrow を指定したあとには、 6つの plot を、横に2列、縦に3段に配置することができる。
par(mfrow=c(3,2))
2段×2列ならばつぎのようにすればよい。
par(mfrow=c(2,2))
2段×2列にして、四季の日照時間のヒストグラムをならべるプログラムをつくった。
【図の表題には、地点名と観測期間を入れた。 ただし、「DJF 2019-2020」は不正確な表現で、 「2018年12月・2019年1月・2月と2019年12月・2020年1月・2月」と具体的に書くか、 「DJFの期間の表示は1月の属する年で代表させている」のようにことわっておくべきである。】
4つの季節の箱ひげ図をまとめてしめすこともできる。(つぎの例文は長い1行である)。
boxplot(subset(df,month %in% c(12,1,2))$sundur, subset(df, month %in% c(3,4,5))$sundur, subset(df, month %in% c(6,7,8))$sundur, subset(df, month %in% c(9,10,11))$sundur, names=c("DJF", "MAM", "JJA", "SON"), las=1, main="日照時間 (時), 2018-12..2019-11, 館林")
気象要素の種類によって、ヒストグラムの形がちがう。 気温は左右対称でわりあい正規分布に近いが、 降水量は値の大きいほうに長いすそをひく。 地上観測の雲量 (全天のうち雲でおおわれた部分のわりあい)は、0と1の両端の頻度が高い。 (ただし、平均期間を長くしていくと、正規分布にちかづく。) 日照時間も、雲量ほど極端ではないが、2つの山をもつ分布をしている。
藤部(2014)の本の図6.7に日ごとのいくつかの気象要素のヒストグラムがある。 【このうち東京の日降水量のヒストグラムで「40 mm」のすぐ右に 背の高い棒があるのは、40mm以上の値をひとつの柱にまとめてしまったための見かけのもので、 ここに頻度分布の山があるわけではない。】 図6.8には、月ごとの値のヒストグラムがある。 同じ図は、藤部(2015)の講義資料の5ページにもある。 (この講義資料は気象につかわれる統計学の専門的内容もふくんでいるが、 わかるところだけ読んでも有用だと思う。)
離れた地点について、同じ気象要素の値を作図して比較してみよう。
世界の気候の分布にとってはこまかい特徴だが、 日本の冬には、冬の季節風による雲・降水が多い地方と、少ない地方がある。 館林は冬の雲・降水が少ない関東地方の特徴をしめす地点である。 ここからわりあい小さな距離でありながら冬の雲・降水が多いところとして、 新潟県の地点をどこかえらんで、冬の降水量や日照時間を比較してみよう。
複数の地点の同時の観測データによって、 気象要素の空間分布を見ることをふくめた実習を予定しています。
2地点のデータを、期間をそろえていっしょにあつかい、 同時の観測値について、散布図をつくったり、 相関係数を計算したりすることも、別の教材ページとして追加しようと思っています。 ただし、学期中に追加できるかどうか、未定です。