Rで解析：地図の描写が面白い「tmap」パッケージ

Rの解析に役に立つ記事
2022.03.06
スポンサーリンク

ggplot2パッケージと同様にレイヤーの概念で数量データを地図に簡単に示すことが可能なパッケージの紹介です。また、分割基準を設定することで地図を簡単に分割して表示することが可能です。

付与するデータ次第で面白い表現が可能と考えます。

パッケージバージョンは3.3-3。実行コマンドはwindows 11のR version 4.1.2で確認しています。

スポンサーリンク

「tmap」パッケージ利用の前準備

・シェープファイルを準備する

「JPNPref.shp」ファイルに気象庁より観測所一覧を取得して「JPNPref_2.shp」ファイルを作成します。

・Rで解析：「leaflet」パッケージで使える日本のシェープファイルの作成例
　https://www.karada-good.net/analyticsr/r-594

コマンド実行に必要な「JPNPref.shp」ダウンロード
コマンド実行で作成できる「JPNPref_2.shp」ダウンロード
#パッケージの読み込み
library("tcltk")

#「JPNPref.shp」シェープファイルを読み込む
#参考:https://www.karada-good.net/analyticsr/r-594
load(paste0(as.character(tkgetOpenFile(title = "ファイルを選択",
                                       filetypes = '{"ファイル" {".*"}}',
                                       initialfile = c("*.*")))))

###気象庁より観測所一覧を取得#####
#Tmpフォルダに保存しデスクトップにファイル展開後に読み込む
#ファイル名が変更されため下記URLよりリンクを確認する
#https://www.jma.go.jp/jma/kishou/know/amedas/kaisetsu.html
temp <- tempfile()
download.file("https://www.jma.go.jp/jma/kishou/know/amedas/ame_master_20211207.zip",
              temp, mode = "wb")
Kansokujyo <- read.csv(unzip(temp))
unlink(temp)

###必要データを日本地図シェープファイルに付与して保存#####
MasterKansokujyo <- Kansokujyo[, c(1, 4, 7:10)]

#北海道の各地方表記を北海道に置換
MasterKansokujyo[, 1] <- c(rep("北海道", 227),
                           as.character(MasterKansokujyo[228:nrow(MasterKansokujyo), 1]))

#緯度経度を変換
LatData <- type.convert(paste0(MasterKansokujyo[, 3] + MasterKansokujyo[, 4]/60),
                        as.is = TRUE)
LonData <- type.convert(paste0(MasterKansokujyo[, 5] + MasterKansokujyo[, 6]/60),
                        as.is = TRUE)

#シェープファイルに付与
#文字エンコードをUTF-8にする
if(!require("stringi", quietly = TRUE)){
  install.packages("stringi");require("stringi")
}
attributes(JPNPref)$Kansokujyo <- data.frame("都道府県" = stri_encode(MasterKansokujyo[, 1], "", "utf-8"),
                                             "観測所名" = stri_encode(MasterKansokujyo[, 2], "", "utf-8"),
                                             "経度" = LatData, "緯度" = LonData)
#シェープファイルを保存
save(JPNPref, file = "JPNPref_2.shp")

・気象庁より最高気温最新を取得する

コマンドを実行すると前日の最高気温、観測所の緯度経度が付与された「JPNPref_3.shp」と「MaxSpDF.shp」が作成できます。

コマンド実行で作成できる「JPNPref_3.shp」ダウンロード
コマンド実行で作成できる「MaxSpDF.shp」ダウンロード
#パッケージの読み込み
library("tcltk")
if(!require("tidyverse", quietly = TRUE)){
  install.packages("tidyverse");require("tidyverse")
}

#「JPNPref_2.shp」シェープファイルを読み込む
#参考:https://www.karada-good.net/analyticsr/r-594
load(paste0(as.character(tkgetOpenFile(title = "ファイルを選択",
                                       filetypes = '{"ファイル" {".*"}}',
                                       initialfile = c("*.*")))))

###気象庁より最高気温最新を取得#####
#http://www.data.jma.go.jp/obd/stats/data/mdrr/docs/csv_dl_readme.html
MaxTemp <- read.csv("http://www.data.jma.go.jp/obd/stats/data/mdrr/tem_rct/alltable/mxtemsadext00_rct.csv")

#地点データのカッコを削除;カッコに表記ゆれがある
MaxTemp$地点 <- str_replace_all(MaxTemp$地点, pattern = "\\（.+?\\）", replacement = "")
MaxTemp$地点 <- str_replace_all(MaxTemp$地点, pattern = "\\(.+?\\)", replacement = "")
MaxTemp$地点 <- str_replace_all(MaxTemp$地点, pattern = "\\）", replacement = "")

#各都道府県の最高気温を抽出
#データ保管用変数
NewMaxTemp <- data.frame()
#処理
for(n in 1:47){
  #都道府県を抽出
  GetPrefData <- MaxTemp[which(MaxTemp[, 2] %in%
                                 grep(JPNPref@data$name_local[n], MaxTemp[, 2], value = TRUE)),]
  #最高気温を降順で並び替え
  GetPrefData <- GetPrefData[order(GetPrefData[, 10], decreasing = TRUE),]
  #観測所の対象都道府県
  GetKansokujyo <- JPNPref@Kansokujyo[which(JPNPref@Kansokujyo[, 1] %in%
                                              grep(substr(JPNPref@data$name_local[n], 1, 2),
                                                   as.character(JPNPref@Kansokujyo[, 1]), value = TRUE)),]
  #観測所の経度緯度を取得
  LatLonData <- GetKansokujyo[which(GetKansokujyo[, 2] %in% GetPrefData[1, 3]),][1, 4:3]
  #1位の地点名と最高気温,経度緯度を結合
  NewMaxTemp <- rbind(NewMaxTemp, cbind(JPNPref@data$name_local[n], GetPrefData[1, c(3, 10)], LatLonData))}

#1列目と2列目を結合後に1列目を置き換え
NewMaxTemp[, 1] <- paste0(NewMaxTemp[, 1], ":", NewMaxTemp[, 2])

#2列目を削除
NewMaxTemp <- NewMaxTemp[, -2]

#データ名を変更
colnames(NewMaxTemp) <- c("都道府県:地点名", "気温", "lat", "lon")

#JPNPref@dataへ追加
JPNPref@data <- cbind(JPNPref@data, NewMaxTemp[, 1:2])

#JPNPrefファイルを保存
save(JPNPref, file = "JPNPref_3.shp")

#SpatialPointsDataFrameを作成
MaxSpDF <- SpatialPointsDataFrame(SpatialPoints(NewMaxTemp[, 3:4],
                                                proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")),
                                  data = data.frame(id = as.factor(NewMaxTemp[, 1])))

#MaxSpDFファイルを保存
save(MaxSpDF, file = "MaxSpDF.shp")

パッケージのインストール

下記コマンドを実行してください。環境によっては「terra」パッケージのバージョンが古く、「tmap」パッケージがうまく起動しないかもしれません。その場合は「install.packages(“terra”)」を実行してください。

#パッケージのインストール
install.packages("tmap")
#必要に応じて下記実行
install.packages("terra")

コマンドの紹介

詳細はコマンド、パッケージのヘルプを確認してください。また、記事内の「tmap」パッケージ利用の前準備を参考に「JPNPref_3.shp」と「MaxSpDF.shp」を準備してください。

#パッケージの読み込み
library("tcltk")
library("tmap")

#作成した「JPNPref_3.shp」シェープファイルを読み込む
load(paste0(as.character(tkgetOpenFile(title = "ファイルを選択",
                                       filetypes = '{"ファイル" {".*"}}',
                                       initialfile = c("*.*")))))

#作成した「MaxSpDF.shp」シェープファイルを読み込む
load(paste0(as.character(tkgetOpenFile(title = "ファイルを選択",
                                       filetypes = '{"ファイル" {".*"}}',
                                       initialfile = c("*.*")))))

#プロット出力を選択する:tmap_modeコマンド
#プロット作成前に実行することで出力方法を選択することができます
#ggplot2でプロット
#オプションをplotにする
tmap_mode("plot")
#leafletでプロット
#オプションをviewにする
tmap_mode("view")

#簡易プロット:qtmコマンド
qtm(JPNPref, fill = "気温", text = "都道府県:地点名", style = "gray",
    text.root = 3, legend.show = TRUE)

#基本はtm_shapeコマンドにtm_XXXXコマンドを実行して作成する
tm_shape(JPNPref) +
  tm_polygons("気温", id = "都道府県:地点名", legend.show = TRUE) +
  tm_shape(MaxSpDF) +
  tm_bubbles(size = .1, "blue") +
  tm_style("gray")

###例えば気温で分割してプロット###
#気温を区分分け:「fancycutp」パッケージを利用
#https://www.karada-good.net/analyticsr/r-544
if(!require("fancycut", quietly = TRUE)){
  install.packages("fancycut");require("fancycut")
}
#分類
CutLabel <- wafflecut(x = JPNPref@data$気温,
                      intervals = c("[0, 10)", "[10, 15)", "[15, 20)", "[20, 30]"),
                      buckets = c("0-9", "10-14", "15-19", "20-30"),
                      unmatched.bucket = "範囲外")
#結合
JPNPref@data <- cbind(JPNPref@data, CutLabel)

####ggplot2でプロット###
tmap_mode("plot") +
  tm_shape(JPNPref) +
  tm_polygons("気温", id = "都道府県:地点名", legend.show = TRUE) +
  tm_facets(by = "CutLabel") +
  tm_shape(MaxSpDF) +
  tm_bubbles(size = .1, "red") +
  tm_style("gray")

####leafletでプロット###
tmap_mode("view") +
tm_shape(JPNPref) +
  tm_polygons("気温", id = "都道府県:地点名", legend.show = TRUE) +
  #分割:tm_facetsコマンド;
  #syncオプション:TRUEを設定すると各分割が同期します
  tm_facets(by = "CutLabel", sync = TRUE) +
  tm_shape(MaxSpDF) +
  tm_bubbles(size = .1, "red") +
  tm_style("gray")

出力例

・例えば気温で分割してggplot2でプロット

・例えば気温で分割してリーフレットでプロット

拡大・縮小の下にあるレイヤーボタンより表示レイヤーを選択することが可能です。

補足：国土数値情報の郵便局データを利用する

国土交通省で公開している国土数値情報の郵便局データを利用して各都道府県の郵便局の所在地をプロットする方法です。なお、データは全国（ファイル名：P30-13.zip）を使用しています。

国土数値情報の郵便局データ：https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P30.html#prefecture00

#必要パッケージの読み込みとインストール
if(!require("rgdal", quietly = TRUE)){
  install.packages("rgdal");require("rgdal")
}
if(!require("tidyverse", quietly = TRUE)){
  install.packages("tidyverse");require("tidyverse")
}
library("tcltk")
library("tmap")

#ダウンロードしたshpデータを読み込み
#読み込みに時間がかかります
JPNPref <- readOGR(paste0(as.character(tkgetOpenFile(title = "shpファイルを選択",
                                                     filetypes = '{"shpファイル" {".shp"}}',
                                                     initialfile = c("*.shp")))),
                   GDAL1_integer64_policy = TRUE, use_iconv = TRUE,
                   #環境によっては"utf-8"に変更
                   encoding = "shift-jis")

#列名を認識できるように変更
JPNPref@data <- JPNPref@data %>%
  rename(latitude = X, longtude = Y)

#「ggplot2」パッケージでプロット
tmap_mode("plot")
tm_shape(JPNPref) +
  tm_bubbles(size = .0001,
             col = "blue",
             border.col = NA,
             border.lwd = NA,
             alpha = 0.2) +
  tm_style("gray")

#「leaflet」パッケージでプロット
tmap_mode("view")
tm_shape(JPNPref) +
  tm_bubbles(size = .0001,
             col = "blue",
             border.col = NA,
             border.lwd = NA,
             alpha = 0.2) +
  tm_style("gray")

出力例

郵便局の所在位置で日本列島が浮かび上がりました。

少しでも、あなたの解析が楽になりますように！！

タイトルとURLをコピーしました