ケッペンの気候区分を多変量解析で検証

# 目次
0 はじめに
1 設定した問題 (多変量解析によってケッペンの気候区分を再現できるかの検証)
　　1.1 ケッペンの気候区分とは
2 データの出典
　　2.1 高度、気温と降水量のデータ(平年値)
　　2.2 各都市とそれぞれの都市が属するケッペンの気候区分について書かれたもの
3 仮説
4 結果と考察
　　4.1 クラスター分析
　　　　4.1.1 各月の気温と降水量のデータを使ったとき
　　　　4.1.2 高度、各月の気温と降水量のデータを使ったとき
　　4.2 主成分分析
　　4.3 判別分析
5 結論

# 0 はじめに
　2019年1月に大学の講義の期末レポートで、Rで多変量解析する自由課題を行った。使ったファイルは以下の通り。https://github.com/m-kunugi/R/tree/master/koeppen

• koeppen.r(解析に使った R スクリプト)
• world_temp_prec.csv(データ)

以下はデータの一部

| 都市名 | class | height | t_year | p_year | t1 | t2 | … | t12 | p1 | p2 | … | p12 |
|:--|:--|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|
| Afシンガポール（シンガポール） | A | 5 | 27.6 | 2199 | 26.6 | 27.2 | … | 26.6 | 246.3 | 114.1 | … | 314.8 |
| Amポートハーコート(ナイジェリア) | A | 18 | 26.5 | 2258.8 | 26.8 | 27.7 | … | 26.2 | 22.6 | 65.7 | … | 25.2 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
>データ各列の説明
1 列目→都市名(先頭のアルファベットは実際のケッペンの気候区分)、
2 列目→気候区分(ABCDEH)、3 列目→高度(m)、
4 列目→年平均気温(°C)、5 列目→年降水量(mm)、
6~17 列目(t1~t12)→各月平均気温(°C)、
18~29 列目(p1~p12)→各月降水量(mm)

#1 設定した問題 (多変量解析によってケッペンの気候区分を再現できるかの検証)
　日本国内の気候区分を、多変量解析によって詳細に考察する研究は既にあった[^1]。しかし、世界各地の気候区分を多変量解析したものは見当たらなかった(少なくとも日本語文献では)。今回は世界の各都市の気温、降水量 (、高度) のデータを用いて多変量解析を行い、それにより導かれる気候区分と、有名なケッペンの気候区分(高校地理で習う)とを比較してみた。その上で多変量解析によって、ケッペンの気候区分を（この気候区分が正しいと仮定して）どれだけ再現できるかどうかを検証してみることにする。
[^1]: 気候要素の変動パターンに基づく日本の気候区分
  https://www.chs.nihon-u.ac.jp/institute/nature/kiyou/2012/pdf/2_6.pdf

##1.1 ケッペンの気候区分とは
　ケッペンの気候区分は、気温と降水量の2変数から経験的な判定法により決定される気候区分である。植生の 特徴を反映しているのが長所である。
参考として具体的な判定法は下の参考1に載せた。

<details><summary>参考1: ケッペンの気候区分の判定法</summary><div>
>
(1)植物分布を基礎に区分
　(ア)樹木気候・・・熱帯気候(A)、温帯気候(C)、冷帯気候(D)
　(イ)無樹木気候・・・乾燥帯気候(B)、寒帯気候(E)
(2)気温と降水量を組み合わせて12種類に気候区分
　(ア)気温
　　1　熱帯気候(A)・・・最寒月平均気温18°C以上
　　2　温帯気候(C)・・・最寒月平均気温-3°C以上、18°C未満かつ最暖月が10°C以上
　　3　冷帯気候(D)・・・最寒月平均気温-3°C未満かつ最暖月平均気温10°C以上
　　4　寒帯気候(E)・・・最暖月平均気温10°C未満
　(イ)降水量
　　1 f ・・・一年中降雨
　　2 w ・・・冬季少雨
　　3 s ・・・夏季少雨
　(ウ)A の区分
　　Af(熱帯雨林気候)・・・最少雨月降水量60 mm 以上
　　Am(弱い乾季のある熱帯雨林気候)・・・最少雨月降水量≧100-0.04×年降水量
　　Aw(サバナ気候)・・・最少雨月降水量<100-0.04×年降水量
　(エ)C の区分
　　Cs(地中海性気候)・・・冬の最多雨月降水量が夏の最少雨月降水量の3倍以上かつ最少雨月                降水量が30 mm 未満
　　Cw(温暖冬季少雨気候)・・・夏の最多雨月降水量が冬の最少雨月降水量の10倍以上
　　Cf・・・Cs、Cw でない
　　　Cfa(温暖湿潤気候)・・・最暖月平均気温22°C以上
　　　Cfb(西岸海洋性気候)・・・最暖月平均気温22°C未満
　(オ)D の区分
　　Df(冷帯湿潤気候)・・・夏の最多雨月降水量が冬の最少雨月降水量の10倍未満
　　Dw(冷帯冬季少雨気候)・・・夏の最多雨月降水量が冬の最少雨月降水量の10倍以上
　(カ)E の区分
　　ET(ツンドラ気候)・・・最暖月平均気温0°C以上10°C未満
　　EF(氷雪気候)・・・最暖月平均気温0°C未満
　(キ)B の区分・・・乾燥限界により区分
　　w型(冬季乾燥)・・・乾燥限界=20×(年平均気温+14)
　　f型(年中湿潤)・・・乾燥限界=20×(年平均気温+7)
　　s型(夏季乾燥)・・・乾燥限界=20×年平均気温
　　BS(ステップ気候)・・・乾燥限界の 1/2 ≦年降水量<乾燥限界
　　BW(砂漠気候)・・・年降水量<乾燥限界の 1/2  
2 高山気候(H)
　　(1)ケッペンの気候区分にはない→ケッペンの気候区分を補正し設定
　　(2)海抜高度の高い地域では、気温の逓減率により同緯度の低地より低温となる。
　　　熱帯地方の高山(標高2000~3000m)→常春の気候

参考 1: ケッペンの気候区分の判定法

- ケッペンの気候区分 http://azeta.jp/tiri/4-2/kikouku.html

より転載。定義が紛らわしい部分については、以下のWikipediaを参考に修正を加えた。

- Wikipedia ケッペンの気候区分
https://ja.wikipedia.org/wiki/ケッペンの気候区分
</div></details>

#2 データの出典
　多変量解析をするにあたって、各都市の高度、気温と降水量、そして結果の確認のためにその都市が属するケッペンの気候区分のデータが必要となる。

##2.1 高度、気温と降水量のデータ(平年値)
- データブック オブ・ザ・ワールド  2015 年版 二宮書店

　上記のデータ(高度、気温、降水量)を全て打ち込むのは大変なので以下の気象庁の全く同じデータを適宜コピーしてCSVファイルを作成した。**また、北半球と南半球では季節が逆転する為、北半球のデータのみにした。**

- 気象庁 主な地点の平年値  
https://www.data.jma.go.jp/gmd/cpd/monitor/mainstn/nrmlist.php

##2.2 各都市とそれぞれの都市が属するケッペンの気候区分について書かれたもの
- データブック オブ・ザ・ワールド  2015 年版 二宮書店

　こちらも上記の文献を参考にした。

#3 仮説
　クラスター分析を行えば、おおよそケッペンの気候区分の ABCDE の区分の通りに分類できると予想される。 しかし各区分内のさらに細かい分類(例えば Am と Aw など)は流石にケッペンの気候区分とはあまり一致しないであろうと考えられる。

#4 結果と考察
　クラスター分析、主成分分析、判別分析をそれぞれ試みた。

##4.1 クラスター分析
　ユークリッド距離[^2]を用いた階層的クラスタリングでクラスター分析を行った結果を以下に記す。また k-平均法を使った非階層的クラスタリングも試してみたが、階層的クラスタリングと比べて上手くケッペンの気候区分に一致しなかった。
[^2]: ユークリッド距離
  二点間の「通常の」距離。例えば、点 $ \mathbf{p} $ および $ \mathbf{q} $ の間のユークリッド距離とは、それらをつなぐ線分${\displaystyle {\overline {\mathbf {pq} }}}$ の長さ$ d(\mathbf{p}, \mathbf{q})=\sqrt{\left(q_{1}-p_{1}\right)^{2}+\left(q_{2}-p_{2}\right)^{2}+\cdots+\left(q_{n}-p_{n}\right)^{2}}=\sqrt{\sum_{i=1}^{n}\left(q_{i}-p_{i}\right)^{2}} $
  のこと。他に用いられる距離の例として、ミンコフスキー距離やマンハッタン距離がある。

###4.1.1 各月の気温と降水量のデータを使ったとき
　まず、各月の気温と降水量のデータを用いてクラスター分析した。それと実際の気候区分とを照らし合わせてみた。
![図 1: クラスター分析(各月の気温と降水量のデータを使ったとき)](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/282333/ced5cc49-c995-bccb-1152-133841baa39d.png)
図 1: クラスター分析(各月の気温と降水量のデータを使ったとき)

　以下はクラスター分析の結果の確認である。5つに分けたのはこの時が最もケッペンの気候区分に近かったためである。
　これを見ると、A、C、D を綺麗に分類するのは難しいことが分かる。一方で B は一部、C と D のグループに 入ってしまっている(マシュハドとカシガル)が、比較的分類できていることが分かる。また、Af と E はかなりケッペンの気候区分に近い分類ができている。前者は赤道付近、後者は北極付近で見られる気候である。
　**以上のことから低緯度と高緯度は分類しやすいが、中緯度の分類は難しいことが見て取れる。**
 
```R:koeppen.r
# ソースコード(10~17 行目)
> ## 結果の確認 (各クラスター内の都市名を表示)
> clust <- cutree(est,5) # 気候区分を 5つに分ける。
> perf <- rownames(mydata) # 都市名の取得
> perf[clust==1] # Af(熱帯雨林気候)
[1] "Af シンガポール(シンガポール)" "Af コロンボ(スリランカ)"
> perf[clust==2] # A(Af 以外)と C
[1] "Am ポートハーコート(ナイジェリア)" "Am マイアミ(アメリカ)"
[3] "Aw バンコク(タイ)" "Aw メリダ(メキシコ)"
[5] "Cw ホンコン(中国)" "Cw グワハティ(インド)"
[7] "Cfa 東京(日本)" "Cfa シャンハイ(中国)"
[9] "Cfa ニューオーリンズ(アメリカ)" "H アディスアベバ(エチオピア)"
> perf[clust==3] # BとCとD
[1] "BS マシュハド(イラン)" "BW カシガル(中国)" "Cs ニース(フランス)"
[4] "Cs サンフランシスコ(アメリカ)" "Cw チンタオ(中国)" "Cfb パリ(フランス)"
[7] "Cfb ロンドン(イギリス)" "Dfa 札幌(日本)" "Dfa シカゴ(アメリカ)"
[10] "Dfb モスクワ(ロシア)" "Dfc アルハンゲリスク(ロシア)" "Dw チャンチュン(中国)"
[13] "Dwb ハバロフスク(ロシア)" "Dwc チタ(ロシア)" "H ラサ(中国)"
> perf[clust==4] # BS と BW
[1] "BS ワガドゥーグー(ブルキナファソ)" "BS ンジャメナ(チャド)"
[3] "BW リヤド(サウジアラビア)" "BW カイロ(エジプト)"
> perf[clust==5] # ET(ツンドラ気候)
[1] "ET ディクソン(ロシア)" "ET バロー(アラスカ)"
```

###4.1.2 高度、各月の気温と降水量のデータを使ったとき
　H(高山気候)は気温や降水量で決定されるものではない。そのため、今度は高度のデータを加えて、同じようなことをしてみた。高山気候が分類できることが予想される。
![図 2: クラスター分析(高度、各月の気温と降水量のデータを使ったとき)](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/282333/498857ae-1336-4527-1e26-5713633f2457.png)
図 2: クラスター分析(高度、各月の気温と降水量のデータを使ったとき)

　以下はクラスター分析の結果の確認である。6つに分けたのはこの時が最もケッペンの気候区分に近かったためである。
　これを見ると、先ほどとは違って Af と B と E と H は綺麗に分けられる。中緯度の分類はまだ難しいものの、 A(Af 以外)と C、C と D とに分類はできている。
　**以上のことからかなりケッペンの気候区分に近い分類ができていると分かる。**
　ケッペンの気候区分では高山気候の分類にしか使わない、高度のデータも使うと精度が上がるのが意外だった。
 
```R:koeppen.r
# ソースコード(25~33 行目)
> ## 結果の確認 (各クラスター内の都市名を表示)
> clust <- cutree(est,6) # 気候区分を 6つに分ける。
> perf <- rownames(mydata) # 都市名の取得
> perf[clust==1] # Af(熱帯雨林気候)
[1] "Af シンガポール(シンガポール)" "Af コロンボ(スリランカ)"
> perf[clust==2] # A(Af 以外)と C
[1] "Am ポートハーコート(ナイジェリア)" "Am マイアミ(アメリカ)" "Aw バンコク(タイ)"
[4] "Aw メリダ(メキシコ)" "Cw ホンコン(中国)" "Cw グワハティ(インド)"
[7] "Cfa 東京(日本)" "Cfa シャンハイ(中国)" "Cfa ニューオーリンズ(アメリカ)"
> perf[clust==3] # B(乾燥帯気候)
[1] "BS マシュハド(イラン)" "BS ワガドゥーグー(ブルキナファソ)"
[3] "BS ンジャメナ(チャド)" "BW カシガル(中国)"
[5] "BW リヤド(サウジアラビア)" "BW カイロ(エジプト)"
> perf[clust==4] # C と D
[1] "Cs ニース(フランス)" "Cs サンフランシスコ(アメリカ)" "Cw チンタオ(中国)"
[4] "Cfb パリ(フランス)" "Cfb ロンドン(イギリス)" "Dfa 札幌(日本)"
[7] "Dfa シカゴ(アメリカ)" "Dfb モスクワ(ロシア)" "Dfc アルハンゲリスク(ロシア)"
[10] "Dw チャンチュン(中国)" "Dwb ハバロフスク(ロシア)" "Dwc チタ(ロシア)"
 > perf[clust==5] # ET(ツンドラ気候)
[1] "ET ディクソン(ロシア)" "ET バロー(アラスカ)" > perf[clust==6] # H(高山気候)
[1] "H ラサ(中国)" "H アディスアベバ(エチオピア)"
```

　他にも様々なデータの組み合わせ(例えば、高度と年平均気温、年降水量のみで解析など)を試してみた。しかし、高度、各月の気温と降水量のデータを使った今回がケッペンの気候区分に一番適合しそうなので、この組み合わせで次に主成分分析を行ってみることにする。

##4.2 主成分分析
　高度、各月の気温と降水量のデータを使って、主成分分析を行った。横軸に第一主成分、縦軸に第二主成分をプロットした図が以下の通り。heightは高度、t1~12 は 1~12 月の気温、p1~12 は 1~12 月の降水量をそれぞれ表している。
　第一主成分は左にシンガポールなどの低緯度地域、右にロシアやアラスカなどの高緯度地域が集まっていることを考えると、**緯度の高さ(あるいは年平均温度の低さ、年降水量の低さ)**を示していると考えられる。また、第二主成分はこの図だけでは断定しにくいが、**冬の降水量の多さ**を示していると考えられる。
![図 3: 主成分分析(第一主成分と第二主成分)](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/282333/cb9c38dd-14a6-c648-a898-3020a587cd55.png)
図 3: 主成分分析(第一主成分と第二主成分)

　さらに、横軸に第二主成分、縦軸に第三主成分をプロットした図が以下の通り。第三主成分は、**夏の降水量の多さ**を示していると考えられる。
![図 4: 主成分分析(第二主成分と第三主成分)](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/282333/a5bce2e5-8c02-0dfa-c531-10cbc0fda152.png)
図 4: 主成分分析(第二主成分と第三主成分)

##4.3 判別分析
　E と H はデータが2個ずつしかないので A、B、C、D のみで判別分析を行うことにする。**しかし、A、B、C、D もそれぞれ数個ずつしかデータが無かったので、判別分析は上手くいかなかった。線形判別は過学習し、データ数が足りないため2次判別はそもそも出来なかった。**

```R:koeppen.r
# ソースコード(56~57 行目)
> table(true=abcd.train$class,est=predict(myld)$class) # 訓練誤差 est
true A B C D
   A 1 0 0 0
   B 0 3 0 0
   C 0 0 5 0
   D 0 0 0 5
> table(true=abcd.train$class,est=myldloo$class) # CV による予測誤差 est
true A B C D
   A 0 0 0 0
   B 2 0 0 1
   C 0 1 1 3
   D 2 1 1 1
> ## 線形判別は過学習している
```

　手軽に使える気候区分のデータをこれ以上集められそうに無かったので、これ以上の分析はできなかった。しかし、南半球のデータを季節を反転させて、北半球に直すなどすれば十分なデータ数に増やすことが可能かと考えられる。

###補足　南半球のデータを季節逆転して追加
時間があれば加筆。

#5 結論
- クラスター分析では、高度、各月の気温と降水量のデータを使ったとき、ケッペンの気候区分に近い分類ができていると分かった。ケッペンの気候区分では高山気候の分類にしか使わない、高度のデータも使うと上手くいくのが意外だった。
- 主成分分析では、第一主成分は緯度の高さ(あるいは年平均温度の低さ、年降水量の低さ)、第二主成分は冬の降水量の多さ、第三主成分は夏の降水量の多さをそれぞれ示していると考えられる。
- 判別分析では、A、B、C、D もそれぞれ数個ずつしかデータが無かったので、上手くいかなかった。線形判別は過学習し、データ数が足りないため2次判別はそもそも出来なかった。