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

ggplot2パッケージと同様レイヤーの概念を持ち、数量データを示す地図を作成するのに便利なパッケージの紹介です。データをもとに分割して地図を作成することが可能です。

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

パッケージバージョンは1.10。windows 10のR version 3.3.3で動作を確認しています。

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

先日紹介した「Rで解析：「leaflet」パッケージで使える日本のシェープファイルの作成例」で作成したシェープファイルに気象庁より観測所一覧を取得してシェープファイルに付与します。

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

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

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

###気象庁より観測所一覧を取得#####
#Tmpフォルダに保存し解凍後に読み込む
temp <- tempfile()
download.file("http://www.jma.go.jp/jma/kishou/know/amedas/ame_master.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))
LonData <- type.convert(paste0(MasterKansokujyo[, 5] + MasterKansokujyo[, 6]/60))
#シェープファイルに付与
library("stringi")
attributes(JPNPref)$Kansokujyo <- data.frame("都道府県" = stri_encode(MasterKansokujyo[, 1], "", "utf-8"),
                                             "観測所名" = stri_encode(MasterKansokujyo[, 2], "", "utf-8"),
                                             "経度" = LatData, "緯度" = LonData)
#シェープファイルを保存
save(JPNPref, file = "JPNPref.shp")

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

作成した日本地図シェープファイルを読み込み後、以下のコードを実行します。例では日本地図シェープファイルのデータ名は「JPNPref」です。

###気象庁より最高気温最新を取得#####
#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")

#各都道府県の最高気温を抽出
#データ保管用変数
NewMaxTemp <- data.frame()
#処理
for(n in 1:47){
#都道府県を抽出
GetPrefData <- MaxTemp[which(MaxTemp[, 2] %in% grep(JPNPref@data$gns_name[n], MaxTemp[, 2], value = TRUE)),]
#最高気温を降順で並び替え
GetPrefData <- GetPrefData[order(GetPrefData[, 10], decreasing = TRUE),]
#観測所の対象都道府県
GetKansokujyo <- JPNPref@Kansokujyo[which(JPNPref@Kansokujyo[, 1] %in%
                                            grep(substr(JPNPref@data$gns_name[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$gns_name[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])
#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,NewMaxTemp以外を削除
rm(GetPrefData, MaxTemp, n)

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

下記コマンドを実行してください。

#パッケージのインストール
install.packages("tmap")

コマンドの紹介

詳細はコマンド、パッケージのヘルプを確認してください。

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

#プロット出力を選択する: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 = .01, "red") +
  tm_style_gray()

#例えば気温で分割してプロット
#気温を区分分け:「fancycutp」パッケージを利用
#https://www.karada-good.net/analyticsr/r-544
#install.packages("fancycut")
library("fancycut")
#分類
CutLabel <- fancycut(x = JPNPref@data[, 5],
                     intervals = c("[10, 20)", "[20, 25)", "[25, 30)", "[30, 35]"),
                     buckets = c("10-20", "20-25", "25-30", "30-35"),
                     unmatched.bucket = "範囲外")
#結合
JPNPref@data <- cbind(JPNPref@data, CutLabel)
#プロット
tm_shape(JPNPref) +
  tm_polygons("気温", id = "都道府県:地点名", legend.show = TRUE) +
  tm_facets(by = "CutLabel") +
  tm_shape(MaxSpDF) +
  tm_bubbles(size = .01, "red") +
  tm_style_gray()

出力例

2017/05/13 13：00頃に取得したデータです。
・簡易プロット:qtmコマンド

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

・例えば気温で分割してインタラクティブにプロット
いろいろいじってみてください。
　・別ウィンドウで開く
　https://www.karada-good.net/wp/wp-content/uploads/2017/05/tmapfacet.html

少しでも、あなたの解析が楽になりますように！！けいおんは素晴らしい。