本文へスキップ
からだにいいもの

Rのトピックスを中心に『まだ、まだ、知らない、役に立つ情報？』を発信します。

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

Rの解析に役に立つ記事

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

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

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

＜おすすめのRに関する書籍です＞

言葉と数式で理解する多変量解析入門 | 小杉 考司
言葉と数式で理解する多変量解析入門 | 小杉 考司
Amazonで小杉 考司の言葉と数式で理解する多変量解析入門。アマゾンならポイント還元本が多数。小杉 考司作品ほか、お急ぎ便対象商品は当日お届けも可能。
www.amazon.co.jp

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

出力例

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

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

＜おすすめのRに関する書籍です＞

サラっとできる!フリー統計ソフトEZR(Easy R)でカンタン統計解析 | 神田 善伸
サラっとできる!フリー統計ソフトEZR(Easy R)でカンタン統計解析 | 神田 善伸
Amazonで神田 善伸のサラっとできる!フリー統計ソフトEZR(Easy R)でカンタン統計解析
www.amazon.co.jp

SHARE
ポストはてブPocket
スポンサーリンク

関連記事

Rで解析：HTMLコードを依存ゼロで作成「hypertext」パッケージの紹介 Rで解析：ggplot2にGraphPad Prism風のテーマを適応「ggprism」パッケージ Rで解析：ggplot2の図の体裁操作がちょっとは楽になるかも「ggeasy」パッケージ Rで解析：MATLABのコマンドが使えます「matlab」パッケージ
Prices and shipping availability may change. Please refer to the product page at time of purchase.
Content displayed on this site is provided by Amazon and may be updated or removed.
Amazon Associate, karada-good earns income through qualifying sales.