日本の年平均気温図を作ってみた話

この記事は[R Advent Calendar 2020](https://qiita.com/advent-calendar/2020/rlang)の7日目です。筆者の[Japan.R 2020](https://japanr.connpass.com/event/194867/)での[LT](https://ytknzw.github.io/TokyoR/%E6%97%A5%E6%9C%AC%E3%81%AE%E5%B9%B4%E5%B9%B3%E5%9D%87%E6%B0%97%E6%B8%A9%E5%9B%B3%E3%82%92%E4%BD%9C%E3%81%A3%E3%81%A6%E3%81%BF%E3%81%9F%E8%A9%B1_JapanR_2020_20201205.pdf)の内容を補完するもので、時間とプレゼンテーションの兼ね合いで割愛した部分（特にコード）を記載します。

## TL;DR
地理データを用いて、線形回帰により日本の1kmメッシュの年平均気温を推定しました。モデルの精度について懸念点はあるものの、全体としては違和感のない結果を得られました。

※結果だけ見たい方は[こちら](#年平均気温の推定結果)。

## データ
- 気象データ
    - 気象庁：[過去の気象データ・ダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)
    - 参考：[アメダス観測所一覧](http://washitake.com/weather/amedas/obs_stations.md)
        - このサイトの記述をもとに観測地点のリストを作成しました。
- 地理データ
    - 国土交通省（国土数値情報）：[標高・傾斜度3次メッシュデータ](http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-G04-a.html)
    - 国土交通省（国土数値情報）：[土地利用3次メッシュデータ](http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-a.html)
- いずれもオープンデータとは程遠い状態[^1]ですので、入手と加工に手間が掛かります。
- これらのデータについて、[Tokyo.R](https://tokyor.connpass.com/)にてLTをしました。
    - [国土数値情報APIのススメ](https://ytknzw.github.io/TokyoR/%E5%9B%BD%E5%9C%9F%E6%95%B0%E5%80%A4%E6%83%85%E5%A0%B1API%E3%81%AE%E3%82%B9%E3%82%B9%E3%83%A1_TokyoR_20200404.pdf)[^1]
    - [日本の気象観測データの（オープンなようでオープンでない）話](https://ytknzw.github.io/TokyoR/%E6%97%A5%E6%9C%AC%E3%81%AE%E6%B0%97%E8%B1%A1%E8%A6%B3%E6%B8%AC%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%EF%BC%88%E3%82%AA%E3%83%BC%E3%83%97%E3%83%B3%E3%81%AA%E3%82%88%E3%81%86%E3%81%A7%E3%82%AA%E3%83%BC%E3%83%97%E3%83%B3%E3%81%A7%E3%81%AA%E3%81%84%EF%BC%89%E8%A9%B1_TokyoR_20200627.pdf)

### 気象データの緯度、経度のメッシュコードへの変換
`jpmesh`ライブラリーの`coords_to_mesh`関数を使用すると、緯度と経度の組み合わせをそれが表す地点を含むメッシュのコードに変換できます。緯度と経度は`35.69167`や`139.75`のように度を単位とする数値として与えます（**経度が第1引数**です）。

```R
jpmesh::coords_to_mesh(経度, 緯度)
```

### 地理データのシェープファイル（拡張子`shp`）の解凍、読込、縦結合
国土数値情報の地理データは約80km四方の1次メッシュごとにzipファイルが用意されています。`unzip`関数で解凍し、`sf`ライブラリーの`read_sf`関数を使ってシェープファイルをdataframeとして読み込みます。
例えば、標高・傾斜度3次メッシュデータの場合、1次メッシュコード`3036`のデータのzipファイル名は`G04-a-11_3036-jgd_GML.zip`です[^2]。その中にある`G04-a-11_3036-jgd_ElevationAndSlopeAngleTertiaryMesh.shp`が目的のシェープファイルです。

```R
f = "G04-a-11_3036-jgd_GML"
unzip(paste0(f, ".zip"), exdir = f)
shapefile_path <- paste(f,
 "G04-a-11_3036-jgd_ElevationAndSlopeAngleTertiaryMesh.shp", sep = "/")
map <- sf::read_sf(shapefile_path, as_tibble = FALSE)

```
地理データの作成年度によっては、シェープファイルを読み込む際に以下のようなエラーが出ることがあります。

```R
 make.names(vnames, unique = TRUE) でエラー: 
   '<83><81><83>b<83>V<83><85>' に不正なマルチバイト文字があります
```
CP932でエンコードされていることが原因のようですので、`options=c("ENCODING=CP932")`を追加してみて下さい[^3]。

```R
map <- sf::read_sf(shapefile_path, as_tibble = FALSE, options=c("ENCODING=CP932"))
```
最後に読み込んだシェープファイルを、これまでに読み込んだもの（ここでは`map_all`とします）と縦結合する必要があります。`rbind`でもできるのですが処理が遅いので、`data.table`ライブラリーの`rbindlist`関数と`sf`ライブラリーの`st_as_sf`関数の組み合わせをお勧めします[^4]。`map_all`と`map`をlistとして渡します。

```R
map_all <- data.table::rbindlist(list(map_all, map)) %>% sf::st_as_sf()
```

## モデル
### 変数
- 目的変数：年平均気温（摂氏）
- 説明変数：
    - 1kmメッシュの中心の緯度と経度
    - 地形
        - 最高標高、最低標高、平均標高
        - 最大傾斜角度、最小傾斜角度、平均傾斜角度
        - 最大傾斜方向、最小傾斜方向
    - 土地利用種別（面積）
        - 田、その他の農用地
        - 森林、荒地、建物用地、道路、鉄道、その他の用地
        - 河川地及び湖沼、海浜、海水域
        - ゴルフ場

### 地理データのメッシュコードから緯度、経度の取得
メッシュの中心の緯度と経度は`jpmesh`ライブラリーの`mesh_to_coords`関数を使うと得られます。メッシュコードは文字列として与えます。

```R
jpmesh::mesh_to_coords(メッシュコード)
```

### 傾斜方向のOne-hot encoding
傾斜方向は、8方位が北を1として時計回りに1から7までの数字、方向なしを8として表されるカテゴリー変数です。これを`recipes`ライブラリーの`step_dummy`関数を用いてOne-hot encoding[^5]してダミー変数にします。最大傾斜方向と最小傾斜方向はそれぞれ`G04a_007`と`G04a_009`というカラムに格納されています。

```R
# （前後の処理は省略）
recipes::step_dummy(G04a_007) %>%
recipes::step_dummy(G04a_009) %>%
```

### Cross validation（交差検証）による性能の推定
気象データ（2019年の観測地点ごとの年平均気温）と地理データが全て紐づく928件（ここでは`learn`とします）をモデル構築に使用し、線形回帰モデルを採用しました。線形回帰の性能を`rsample`ライブラリーの`vfold_cv`関数を用いてCVにより推定しました[^6]。ここでは、10-fold CVを10回繰り返し、メッシュの中心緯度を4等分してstratify[^7]しています。使用した指標はRMSE、MAE、MAPE、R<sup>2</sup>、CCC[^8]です。目的変数の年平均気温は`avg_temp`というカラムに入っています。

```R
# CVの設定
set.seed(123)
cv <- learn %>% rsample::vfold_cv(v = 10, repeats = 10,
 strata = "lat_cent", breaks = 4)

# 線形回帰モデルの用意
lr_mod <- linear_reg(mode = "regression") %>% set_engine(engine = "lm")

# CVの中で使う関数の定義
cv_fit <- function(splits, mod, ...) {
  res_mod <-
    fit(mod, avg_temp ~ ., data = analysis(splits), family = gaussian)
  return(res_mod)
}
cv_pred <- function(splits, mod){
  holdout <- assessment(splits)
  pred_assess <- bind_cols(truth = holdout$avg_temp,
                           predict(mod, new_data = holdout))
  return(pred_assess)
}

# CVの実行
res_cv_train <- cv %>% 
  mutate(res_mod = map(splits, .f = cv_fit, lr_mod),
         res_pred = map2(splits, res_mod, .f = cv_pred))

# 指標の指定
multi_metric <- yardstick::metric_set(yardstick::rmse,
                                      yardstick::mae,
                                      yardstick::mape, rsq, ccc)
# 指標の計算
res_metric <- res_cv_train %>%
  mutate(metrics = map(res_pred, multi_metric, truth = truth, estimate = .pred)) %>% 
  unnest(metrics)
```

結果は以下の通りですが、線形回帰と簡単な説明変数だけにしては良すぎる結果なので、今後精査します。

|指標|平均|標準偏差|
|----|----|----|
|RMSE|0.705|0.0595|
|MAE|0.546|0.0431|
|MAPE|4.89|0.472|
|R<sup>2</sup>|0.973|0.00589|
|CCC|0.975|0.00306|

### モデル構築
前述の懸念点はあるものの、線形回帰を用いて年平均気温を推定するモデルを作れることは分かったので、実際にモデルを作って推定してみます。推定対象の全データ（構築用データ含む）をここでは`dat`とし、推定結果を`pred`というカラムに入れます。

```R
mdl <- linear_reg(mode = "regression") %>%
  set_engine(engine = "lm") %>%
  fit(avg_temp ~ ., data = learn, family = gaussian)
dat <- dat %>% mutate(pred = predict(mdl, .)$.pred)
```

## 年平均気温の推定結果
`geom_pred`を用いてコロプレス図[^9]を描画します。図形の元になるポリゴンデータは`geometry`カラムに入っています。各メッシュの年平均気温の推定値`pred`に応じて色を塗ります。

```R
ggplot(dat) +
  geom_sf(aes(fill = pred, geometry = geometry), colour = NA) +
  labs(title = "Estimated Avg Temperature") +
  coord_sf() +
  scale_fill_gradientn(colours = rev(rainbow(5)),
                       limits = c(-10, 35),
                       name = "Estimated\navg temp (C)")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412319/e82c18ec-6389-ea9d-fa87-817da656c470.png)
日本全体で見ると、概ね経験に近い結果です。一方で、海流や季節風、日照など気候の影響が反映されていないようです（地形や土地利用といった地理データしか使っていないので当然と言えば当然ですが）。

次に、関西・中部と関東を拡大してみます。分かりやすいようにグラデーションを変更しました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412319/54c36648-5d8c-6ef9-0fb3-a44d04c30dc9.png)
京都、大阪、奈良、名古屋といった都市部や平野部と琵琶湖が周辺より年平均気温が高いと推定されていて、違和感はありません。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/412319/7519db72-085f-8bac-d9fc-caff80613a37.png)
それに比べると関東は、平野部（＝都市部）が一様に高いと推定されています。一方で、箱根や丹沢、群馬県北西部の山地は低いと推定されていて、全体としては違和感はありません。

## まとめ

地理データと線形回帰から、日本の1kmメッシュの年平均気温を推定することができました。ただし、モデルの精度が良すぎるように思われるので調査が必要です。一方で、今後は、他の推定手法を試したり、説明変数（データ）を追加したり、時系列予測をしてみたいと思います。結果はTokyo.RのLTにて発表していきます。

Enjoy!

[^1]: 国土数値情報については、2020年6月までは「国土数値情報API」という公式のAPIサービスがベータ版という位置付けながら提供されていました。
[^2]: 「11」は、今回使用したデータが作成された2011年度を指すようです（測量自体は2009年度）。
[^3]: https://ksmzn.hatenablog.com/entry/sf_mapview_001
[^4]: https://github.com/r-spatial/sf/issues/798#issuecomment-405171212
[^5]: ここでは、線形独立のためにカテゴリーが1つ除かれます。
[^6]: https://stateofther.github.io/finistR2019/s-tidymodels.html
[^7]: 日本語では「層化抽出法」と呼ばれ、指定した変数の分布を分割データ（サンプル）間で同じにする手法です。
[^8]: Concordance correlation coefficient
[^9]: 正確な地図ではないですが、色を塗り分けるのでコロプレス図と呼んでもよいと思います。