Rを使って鳥類などの分布図を作成する（その2：地域メッシュと年別マップ）
前回「Rを使って鳥類などの分布図を作成する（その1：緯度経度と市町村別マップ）」としてコードをコピペするだけで地図を描画できる方法を紹介しました。今回はその続きでRを使った地域メッシュと年別の地図描画を紹介します。
左の緯度経度を示した地図は点データなので高解像度ですが、どのくらいの数が集まっているのか分かりにくく、また緯度経度の数値は少しでも場所がズレると変わるため、データ集計が煩雑になってしまいます。
真ん中の市町村別の地図では、面積がバラバラになるので、単純なデータ比較ができません。また生物は行政区域を基準に分布しているわけではないので、どんな環境にいるのかが分かりにくくなります。
しかし、右の地域メッシュ地図（英語ではmeshではなくgrid mapなどと言う）は決められた範囲の面データであるため（図は10km四方の2次メッシュ。3次メッシュなら1km四方）、一定範囲内の集計をしたい場合や、異なる場所間の比較が容易になります。
以下ではそんな便利な地域メッシュを地図に描画する方法を紹介します。RとRStudioはインストール済み、また前回の愛媛県の（もしくはご自身の地域の）サンプルデータを使った妖怪3種の緯度経度と地域別の分布図は既に描画できているものとします。データは前回と同様のサンプルデータを使います。
メッシュを描画する記述を追記する
前回の「3. 地図を描くためのコードをコピペして実行」の地図を表示するコードまでは一緒です。
library(openxlsx)
library(sf)
library(rmapshaper)
library(dplyr)
library(ggplot2)
setwd("C:/Rdata")
kansatsu <- read.xlsx("data.xlsx")
map_shape <- read_sf("愛媛", options = "ENCODING=SHIFT-JIS") %>%
rename(city = N03_004)
mapdayo <- aggregate(map_shape, list(map_shape$city), unique) %>%
select(city, geometry) %>%
ms_simplify(keep = 0.005)
plot(mapdayo)この後に、新たに以下のコードをコピペして全て実行してください。前回と異なる部分にだけコメントを付けています。
#jpgridパッケージの方が高速で海上のメッシュも描画できますが、今回はこちらを使います。
#前回同様install.packagesの行は初回実行時のみでOKです。
install.packages("jpmesh")
library(jpmesh)
#2次メッシュ（10km四方）のデータを作成します。codeの数字は都道府県コードです。
#愛媛にしたいので38、メッシュサイズ(km)は10を指定します。
mesh10map <- administration_mesh(code = 38, to_mesh_size = 10)
#サンプルデータのメッシュは数値、mesh10mapのメッシュは文字列なので結合できるよう変換します。
kansatsu$メッシュ <- as.character(kansatsu$メッシュ)
#mesh10mapとkansatsuデータを、メッシュ番号を基準に結合してkansatsu_meshを作ります。
kansatsu_mesh <- left_join(kansatsu,mesh10map, by = c("メッシュ" = "meshcode"))
#kansatsu_meshから地図を描画します。
kansatsu_mesh %>%
#種名とメッシュの列でグループ化（これは元データの列名。適宜変更してください）
group_by(種名,メッシュ)%>%
#種ごとに、各メッシュの記録が何回あるかRでカウントし、新たなデータ列（メッシュ数）を追加します。
mutate(メッシュ数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo, fill = "white") +
#mesh10mapのメッシュを描画します。alpha=0で透過させ、線はグレーにします。
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
#先ほど追加したメッシュ数の列のデータを基に色付け(fill)をします。
geom_sf(aes(fill = メッシュ数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
#種名ごとに地図を分け、1列で表示します
facet_wrap(~ 種名, nrow = 1) +
theme_void()都道府県コードは国土交通省のページを参照してください。
どうでしょう。前回のコードを少し変えるだけで、地域メッシュを使った分布図がコピペで簡単に描画できたと思います。緯度経度や市町村別のコロプレスマップよりも、どの場所に多いかが分かりやすいですね。
サンプルデータは2次メッシュ（数字6桁）ですが、もちろん3次メッシュ（8桁）でもできます（3次メッシュの場合は数が多いので描画にちょっと時間がかかります）。
地図を種別かつ年別に表示する
分布が年々拡大しているor縮小している種の場合、分布変化を年毎に見たい場合があると思います。その場合は、先ほどのコードの19行目のgroup_byと30行目のfacet_wrapの記述に年（元データの年情報が入っている列名）を以下のように追加すればよいだけです。
kansatsu_mesh %>%
#種名+メッシュ+年でグループ化してカウント
group_by(種名, メッシュ, 年) %>%
mutate(メッシュ数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo, fill = "white") +
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
geom_sf(aes(fill = メッシュ数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
#facet_wrap(年 ~ 種名)でもOK。グラフのタイトル順（優先カテゴリ）が変わります。
#group_byで細かく年+月別にカウントした場合は 種名~年+月 などと + を付けて条件を増やせます。
facet_wrap(種名 ~ 年) +
theme_void()
ただこれだけだと全ての種で年毎の地図が表示されて多くなるので、種名を特定の種のみにフィルターするコードをkansatu_mesh %>% の後に追加して、行数も調整してみます。
kansatsu_mesh %>%
#kansatsu_meshのデータからカッパのみ抽出します。
filter(grepl("カッパ",種名)) %>%
group_by(種名,メッシュ,年) %>%
mutate(メッシュ数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo, fill = "white") +
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
geom_sf(aes(fill = メッシュ数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
#nrowを指定して2行表示にします。ncolを指定すると列数も指定できます。
facet_wrap(種名 ~ 年, nrow = 2) +
theme_void()
これで、カッパだけの年毎の分布図が作成できました。
応用して、group_byとfacet_wrapの箇所を月や市町に変えれば、月ごと・市町村ごとの分布図も作成できます。記録されたメッシュの場所だけ表示し、県全体の灰色の空白枠を非表示にしたい場合は geom_sf(data = mesh10map, alpha = 0, color = “grey”) + の行をコメントアウトします。
ここでは種類をフィルターしましたが、数百種類の結果を一気に出力したい場合もあると思います。次回は、ggforceパッケージを利用した複数ページの一括出力を紹介しようと思います。
地域メッシュと年別マップを表示するRスクリプトまとめ
library(openxlsx)
library(sf)
library(rmapshaper)
library(dplyr)
library(ggplot2)
setwd("C:/Rdata")
kansatsu <- read.xlsx("data.xlsx")
map_shape <- read_sf("愛媛", options = "ENCODING=SHIFT-JIS") %>%
rename(city = N03_004)
mapdayo <- aggregate(map_shape, list(map_shape$city), unique) %>%
select(city, geometry) %>%
ms_simplify(keep = 0.005)
plot(mapdayo)
#ここまで前回のコードと同じ。地図が表示されればOK
#jpmeshパッケージをインストール。初回実行時のみでOK
install.packages("jpmesh") library(jpmesh)
#パッケージを読み込み
library(jpmesh)
#愛媛県の2次メッシュ（10km四方）のデータを作成
mesh10map <- administration_mesh(code = 38, to_mesh_size = 10)
#データ型を変換して結合しkansatsu_meshというデータを作成
kansatsu$メッシュ <- as.character(kansatsu$メッシュ)
kansatsu_mesh <- left_join(kansatsu,mesh10map, by=c("メッシュ"="meshcode"))
#種別のメッシュ地図を描画する場合はここから48行目までを実行
kansatsu_mesh %>%
#種名とメッシュの列でグループ化
group_by(種名,メッシュ)%>%
#種ごとに、各メッシュの記録が何回あるかRでカウントし、新たなデータ列（メッシュ数）を追加します。
mutate(メッシュ数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo, fill = "white") +
#mesh10mapのメッシュを描画します。alpha=0で透過させ、線はグレーにします。
geom_sf(data=mesh10map,alpha=0,color="grey")+
#先ほど追加したメッシュ数の列のデータを基に色付け(fill)をします。
geom_sf(aes(fill = メッシュ数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
#種名ごとに地図を分け、1列で表示します
facet_wrap(~種名, nrow = 1) +
theme_void()
#カッパのみ、年別のメッシュ地図を描画する場合は、31行目の後に以下を実行。
kansatsu_mesh %>%
filter(grepl("カッパ", 種名)) %>%
group_by(種名, メッシュ, 年) %>%
mutate(メッシュ数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo, fill = "white") +
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
geom_sf(aes(fill = メッシュ数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
facet_wrap(種名 ~ 年, nrow = 2) +
theme_void()
