以下の記事がとても興味深かったので、記事を参考にふるさと納税の納税額や使用用途をマップ上に表示するShinyアプリを作成しました。
以下のコードを実行することでShinyをローカル環境で実行することができます。
if(!require(shiny)){install.packages("shiny")};runGitHub("furusato_map", "sb8001at")
地図情報の取得
日本の行政区分の地理情報として、国交省情報活動推進課が運営している国土数値情報ダウンロードサイトからgeojsonファイルをダウンロードします。
出典：国土交通省国土数値情報ダウンロードサイト 行政区域データ(2024年) 全国データ
geojsonファイルをRで利用する場合には、sf::st_readで読み込みます。st_readでgeojsonを読み込むことでgeojsonのデータをsfクラスのオブジェクトとして利用することができます。
pacman::p_load(tidyverse, sf, leaflet, readxl)
geo <- st_read("N03-20240101.geojson")
geoの中身は以下のようなファイルになっています。sfはデータフレームと同様に取り扱うことができるオブジェクトで、地理情報はgeometryの列に記録されています。このgeometryの列はリストのリスト（sfc）になっていて、個別のデータはさらにリスト（sfg）となりますが、地図を描画するだけならデータ型は特に気にすることはないでしょう。
> geo
Simple feature collection with 124134 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 122.9326 ymin: 20.42275 xmax: 153.9867 ymax: 45.55724
Geodetic CRS: WGS 84
First 10 features:
N03_001 N03_002 N03_003 N03_004 N03_005 N03_007 geometry
1 北海道 石狩振興局 <NA> 札幌市 中央区 01101 POLYGON ((141.2569 42.99782...
2 北海道 石狩振興局 <NA> 札幌市 北区 01102 POLYGON ((141.3333 43.07505...
3 北海道 石狩振興局 <NA> 札幌市 東区 01103 POLYGON ((141.3734 43.0684,...
4 北海道 石狩振興局 <NA> 札幌市 白石区 01104 POLYGON ((141.382 43.04832,...
5 北海道 石狩振興局 <NA> 札幌市 豊平区 01105 POLYGON ((141.3637 42.94154...
6 北海道 石狩振興局 <NA> 札幌市 南区 01106 POLYGON ((141.2354 42.82273...
7 北海道 石狩振興局 <NA> 札幌市 西区 01107 POLYGON ((141.169 43.0829, ...
8 北海道 石狩振興局 <NA> 札幌市 厚別区 01108 POLYGON ((141.5048 43.0211,...
9 北海道 石狩振興局 <NA> 札幌市 手稲区 01109 POLYGON ((141.227 43.08302,...
10 北海道 石狩振興局 <NA> 札幌市 清田区 01110 POLYGON ((141.4138 42.89511...
geoは124134行、元のjsonファイルは463MB、オブジェクトサイズは334MBと巨大で、取り扱いが難しいので、要約しファイルサイズを調整します。
都道府県＋市町村の列（city）を追加し、aggregate.sf関数（aggregate関数のジェネリック）を用いて市町村ごとにデータをまとめます。データをまとめるための列としてnullcolという列を一時的に追加しています（N03_007の列を使っても良いかもしれません）。
geo$city <- paste0(geo$N03_001, " ", geo$N03_004, geo$NO3_005)
geo$nullcol <- 0
geo2 <- geo[,"nullcol"] |> aggregate(by=list(geo$city), sum)
geo2$nullcol <- NULL
colnames(geo2) <- c("city", "geometry")
これでgeo2は1753行のデータとなります。かなり圧縮されたことがわかります。しかし、オブジェクトサイズはgeoが334MB、geo2が280MBとまだ巨大です。これは、地図情報の解像度が高すぎる（POLYGONの点が多すぎる）ためです。コロプレス図を作るのにこれほど高解像度のデータは必要ないため、sf_simplify関数を用いて解像度の低いPOLYGONに変換します。preserveTopology=TRUEで形状をあまり変えないようにし、dToleranceで解像度の下げ具合を指定します。dToleranceはメートル単位とされていますので、以下の例では300m単位まで地理情報を粗くしていると思われます。
geo2 <- st_simplify(geo2, preserveTopology = TRUE, dTolerance = 300)
これでgeo2は3.5MBのオブジェクトとなり、取り扱いが簡単になります。ただし、後で述べるように4行のデータのgeometryが失われる（with 4 geometries empty）ので、後ほどこの行を取り除く必要があります。
国勢調査データとマージする
次に、人口統計とデータを比較するために、以下の政府統計の総合窓口（e-Stat）から国勢調査のデータを拾ってきます。市町村別の人口統計は（私が調べた限りでは）国勢調査時のものしか無いようですので、人口は5年おき（2015年、2020年）のものしか（おそらく）利用できません。
出典：調査項目を調べる－国勢調査（総務省）令和２年国勢調査 人口等基本集計 （主な内容：男女・年齢・配偶関係，世帯の構成，住居の状態，母子・父子世帯，国籍など）
ダウンロードしたデータを整形します。データを読み込むと数値が文字列になりがちですので、一通り数値に置き換えます。NAが生じますがあまり気にしないことにします。
d <- read_excel("b01_01.xlsx", skip=14)
d <- d |> select(c(5, 7:13, 16, 17))
colnames(d) <-
c("pref", "townname", "pop", "pop_M", "pop_F",
"pop_2015", "pop_change","pop_change_rate", "pop_density", "household")
d$pref <- str_sub(d$pref, start=4, end=1000)
d$townname <- str_sub(d$townname, start=6, end=1000)
d$pop <- as.numeric(d$pop)
d$pop_M <- as.numeric(d$pop_M)
d$pop_F <- as.numeric(d$pop_F)
d$pop_2015 <- as.numeric(d$pop_2015)
d$pop_change_rate <- as.numeric(d$pop_change_rate)
d$household <- as.numeric(d$household)
d <- d |> filter(!str_detect(townname, "（旧："))
d$city <- paste(d$pref, d$townname)
整形時に地理空間情報と同じく都道府県＋市町村の列（city）を追加しているので、この列を参照し、地理空間情報と国勢調査のデータをマージします。マージにはdplyr::left_joinを用いると良いでしょう（左側がsfでないと、地理情報を失った行が生じます）。マージしたデータには色々用途がありそうなので、とりあえずst_write関数でgeojsonとして保存しておきます。
geo2 <- left_join(geo2, d, by=join_by(city))
st_write(geo2, "census2020_geodata.geojson")
ふるさと納税のデータをマージする
上で準備したsfオブジェクトにふるさと納税のデータをさらにマージしていきます。
出典：総務省 ふるさと納税ポータルサイト 関連資料 受入額の実績等
データのヘッダーと取り込む列は年度により異なるので、read_excel関数のskipの数値・列選択（select）を調整して、同じ形のデータを取り込めるようにします。地理空間情報とのマージには上と同じくleft_join関数を用います。また、市町村外から受け入れた寄付金額、寄付件数を人口当たりとして計算します（ただし、人口は2020年度のものなので、2020年度の結果以外は正確ではありません）。
geo <- st_read("census2020_geodata.geojson")
d <- read_excel("results20170704-02.xlsx", skip=5)
d <- d |> select(c(2:5, 7, 8, 25:31))
colnames(d) <- c(
"pref", "townname", "Number_donations_accepted",
"Donation_amount_accepted", "Number_donations_from_outside", "Donation_amount_from_outside",
"Costs_procure_gifts","Costs_return_gifts", "Costs_PR", "Costs_settlement", "Costs_Administration", "Costs_others", "Total_costs")
d$city <- paste(d$pref, d$townname)
d$pref <- NULL
d$townname <- NULL
geo2 <- left_join(geo, d, by=join_by(city))
geo2$Donation_amount_per_pop <- geo2$Donation_amount_accepted / geo2$pop
geo2$Donation_number_per_pop <- geo2$Number_donations_from_outside |> as.numeric() / geo2$pop
st_write(geo2, "furusato2017.geojson")
ここまでで、地理空間情報の準備は完了です。
leafletでコロプレス図を作成する
この地理空間情報のgeojsonを用いて、コロプレス図を作成していきます。地理空間情報の単位（市町村）が細かいので、静的な図では情報が読み取りにくくなります。ですので、leafletパッケージを用いてインタラクティブなコロプレス図を作成することにします。
leafletパッケージにはいろいろな関数が登録されており、いじるところはたくさんあるのですが、とりあえずはleafletのリファレンスのページに記載されているコードをそのまま使うことにします。Shinyに乗せることを考慮して、年と選択する列は別途オブジェクトとしておきます。
ただし、上で述べた通り、地理空間情報の無い行が存在するため、あらかじめst_is_empty関数で地理空間情報のない行を調べ、取り除いておきます。plot.sfは地理空間情報の無い行を勝手に取り除いてグラフにしてくれますが、leafletはきっちりエラーを出すので、行の取り除きは必須です。
options(scipen=100) # 桁の大きい数値を指数表示にしない
selected_col_number <- 20
year <- 2020
geofile <- paste0("furusato", year, ".geojson")
vec <-
c("pop", "pop_M", "pop_F", "pop_2015", "pop_change", "pop_change_rate",
"pop_density", "household", "Number_donations_accepted", "Donation_amount_accepted",
"Number_donations_from_outside", "Donation_amount_from_outside", "Costs_procure_gifts",
"Costs_return_gifts", "Costs_PR", "Costs_settlement", "Costs_Administration",
"Costs_others", "Total_costs", "Donation_amount_per_pop", "Donation_number_per_pop")
names(vec) <-
c("人口（2020年）", "男性人口", "女性人口", "人口（2015年）", "5年間の人口増減数",
"5年間の人口増減率", "人口密度", "世帯数", "受け入れた寄付件数", "受け入れた寄附金額",
"市町村外から受け入れた寄付件数", "市町村外から受け入れた寄付金額", "返礼品等の調達に係る費用",
"返礼品等の送付に係る費用", "広報に係る費用", "決済等に係る費用", "事務に係る費用",
"その他の費用", "費用合計", "一人当たりの市町村外から受け入れた寄付金額",
"一人当たりの市町村外から受け入れた寄付件数")
vec_nin <- c("pop", "pop_M", "pop_F", "pop_2015", "pop_change", "pop_density")
vec_ken <- c("Number_donations_accepted", "Number_donations_from_outside", "Donation_number_per_pop")
vec_yen <- c(
"Donation_amount_accepted", "Donation_amount_from_outside", "Costs_procure_gifts",
"Costs_return_gifts", "Costs_PR", "Costs_settlement", "Costs_Administration",
"Costs_others", "Donation_amount_per_pop")
geo <- st_read(geofile)
# st_is_empty(geo) # geometryが空（empty）の行がある
geo <- geo |> filter(!st_is_empty(geo)) # geometryにemptyがあるとleafletで取り扱いできない
col_selected <- vec[selected_col_number]
val <- geo[col_selected] |> st_drop_geometry() |> unlist()
if(col_selected %in% vec_nin){
char_val <-
paste0(format(val, big.mark = ",", scientific = F) |> str_trim(), "人")
} else if(col_selected %in% vec_ken){
char_val <-
paste0(format(val |> round(2), big.mark = ",", scientific = F) |> str_trim(), "件")
} else if(col_selected %in% vec_yen){
char_val <-
paste0(format(val |> round(2), big.mark = ",", scientific = F) |> str_trim(), "円")
} else if(col_selected == "pop_change_rate"){
char_val <-
paste0(val |> round(3), "%")
} else if(col_selected == "household"){
char_val <-
paste0(format(val, big.mark = ",", scientific = F) |> str_trim(), "軒")
}
bins <- quantile(val |> na.omit(), probs=seq(0, 1, by=0.1)) |> round(3)
pal <- colorBin("Spectral", domain = val, bins = bins, reverse=TRUE)
geo |>
leaflet() |>
addProviderTiles(providers$OpenStreetMap) |>
setView(137.5, 37.5, zoom = 6) |>
addPolygons(
fillColor = ~pal(val),
weight = 1,
opacity = 1,
color = "white",
dashArray = "3",
fillOpacity = 0.5,
highlightOptions = highlightOptions(
weight = 5,
color = "#666",
dashArray = "",
fillOpacity = 0.3,
bringToFront = TRUE),
label = map2_chr(geo$city, char_val, paste),
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")) |>
addLegend(pal = pal, values = ~val, opacity = 0.7, title = names(vec)[selected_col_number],
position = "bottomright")
選択した市町村の情報をポップアップで得られるようにするため、数値を処理（コンマと単位をつける）しています。また、色分けは10% quantileごとに分けられるようにしています（一人当たりの寄付件数では0付近のデータが極めて多いので、小数点以下3桁まで用いるようにしています）。ここまでできればShinyに乗せるのはそれほど難しくありません。
Shinyでコロプレス図を表示する
Shinyでは、上のコードをrenderLeaflet関数の中身としてserver側で呼び出し、ui側ではleafletOutput関数で結果を呼び出すだけです。ただし、地図データを大きく表示するためには少し工夫が必要です（このページを参照）。
div(leafletOutput("leafletPlot", height="100%"), style = "height: 73vh")
フッターをつけるために"height: 73vh"と少し値をいじっています。Shinyアプリはgithubで公開しているので、以下のコードで実行可能です。
if(!require(shiny)){install.packages("shiny")};runGitHub("furusato_map", "sb8001at")
うまくShinyが走ると以下の図のようなアプリが表示されます（その他の費用だけは0過剰データすぎてquantileに0が並んでしまうため、表示されません）。
感想
ちょうどsfやstarsパッケージの勉強をしていたので、データの取り扱いを学ぶのに良いデータだったと思います。
政府統計のデータは各種取り揃えられて公開されている状況になっており、地理データは上記の国土数値情報ダウンロードサイトからダウンロードできます。10～20年前にこのような図を作ろうと思えば各種データをあちこちに問い合わせ、専用のソフトウェアを準備してやっと作れるかどうか、インタラクティブなものにするのはまず無理といった状況だったと思います。
現在ではRではsfパッケージ、starsやterraなどのラスタデータを取り扱うパッケージ、leafletやtmapなどのインタラクティブ・静的マップの表示に関するパッケージが揃い、前述の通り政府統計へのアクセスも非常に簡単になり、地理空間情報というハードルが高そうなものも簡単に取り使うことができるようになっています。データいじりをするにはいい時代になったものです。
Shinyとして公開するのも、shiny::runGitHub関数で共有すればそれほど難しくありません（shinyapps.ioも使いやすいのですが、ちょっと計算が重そうな地理データなどは避けたほうがよいでしょう）。インタラクティブなグラフや地図を用いたアプリをgithub経由で共有することができます。Shinyの共有にはexe化やRportableと一緒に配布する方法もありますが、githubを使って共有するのが比較的簡単です。
ふるさと納税のデータに関しては、市町村ごとで力の入れようが異なるのがわかります。北海道の市町村は都心を除き一人当たり寄付金が多め、東北にもぽつぽつと一人当たり寄付金が多めの場所があります。一方で山間部の田舎（例えば奈良南部など）では寄付金が少なく、下位から数えたほうが圧倒的に早い結果になっています。おそらく返礼品を準備するだけの人もパワーも土地も無いからではないでしょうか。人口の多い都市でも、例えば泉佐野市（大阪）や高山市（岐阜）などは寄付金が多く、「人口の多いところから少ないところへ」お金が移転しているわけではないようです。また、返礼品の輸送や返礼品に関する広告等の料金もバカにならず、輸送・広告を（何なら返礼品の企画開発も）担っているのは東京などの大都市に本社機能のある大会社です。どちらかというとお金を集めやすい人・人口の多いところへお金が動いているような感があります。
どのような経緯でふるさと納税が作られたのかよくわかりませんが、国が税なり何なりでお金を集め、使うのであれば、「機会の平等」・「結果の平等」・「公共のため」のいずれかを達成するためである方がよいと思います。機会は平等ではなさそうですし、田舎からお金を巻き上げている要素もあるため結果の平等もありません。公共のためになっているのかも疑わしいです。正直あんまり良い制度ではないな、というのが解析結果に対する感想です。