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

#パッケージの読み込み
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")

#パッケージの読み込み
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")

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

#パッケージの読み込み
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")

#必要パッケージの読み込みとインストール
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")