More than 3 years have passed since last update.

RAdvent Calendar 2020

@yutakanzawa(Yuta Kanzawa)

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

Last updated at 2020-12-07Posted at 2020-12-07

この記事はR Advent Calendar 2020の7日目です。筆者のJapan.R 2020でのLTの内容を補完するもので、時間とプレゼンテーションの兼ね合いで割愛した部分（特にコード）を記載します。

TL;DR

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

※結果だけ見たい方はこちら。

データ

気象データ
- 気象庁：過去の気象データ・ダウンロード
- 参考：アメダス観測所一覧
  - このサイトの記述をもとに観測地点のリストを作成しました。
地理データ
- 国土交通省（国土数値情報）：標高・傾斜度3次メッシュデータ
- 国土交通省（国土数値情報）：土地利用3次メッシュデータ
いずれもオープンデータとは程遠い状態¹ですので、入手と加工に手間が掛かります。
これらのデータについて、Tokyo.RにてLTをしました。
- 国土数値情報APIのススメ¹
- 日本の気象観測データの（オープンなようでオープンでない）話

気象データの緯度、経度のメッシュコードへの変換

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

jpmesh::coords_to_mesh(経度, 緯度)

地理データのシェープファイル（拡張子`shp`）の解凍、読込、縦結合

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

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)

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

 make.names(vnames, unique = TRUE) でエラー: 
   '<83><81><83>b<83>V<83><85>' に不正なマルチバイト文字があります

CP932でエンコードされていることが原因のようですので、options=c("ENCODING=CP932")を追加してみて下さい³。

map <- sf::read_sf(shapefile_path, as_tibble = FALSE, options=c("ENCODING=CP932"))

最後に読み込んだシェープファイルを、これまでに読み込んだもの（ここではmap_allとします）と縦結合する必要があります。rbindでもできるのですが処理が遅いので、data.tableライブラリーのrbindlist関数とsfライブラリーのst_as_sf関数の組み合わせをお勧めします⁴。map_allとmapをlistとして渡します。

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

モデル

変数

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

地理データのメッシュコードから緯度、経度の取得

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

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

傾斜方向のOne-hot encoding

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

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

Cross validation（交差検証）による性能の推定

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

# 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²	0.973	0.00589
CCC	0.975	0.00306

モデル構築

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

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を用いてコロプレス図⁹を描画します。図形の元になるポリゴンデータはgeometryカラムに入っています。各メッシュの年平均気温の推定値predに応じて色を塗ります。

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)")

日本全体で見ると、概ね経験に近い結果です。一方で、海流や季節風、日照など気候の影響が反映されていないようです（地形や土地利用といった地理データしか使っていないので当然と言えば当然ですが）。

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

京都、大阪、奈良、名古屋といった都市部や平野部と琵琶湖が周辺より年平均気温が高いと推定されていて、違和感はありません。

それに比べると関東は、平野部（＝都市部）が一様に高いと推定されています。一方で、箱根や丹沢、群馬県北西部の山地は低いと推定されていて、全体としては違和感はありません。

まとめ

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

Enjoy!

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

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up