市区町村役所間の距離行列を求める~pdftoolsでのデータ抽出とsfによる算出編~
はじめに
Qiitaのなかで@3tkyさんが、都道府県庁間の距離行列をRで求める記事を書かれていました。
そして、挑戦状のような一文が残されています。
市区町村の役所の緯度経度は 国土地理院 がまとめているので、PDFを整形し、60進法表記を @uri さんのkuniezuパッケージ や pazerパッケージを活用して10進法に変換後、上のコードから市区町村役所間の距離行列が作成できそう。
これはやらねばなるまい（謎の使命感）。 本記事では、[@3tky](https://qiita.com/3tky) さんの書かれていた都道府県庁間の距離行列を算出する方法の別解と同様に市区町村役所の距離行列を求める方法を書きます。 大元の市区町村の役所の緯度経度データは国土地理院が公開するPDFデータです。そのためPDFからデータを抽出、加工する方法についても言及します。
記事の中で以下のパッケージを利用します。予め読み込んでおきましょう。
library(sf) library(ggplot2) library(rvest) # library(units) # library(tidyverse)
都道府県庁間の距離行列
記事を見たあと、CRANにリリース準備だった kuniezu パッケージに都道府県の県庁位置のデータを収納しました。 このデータの取得・整形するコードは
https://github.com/uribo/kuniezu/blob/master/data-raw/office_locaiton.R
に置いています。
data("jp47prefectural_offices", package = "kuniezu") # 47都道府県の県庁位置データ jp47prefectural_offices
## Simple feature collection with 47 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 127.681114197 ymin: 26.2124996185 xmax: 141.346939087 ymax: 43.0641670227
## CRS: EPSG:4326
## # A tibble: 47 x 2
## office geometry
## <chr> <POINT [°]>
## 1 北海道庁 (141.346939087 43.0641670227)
## 2 青森県庁 (140.740005493 40.8244438171)
## 3 岩手県庁 (141.152496338 39.7036094666)
## 4 宮城県庁 (140.871948242 38.2688903809)
## 5 秋田県庁 (140.102493286 39.7186126709)
## 6 山形県庁 (140.363327026 38.2405548096)
## 7 福島県庁 (140.467773438 37.75)
## 8 茨城県庁 (140.446670532 36.3413887024)
## 9 栃木県庁 (139.883605957 36.5658340454)
## 10 群馬県庁 (139.060836792 36.3911094666)
## # … with 37 more rows
まずは対角要素が0となる距離行列を作成します。sfオブジェクトでの地物間の距離は st_distance() を使って求められます。
dist_pref <- jp47prefectural_offices %>% st_distance()
colnames(dist_pref) <- jp47prefectural_offices$office rownames(dist_pref) <- jp47prefectural_offices$office dist_pref[seq_len(6), seq_len(6)]
## Units: [m]
## 北海道庁 青森県庁 岩手県庁 宮城県庁
## 北海道庁 0.000000 253808.763222 373582.0205556 534013.529177
## 青森県庁 253808.763222 0.000000 129308.2699883 283959.259182
## 岩手県庁 373582.020556 129308.269988 0.0000000 161119.538069
## 宮城県庁 534013.529177 283959.259182 161119.5380689 0.000000
## 秋田県庁 385850.690256 134229.605213 90055.4114243 174198.478652
## 山形県庁 542058.359064 288699.213957 176229.8003421 44629.402174
## 秋田県庁 山形県庁
## 北海道庁 385850.6902557 542058.359064
## 青森県庁 134229.6052134 288699.213957
## 岩手県庁 90055.4114243 176229.800342
## 宮城県庁 174198.4786523 44629.402174
## 秋田県庁 0.0000000 165635.769370
## 山形県庁 165635.7693696 0.000000
また、距離の単位はmですが units パッケージの関数を使って任意の単位に変更可能です。
dist_pref <- dist_pref %>% units::set_units(km) dist_pref[seq_len(6), seq_len(6)]
## Units: [km]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000000 253.808763222 373.5820205556 534.013529177 385.8506902557
## [2,] 253.808763222 0.000000000 129.3082699883 283.959259182 134.2296052134
## [3,] 373.582020556 129.308269988 0.0000000000 161.119538069 90.0554114243
## [4,] 534.013529177 283.959259182 161.1195380689 0.000000000 174.1984786523
## [5,] 385.850690256 134.229605213 90.0554114243 174.198478652 0.0000000000
## [6,] 542.058359064 288.699213957 176.2298003421 44.629402174 165.6357693696
## [,6]
## [1,] 542.058359064
## [2,] 288.699213957
## [3,] 176.229800342
## [4,] 44.629402174
## [5,] 165.635769370
## [6,] 0.000000000
この値が @3tky さんの結果や国土地理院の公表データと一致していることを確認します。
市区町村役所の距離行列
続いて市区町村役所の距離行列を求めます。 こちらのデータも国土地理院が位置情報を整理しています。 ですがPDFなので、Rで扱う際にはデータを抽出する必要が生じます。
PDFファイルのダウンロードは以下のコードで行います。
df_link <- read_html("https://www.gsi.go.jp/KOKUJYOHO/center.htm") %>% html_nodes(css = "div.base_txt > div:nth-child(5) > table > tbody > tr > td > a") %>% { tibble::tibble( name = stringr::str_remove(html_text(.), "\\[.+\\]"), link = html_attr(., "href")) } df_link$link %>% purrr::walk( ~ download.file(url = .x, destfile = basename(.x))
RでPDFのデータ抽出を行うパッケージはいくつかありますが、今回は
pdftools を使いました。
pdftools::pdf_text()によりテキストを抽出、若干血生臭い文字列処理を行いデータフレーム化します。
これらの処理を一括で実行する関数を書きました。
リンク先のコード（関数部分）をコピーペーストで読み込んでください。 最終的には各都道府県のPDFを引数に与える
gsi_office_extract() を実行することでデータ抽出が完了します。
# 茨城県の市区町村役所の位置情報 sf_pref08office <- gsi_office_extract("ibaraki_heso.pdf")
sf_pref08office
## Simple feature collection with 45 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 139.745285034 ymin: 35.8577766418 xmax: 140.751113892 ymax: 36.8019447327
## CRS: EPSG:4326
## # A tibble: 45 x 2
## office geometry
## <chr> <POINT [°]>
## 1 茨城県庁 (140.446670532 36.3413887024)
## 2 水戸市役所 (140.471664429 36.3658332825)
## 3 日立市役所 (140.650558472 36.5988883972)
## 4 土浦市役所 (140.204162598 36.0783348083)
## 5 古河市役所 (139.755004883 36.1783332825)
## 6 石岡市役所 (140.286941528 36.1905555725)
## 7 結城市役所 (139.876663208 36.305557251)
## 8 龍ケ崎市役所 (140.182220459 35.9116668701)
## 9 下妻市役所 (139.967498779 36.1844444275)
## 10 常総市役所 (139.993896484 36.0236129761)
## # … with 35 more rows
この位置情報データに対して、都道府県庁間の距離行列を求めた時と同じコードを実行することで市区町村役所の距離行列も求められます。 これらの処理も関数化しておきました。可視化と合わせてどうぞ。
# 距離行列を作成 st_distmatrix <- function(data, var) { res <- data %>% st_distance() %>% units::set_units(km) vars <- data %>% purrr::pluck(var) colnames(res) <- vars rownames(res) <- vars res } # ggplot2で描画するためのデータ整形。三角行列に変換して縦長のデータにします matrix_to_longer <- function(data) { res_mt <- as.matrix(units::drop_units(data)) res_mt[lower.tri(res_mt)] <- NA res_mt %>% as.data.frame() %>% tibble::rownames_to_column(var = "from") %>% tibble::as_tibble() %>% tidyr::pivot_longer(cols = seq.int(2, ncol(.)), names_to = "to", values_to = "dist") %>% dplyr::filter(!is.na(dist)) %>% dplyr::mutate_at(dplyr::vars(from), list(~ forcats::fct_inorder(.))) %>% dplyr::mutate(to = forcats::fct_rev(forcats::fct_inorder(to))) } # ggplot2でのプロット plot_distmatrix <- function(data) { data %>% ggplot(aes(from, to)) + geom_tile(aes(fill = dist), color = "white") + scale_fill_viridis_c() + guides(fill = guide_colorbar(title = "distance (km)")) + theme_bw(base_family = "IPAexGothic") + theme( axis.text.x = element_text(angle = 90), axis.ticks = element_blank(), axis.line = element_blank(), panel.border = element_blank(), panel.grid.major = element_line(color = '#eeeeee') ) }
st_distmatrix(sf_pref08office, "office") %>% matrix_to_longer() %>% plot_distmatrix()
Enjoy!
kuniezu: 日本の国土地理を扱いやすくするRパッケージをCRANに登録しました
はじめに
kuniezuパッケージ (v0.1.0) をCRANにリリースしました。
このパッケージは、私が業務や趣味で日本国内の地理空間データを扱う時に作っていた関数を一つのパッケージに整理したものです。 空間的に世界規模のデータを扱うのではなく、日本国内に限った話であれば、日本に即した仕様や座標参照系を利用した方が良いことがあります。 そうした日本の地理空間データを処理する際に利用することがある機能や、あると便利なデータセットを提供できるように努めています。
ゆるゆると開発するつもりでいましたが、Twitterで開発中であることをつぶやいたところ想定以上に反響が大きかったのでCRAN登録を急ぐことにしました。
仕事でも使うし、誰かの役に立つかと思い、日本の地理空間データのための処理を関数化し、パッケージとして開発始めました。🗾ブログに書いた、度分秒での経緯度の変換や平面直角座標系の特定とか。機能は少ないですが使ってみてください。（フィードバック、機能要望もぜひ） https://t.co/VUaPyC6MaT pic.twitter.com/bwaFL0vJqw
— Uryu Shinya (@u_ribo) 2020年4月27日
この記事では、そんなkuniezuパッケージの機能紹介を行います。
まずはみなさんパッケージがインストールされていないと思うので、利用されたい方は以下のコマンドを実行してインストールを済ませてください。
install.packages("kuniezu")
library(kuniezu) # 使い方を説明するために以下のパッケージも読み込んでおきます library(sf) library(ggplot2) library(leaflet)
使い方
最初のリリースでは以下の機能を実装しました。
- DMS表記を十進数表記に変換
- 日本測地系2011における平面直角座標系の特定
- 南西諸島・小笠原諸島を移動した日本地図の描画
- 地理院タイルをleafletで簡単に利用できるように
- 国土地理に関するいくつかのデータセット
順に解説します。
parse_*_dohunbyo(): DMS表記を十進数表記に変換
地球上の任意の座標を示す際に使われる緯度と経度の表記には、十進数で表される場合と「度(Degree)・分(Minute)・秒
(Second)」を用いて表記される場合があります。
後者はDMS表記と呼ばれ、十進数で139.7414、35.6581と表される座標に対してE139°44′28.8869″、N35°39′29.1572″のように表現します。
更に日本では度分秒の記号に対して漢字使って東経139度44分28.8869秒、北緯35度39分29.1572秒等の表記が用いられることがしばしばあります。
この日本独自のDMS表記を十進数に変換する関数がparse_*_dohunbyo()です（アスタリスクの部分にはlon、latがそれぞれ指定できます）。
lonは東経、西経を示す(180~-180）経度、latが北緯、南緯（90~-90）の緯度を扱うために使い分けます。
# 日本経緯度原点（東京都港区麻布台二丁目）の座標 parse_lon_dohunbyo("東経139度44分28.8869秒")
## [1] 139.7414
parse_lat_dohunbyo("北緯35度39分29.1572秒")
## [1] 35.6581
これは過去のブログ記事に書いた関数をパッケージに含めたものです。
加えて、十進数に変換する前段階として漢字を記号に処理するreplace_dohunbyo_kanji()も用意しています。
replace_dohunbyo_kanji("東経139度44分28.8869秒")
## [1] "E139°44’28.8869."
日本測地系2011における平面直角座標系の特定
日本の現行の測地系である日本測地系2011には、19の平面直角座標系の区分がなされています。
平面直角座標系は狭い範囲の測量や距離の算出に用いるのに適しています。
任意の経緯度を平面直角座標に変換するには、まず座標が平面直角座標系のどの系番号に含まれるのかを知っておく必要が生じます。
これを求めるのがst_detect_jgd2011()です。
sfパッケージでおなじみのst_*()が関数名につくように、関数にはsfcオブジェクトを与えて実行します。
# 日本経緯度原点の座標 st_detect_jgd2011(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326))
## [1] 6677
move_jpn_rs(): 南西諸島の一部・小笠原諸島を移動した日本地図の描画
日本だけを地図上に描画する際、表示する画面の都合からか、南西諸島の一部や小笠原諸島の位置を海上に移動することがあります（例えばNHKの気象情報のページや白地図）。 これを実現するには、予め移動対象の地物に対する操作が必要となります。
move_jpn_rs()はこれらの島々の位置を変更するための関数です。
地図描画のために奄美群島と沖縄県、小笠原諸島の位置を変更する処理が実行されます。
また、この地図化もパッケージの関数で行えるようにggplot2ベースのgeom_jpsegment()を用意しました。
使い方は次のように、まず都道府県や市区町村レベルでポリゴンが分かれた地理空間データに対してmove_jpn_rs()を適用、その後マッピングを行います。
move_jpn_rs(jgd2011_bbox) %>% ggplot() + geom_sf() + geom_jpsegment()
デフォルトの設定で移動した地物がわかるように境界線が引かれます。
この関数の元も、かつて記事に書いたものです。
少し前に @shoei05 さんがブラッシュアップしてくれたのをきっかけにパッケージへ含めました。
@u_ribo いつも勉強させていただいています。以前の記事（https://t.co/iJxfDxFhJ3 ）に関連して、コロプレスマップなど作図する場合は、sf_ja_omit47 とするのではなく、やや強引ですが同一セット内の $geometry に座標を足してあげると良さそうだったので、少しスクリプトをいじってみました。 pic.twitter.com/EI9NvHMpnD
— Shoei/将英 (@shoei05) March 26, 2020
一方で、まだこの関数はさらなる磨き上げが必要で、解決策がわかる人はぜひ教えてください。
地理院タイルをleafletで簡単に利用できるように
インタラクティブなウェブ地図を利用したことがある方ならわかるかと思いますが、leafletでは任意の地図タイルを背景画像として扱うことができます。 国土地理院が発行する地理院タイルもまた、この背景画像として利用可能で、 一覧ページに示されるように多種多様な地図タイルが整備されています。
kuniezuパッケージではこの地理院タイルをleafletで簡単に導入できるようになっています。 これらはgsi_tilesオブジェクトを通して呼び出します。
現在のバージョンでは「近年の災害」をはじめとして、すべての地理院タイルに対応しているわけではありません。 kuniezuで利用可能な地理院タイルの一覧は以下のコードで確認できます。
names(gsi_tiles)
## [1] "standard" "pale" "english"
## [4] "lcm25k_2012" "lcm25k" "ccm1"
## [7] "ccm2" "vbm" "vbmd_bm"
## [10] "vbmd_colorrel" "vbmd_pm" "vlcd"
## [13] "lum200k" "lake1" "lakedata"
## [16] "blank" "seamlessphoto" "ortho"
## [19] "airphoto" "gazo4" "gazo3"
## [22] "gazo2" "gazo1" "ort_old10"
## [25] "ort_USA10" "ort_riku10" "relief"
## [28] "anaglyphmap_color" "anaglyphmap_gray" "hillshademap"
## [31] "earthhillshade" "slopemap" "slopezone1map"
## [34] "afm" "lcmfc2" "lcmfc1"
## [37] "swale" "cp" "jikizu2015_chijiki_d"
## [40] "jikizu2015_chijiki_i" "jikizu2015_chijiki_f" "jikizu2015_chijiki_h"
## [43] "jikizu2015_chijiki_z" "jikizu_chijikid" "jikizu_chijikii"
## [46] "jikizu_chijikif" "jikizu_chijikih" "jikizu_chijikiz"
使い方は次のようになります。leafletおよびmapviewでの例を紹介します。
# 出力結果は省略します gsi_tiles$standard %>% addCircles( data = sf::st_transform(extreme_points %>% purrr::reduce(c), crs = 4326)) mapview::mapview(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326), map = gsi_tiles$pale)
なお、これらの地理院タイル利用する際はタイルごとに設定されている利用規約に従ってください。
国土地理に関するデータセット
これまでの例で出てきたものがありますが、日本の国土地理を扱う上であると便利なデータ等をパッケージで提供するようにしました。 現在、次のようなデータがあります。
- extreme_points… 東西南北端点の経度緯度データ
- jgd2011_bbox… 日本測地系2011の平面直角座標系の区域を示すポリゴンデータ
- jp47prefectural_offices… 47都道府県の庁舎の位置データ
こんなデータや機能が欲しいという方はぜひ issue に書き込んでください。 バグ報告もありがたいです。 GitHub以外ではTwitterでコメントいただいても対応します。
Enjoy!
ある座標からの指定半径に含まれるメッシュコードを知る
新型コロナウイルスのデータを扱う際に、メッシュコード（標準地域メッシュ）が利用されることがあります。 特にNTTドコモ「モバイル空間統計」分析レポートのデータは、内閣官房の新型コロナウイルス感染症対策のページにも掲載されているように 全国各地の人口変動を分析するのに欠かせないデータとなっています。
人流データを使った分析は、先日発表された「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日） (PDF)の中でも 行われており、本文中に以下の記述があります。
渋谷駅周辺と難波駅周辺から半径 1 ㎞圏内においては、10 歳台および 20 歳台の若者を中心として昼夜問わず接触頻度が 80％以上、減少したことがうかがえる。
また、GitHubにアップロードされている
この結果の補足資料を見るとモバイル空間統計のメッシュデータを使った分析と可視化の事例が確認できます（3. 各エリアの接触頻度と変化率）。
補足資料の図2など、対象の駅周辺でのメッシュごとの接触の変化率を表示する際、駅周辺の半径1km圏内とそれに含まれる500mメッシュの枠が示されています。
緊急事態宣言中の人流データを整理した資料をクラスター対策班がまとめたので共有します。
— 新型コロナクラスター対策専門家 (@ClusterJapan) 2020年5月1日
難しい内容になっていますので、読むのに時間がかかるかもしれませんが、質問いただければ私たちとしてできるだけ回答します。がんばります。
▼緊急事態宣言中の人流データーまとめhttps://t.co/NKiVE8oFYC
この図を見た直後、これはRでできるぞー！と思いました（もともとRでやられている、より効率的な処理を実行されているのかもしれませんが...）。 ブログの主題にあるように「ある座標からの指定半径に含まれるメッシュコードを知る」ことができれば簡単に実行できます。
真似事になりますが、メッシュコードを使った分析例として優れており、応用範囲の広いトピックスですのでRでやる方法を整理しておきます。
Rコード
library(dplyr) library(mapview) library(sf) library(jpmesh)
まず対象の座標からPOINTのsfcオブジェクトを作成します。
sfパッケージではst_point()に任意の座標をベクトルで与えることでPOINTの作成が行われます。
これを測地系WGS84の座標参照系、EPSG:4326として扱えるようにst_sfc()で指定します。
st_point()で定義する座標は、500mメッシュコードに変換した時に 533935961 となる渋谷駅の経緯度です。
p_shibuya_st <- # 目測で座標を決めているので厚生労働省のポイントとはズレがあります st_point(c(139.70172, 35.65846)) %>% st_sfc(crs = 4326)
続いてこの座標点を起点としたバッファ領域（緩衝帯）を生成します。これにはst_buffer()を使いますが、EPSG:4326の地球楕円体モデルは先述の通りWGS84です。
これは水平位置を表す経緯度と垂直位置を表す高度との組み合わせからなる、3次元の地理座標系です。
2地点間の距離やバッファを求める際には地物を平面に投影し、XY座標で表す投影座標系である平面直角座標系を利用すると良いので st_transform()により座標参照系の変更を行います。
ここでは日本測地系2011の平面座標系の一つであるEPSG:6677（福島県、栃木県、茨城県、埼玉県、千葉県、群馬県、神奈川県、島嶼部を除いた東京都が含まれる）を指定しました。
b1km_shibuya_st <- p_shibuya_st %>% st_transform(crs = 6677) %>% st_buffer(dist = units::set_units(1, km)) %>% st_transform(crs = 4326)
上記の例ではバッファの範囲を半径1kmとしました。処理のあと、st_transform()により再び元の座標参照系に戻しています。
一度確認してみましょう。白地図では位置関係が掴みにくいので地図タイルを背景に、データを重ます。
m <- mapview(p_shibuya_st) + mapview(b1km_shibuya_st) m
座標とバッファ領域が表示されました。
続いて、この範囲に含まれる500mメッシュを特定します。 ここではメッシュコードが未知のものとして座標を与えて探し出すところから行う例を示します。 座標の位置する80kmメッシュから、それに含まれる10kmメッシュコードの生成、ポリゴン化を以下のコードで行います。 なお1kmメッシュから探索を開始しないのは、候補となるメッシュの数を減らし、実行時間を短縮するためです。
mesh_candidate <- p_shibuya_st %>% coords_to_mesh(geometry = ., mesh_size = 80) %>% mesh_convert(to_mesh_size = 10) %>% export_meshes() nrow(mesh_candidate)
10kmメッシュのポリゴンが用意できたら、バッファのポリゴンに対する空間関係を調べます。
st_join()は2つの地物の空間関係をもとにデータを結合するために使われます。
今回は2つの地物（10kmメッシュ、バッファポリゴン）に共有部分がある場合にデータを結合するように join = st_intersectsを引数で与えます（st_join()の既定値）。
さらに結合後のデータを共有部分だけのものにしたいのでleft = FALSEを指定しました。
mesh_candidate <- mesh_candidate %>% st_join(st_sf(b1km_shibuya_st), join = st_intersects, left = FALSE) mesh_candidate
実行結果を見ると、10kmメッシュのデータが2件に絞られています。 これはバッファポリゴンと交差するポリゴンを持っている10kmメッシュコードが2つあることを示します。
残りの処理は対象のメッシュに対してスケールダウン、バッファとの空間関係を調べて抽出するを繰り返すだけです。 ただし、10kmから1km、500mのスケールダウンを何度も行うのは手間ですので、ここで10kmメッシュに含まれる500mメッシュを用意して処理を簡略化します。 1つの10kmメッシュコードに含まれる500mメッシュの数は400、そのため800個の500mメッシュから最終的なメッシュコードを抽出することになります。
mesh_target <- mesh_candidate %>% pull(meshcode) %>% purrr::map( ~ mesh_convert(.x, to_mesh_size = 0.5) ) %>% purrr::reduce(c) %>% unique() %>% export_meshes() %>% st_join(st_sf(b1km_shibuya_st), join = st_intersects, left = FALSE) nrow(mesh_target)
最終的な渋谷駅周辺の半径1kmに含まれる500mメッシュの数は22個となりました。 確認のためにもう一度データを表示してみます。
m + mapview(mesh_target)
バッファ領域とそれに含まれるメッシュの抽出ができました。 いくらか手間がかかりますがRでもできます、という話でした。
Enjoy!
おまけ: 出力のための調整
このマップをさらに報告書などのPDFに添付することがある場合を想定し、その際の調整方法を書いておきます。 具体的にはマップに対して
- 凡例の変更... 元のマップには凡例もないので載せる
- ファイルへの出力... 厚生労働省の補足資料では切り抜きが雑なので統一する
を行います。これらもRで実行します。
まずはマッピングするデータを用意します。手元にデータがないので図2 左の値をデータ化します。
df_value <- tibble::tibble( meshcode = as.character(c(533945052, 533945061, 533935953, 533935954, 533935963, 533935964, 533935951, 533935952, 533935961, 533935962, 533935971, 533935853, 533935854, 533935863, 533935864, 533935873, 533935851, 533935852, 533935861, 533935862, 533935754, 533935763)), value = c(-51, -60, -22, -49, -53, -54, -35, -60, -62, -57, -33, -35, -55, -60, -52, -31, -15, -35, -38, -38, -27, -35))
続いて凡例で使われるカラーパレットを定義しておきます。変化率は正負いずれの値も取り得るとし、 原点(0)を中心に赤系（増加）と青系（減少）を示す10刻みのカラーパレットを作成します。
class_int <- classInt::classIntervals(seq.int(-100, 100, by = 10), n = 20, style = "fixed", fixedBreaks = seq.int(-100, 100, by = 10)) # カラーパレットの確認 pals::pal.bands(pals::ocean.balance(20))
df_map <- mesh_target %>% left_join(df_value, by = "meshcode") %>% mutate(value_class = cut(value, class_int$brks, include.lowest = TRUE)) res_map <- df_map %>% mapview(zcol = "value_class", layer.name = "変化率", col.regions = pals::ocean.balance(20), homebutton = FALSE) %>% leafem::addStaticLabels( data = df_map, textsize = "18px", label = paste0(df_map$value, "%"))
mapshot(res_map, file = "out.png", remove_controls = c("zoomControl", "layersControl", "homeButton"))
以上です。
