Rを使って鳥類などの分布図を作成する（その1：緯度経度と市町村別マップ）
（広報担当：上沖正欣）
当会では2024年中の愛媛県鳥類目録の出版を目指しています。利用したデータは文献情報や会員から毎月収集している野鳥情報で、1954年から約360種類18万件以上もあり、種名だけではなく年月日や市区町村などの地点、3次メッシュ情報などが含まれます。「目録」としては種類がリストされてればそれでよいのですが、せっかくなので経年変化や記録の多い地点を可視化したいと思いました。
地図を扱う場合はQGISなどGISソフトのほうがよいのかもしれませんが、他にもデータ整形やグラフ一覧出力などもしたかったので、そういったデータ処理が得意とされるRを使いました（勿論、GISや流行りのPythonでも同じことが実現できます）。生物系で分布情報をRで可視化している例、さらに地域メッシュを利用した例となると意外と少なかったので、需要があるかどうか分かりませんが、ここで紹介していこうと思います（正直なところ、ChatGPTやGoogle GeminiなどのAIに「このデータを地図に可視化するコードを作って」と質問したほうが早いです。実際、サンプルデータを渡してこんなコード書いて、と指示すると、この記事で紹介しているものとほぼ同じコードを出力してくれました・・・）。
Rで野鳥情報を可視化したい
手元にエクセルデータはあって数を集計したりグラフを描くことはできるけど、分布図を描く場合にはどうすれば？と思ったことはないでしょうか。可視化することで数字や文字のデータでは分からなかった情報を読み取ることができるようになります。全国鳥類繁殖分布調査の報告書で使われているような分布図を、自分のデータで描けたら素敵ですよね。
これはRで描いた1993年から2022年までの県内のリュウキュウサンショウクイの報告地点です（会員が住んでいる場所に偏っていることに注意）。左端は緯度経度情報、真ん中と右端はそれぞれ各市町と2次メッシュごとの記録数を、グラデーションで表現しています。真ん中の地図はコロプレスマップといって、コロナの感染者数を地域ごとに可視化するためにもよく使われていたもので馴染みがあるかと思います。
こうした地図を描くには、下図のような形式のデータが必要です。種名の列、および位置情報である市町orメッシュor緯度経度の列は必須です。個体数や日付など他の列は今回関係ないので、あってもなくてもOKです。
ご自身がお持ちのデータでもよいですが、データ型や欠損値などの問題があるとエラーが出てしまうので、まずはこちらのサンプルデータをダウンロードして利用してみてください（ヘッダーの列名は shumei/species などと英字にしたほうがよいのですが、今回は分かりやすくするため日本語にしています）。
以下長々と説明していますが（1）ソフトのインストール（2）データファイルのダウンロード（3）フォルダ作成（4）コピペして実行、とたったの簡単4ステップです。
1. 必要な事前作業：ソフトウェアのインストールと作業用フォルダの作成
Rって、研究者が統計処理に使うものでしょ？難しそう、と思われるかもしれません（私も苦手で使いこなしているわけではありません）。でも、地図を描くだけならモデルを組んだり難しい計算をするわけではないので、意外と簡単です。以下のRを使った操作は、野鳥情報を管理されているような、エクセルを普段使う方でPC操作は結構できる、でもプログラミングの経験はなくRも使ったことがない、という方向けに解説していきます（すごーくニッチですが）。
まず、事前にR（あーる）とRStudio（あーるすたじお）という2つのソフトウェア、また市町村の境界線データのダウンロードが必要です。
Rのダウンロードはこちら Windows版 Mac版
RStudioのダウンロードはこちら RStudio
これらのインストール方法は、一般的なソフトウェアと同じなのでここでは省略します（参考リンク）。インストール中の設定はデフォルトのままで、変更しないでください。
市町村の境界線データ 国土交通省の行政区域データ
リンク先の地図から愛媛県のZIPファイルをダウンロードしてください。ZIPファイルを解凍すると、地図の境界線を描くのに必要な一連のファイルが取得できます（jpndistrictというパッケージを使えば別途ダウンロードせず済むのですが、CRANにないのと市町村名に郡が含まれていて処理が煩雑になってしまうのでここでは使いません。また全国データをダウンロードして都道府県を指定することもでき、こちらの方がスマートかもしれませんがデータが重いためか私の環境ではうまくいきませんでした）。愛媛ではなく違う都道府県のデータをダウンロードし、上記サンプルデータを使う場合は、サンプルデータの市町村名を、ダウンロードした行政区域の市町村名に変更してください。
フォルダ名が「N03-20210101_39_GML」と長いので簡略化のため「愛媛」に変更しておきます。
そうしたら、この行政区域データが入ったフォルダと観察記録のエクセルファイルを、Rで使うことになるフォルダを作って、そこにコピーします。今回はCドライブの直下にRdataというフォルダを作成することにします。作業用フォルダの場所を別に設定したい場合は、ご自身で変更してください。
上図のようにCドライブにRdataというフォルダを新規作成し、その中に先ほどの行政区域データが入った「愛媛」というフォルダと、観察データのエクセルファイル（data.xlsx）を入れています。繰り返しますが、フォルダ名やファイル名はできれば英字の方が好ましいです。
2. RStudioの事前設定：文字コードを変更する
ここまでの事前作業が終わったら、RStudioを起動してください（PCのユーザー名が日本語の場合RStudioが動かないことがあるようです。ここでは対処法には触れませんが、起動できなかった場合はユーザー名を要確認）。起動できれば、もう全行程の8割くらいの作業は終了しています。RStudioを開くと、下記のようなウィンドウが表示されるかと思います（英語ですが気にしない）。
文字コードの変更：上メニューのToolから、Global Options>Code と進んで文字コードをUTF-8に変更して下さい（参考リンク）。文字化けしてしまうので、日本語を扱う上でこの設定は必須です。
3. 地図を描くためのコードをコピペして実行
では、いよいよ地図を描いていきます。と言ってもコピペして一瞬で終わります。
RStudioの上部メニューFile>New File>R Scriptを選択して、表示される右上のUntitled1タブの空白のところ（スクリプトを書くメモ帳のようなものです）に下記のコードを全てコピペして、最初の行からコマンドを実行（全選択して、Ctl+Enterもしくはコードを貼り付けた上右あたりに「Run」があるのでそれをクリック）すれば、ダウンロードした都道府県の地図が出力できるはずです。
#必要なパッケージをインストールします。
#この3行目は、初回実行時のみでOKです。
install.packages(c("openxlsx", "sf", "rmapshaper", "dplyr","ggplot2"))
#インストールしたパッケージを読み込みます。
library(openxlsx)
library(sf)
library(rmapshaper)
library(dplyr)
library(ggplot2)
#作業用フォルダの場所を指定します。今回はCドライブ直下に作ったRdataです。
#違う場所にデータを格納した場合は、ご自身の設定に合わせてフォルダの場所を指定してください。
setwd("C:/Rdata")
#Rdataフォルダの中にある観察データのエクセルを読み込み、kansatsuという名前にします。
#ファイル名はご自身のものに変更してください。
kansatsu <- read.xlsx("data.xlsx")
#「愛媛」フォルダの中にある行政区域データを読み込み、map_shapeという名前にします。
#フォルダ名はご自身のものに変更してください。
map_shape <- read_sf("愛媛", options = "ENCODING=SHIFT-JIS") %>%
#N03_004列に市町村名が入力されているので、分かりやすくcityと名前を変更します。
rename(city = N03_004)
#city名を基準に統合し、map_shapeのデータを「mapdayo」にまとめます。
mapdayo <- aggregate(map_shape, list(map_shape$city), unique) %>%
#必要なcity列とgeometry列のデータのみを選択します。
select(city, geometry) %>%
#rmapshaperの機能で境界線を単純化し、地図の出力を高速化します。値は自由に変えてみてください。
ms_simplify(keep = 0.005)
#地図が表示されれば成功！
plot(mapdayo)ぱっと見、なにやら複雑なことをしているように見えますが、単にファイルを読み込んでその名前をXとしてね、と簡単な指示をしているだけで、計算は全くしていません。17行目などにある記号「<-」は矢印（= イコールとほぼ同義）です。X <- Y とあれば、XにYを代入せよ、ということです。
#の行は単なるコードの説明しているコメントなので、改行など除くと実際には11行のとても短いコードです。3行目でデータ処理に必要なパッケージ（追加機能）をインストールしています。実行するとConsoleのところに何やら赤字で色々表示されて待たされると思いますが、エラーではないので大丈夫です。今回必要なのは下記の5つのパッケージです。RStudioのTools>Install Packagesからでもインストールできます。
（1）エクセルを読み込む openxlsx
（2）シェープファイルを読み込む sf
（3）地図の境界線を単純化する rmapshaper
（4）色々便利機能が使える dplyr
（5）グラフ表示をおこなう ggplot2
もし左下のConsole画面に赤文字で何やらエラーが表示されたら、コピペ間違いか、フォルダの名前違いか、ファイル名が違っている可能性が高いです。そのエラーメッセージをそのままコピーしてChatGPTやGoogle Geminiに何が原因か聞いてみましょう。どこが間違っているか、どう直せばいいか、きっと人間より的確に教えてくれます。
ちなみに26行目のaggregateの処理はsummarizeを使った下記コードでも同じ結果が得られますが、処理が遅かったのでaggregateを使いました。
mapdaoyo <- map_shape %>%
group_by(city) %>%
summarise(geometry = st_union(geometry))境界線のデータは複雑な処理が必要なので、左下のPlotsのところに地図が表示されるまで少し時間がかかりますが、Consoleにエラーが出ていなければ表示されるはずですので、しばらく待ちましょう。
4. 地図に観察データを反映させる
ではいよいよこの地図に観察データを反映させていきます。先ほどのコードの下に、以下のコードをコピペしてください。こちらも、一行ずつコメントを書いていますが、12行ほどの短いコードです。
#先ほどのmapdayoとkansatsuのデータを結合し、それをkansatsu_mapという名前のデータにします。
kansatsu_map <- left_join(kansatsu, mapdayo, by = c("市町" = "city"))
#kansatsu_mapから地図を描画します。
kansatsu_map %>%
#特定の種類のみ表示したい場合はフィルターできます。ここでは全種表示したいので実行しません。
#filter(grepl("カッパ",種名)) %>%
#種名と市町の列でグループ化します（これはサンプルデータの列名。適宜変更してください）
group_by(種名,市町)%>%
#種ごとに、各市町の記録が何回あるかカウントし、新たなデータ列（市町別記録数）を追加します。
mutate(市町別記録数 = n()) %>%
#インストールしたパッケージggplotで地図を描画します。
ggplot() +
#mapdayoのデータで地図を描画します。地図は白色にします。
geom_sf(data = mapdayo, fill = "white") +
#先ほど追加した市町別記録数の列のデータを基に色付け(fill)をします。
geom_sf(aes(fill = 市町別記録数, geometry = geometry)) +
#グラデーションで表現。少なければ黄色（yellow）多ければ深緑（#009933）にする。
#色はカラーコードで細かく指定可能です。
scale_fill_continuous(low = "yellow", high = "#009933") +
#サンプルデータにある緯度経度の情報を表示。実際の経度や緯度の列名を指定してください。
#Xが経度、Yが緯度になります。逆だと変な図が出力されるので注意。
geom_point(aes(x = 経度, y = 緯度),
#地図上の点のサイズと色、重なりを表すため透過率を指定。数字や色は好みの値に変更できます。
size = 2, color = "black", alpha = 0.25) +
#複数種のデータがあるので、種名ごとに地図を分けて表示します。
facet_wrap(~種名) +
#背景に緯度経度が表示されるので、すっきりさせるために白紙にします。
theme_void()先ほどと同じように全て選択して実行すると、このようにサンプルデータにあるカッパ・コダマ・テングの分布図が表示されると思います（これはダミーデータで、実際の分布とは異なります）。
この地図は15行目のgeom_sf（愛媛県全域のデータ）と17行目のgeom_sf（記録のある市町のデータ）、23行目のgeom_point（緯度経度のデータ）の三つのレイヤーが重なって表示されています。なので、緯度経度の点のみを表示したい場合は、上のコードの17-20行目の先頭に#を付けて、コメントアウトしてください。行政区域の色分けデータのみ表示したい場合は、23-25行目をコメントアウトしてください。ちなみに、市町村ではなく全国の都道府県ごとの色分けであれば、エクセルだけで簡単にコロプレスマップが作れます。
どうでしょうか。ソフトウェアをインストールしてコピペするだけで分布図が簡単に出力できた（はず）と思います。次回は、これを応用して、地域メッシュを使った方法と、年代別の地図を表示する方法を解説します。
緯度経度・市町村別の観察データから地図を表示するR Script
全てのコードをまとめ最低限のコメントのみ付けたものは以下になります（GitHubで全3回分のまとめコードを公開しています）。
#必要なパッケージをインストール（初回実行時のみ）
install.packages(c("openxlsx", "sf", "rmapshaper", "dplyr","ggplot2"))
#インストールしたパッケージを読み込み
library(openxlsx)
library(sf)
library(rmapshaper)
library(dplyr)
library(ggplot2)
#データを入れたフォルダの絶対パスを指定
setwd("C:/Rdata")
#Rdataにある観察データのエクセルを読み込み
kansatsu <- read.xlsx("data.xlsx")
#行政区域のシェープファイルをRdataの「愛媛」フォルダから読み込み列名変更
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)
#地図情報と観察データを結合
kansatsu_map <- left_join(kansatsu, mapdayo, by = c("市町" = "city"))
#結合データを利用して、地図を描画
kansatsu_map %>%
group_by(種名,市町) %>%
mutate(市町別記録数 = n()) %>%
ggplot() +
geom_sf(data = mapdayo) +
geom_sf(aes(fill = 市町別記録数, geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933") +
geom_point(aes(x = 経度, y = 緯度), size=2, color = "black", alpha = 0.25) +
facet_wrap(~種名) +
theme_void()最後に：会員から収集した生物データは広く公開することが大切
野鳥の会だけでなく、各自然団体では独自の記録表を作って会員から情報収集している所もあると思います。しかし、そんな貴重なデータがあっても、例え会報やウェブサイトで公開されていたとしても、誰にも見られることなく情報が埋まってしまえば、その地域の生物多様性が正しく評価されないことにつながってしまう可能性があります。
会員から収集したデータを有効活用する最善の方法は、より大規模な生物データベース（DB）に登録することで、利用される機会をなるべく増やすことです。例えば生物全般では環境省の生きものログや米国のNPO団体が運営するiNaturalist、鳥類ではバードリサーチのフィールドノートやeBirdが挙げられます。
人々がスマホを持つようになった2000年以降のここ10-20年の間にこうした大規模なDBが構築されているので、積極的に利用することをお勧めします。当会も将来的には独自に収集している野鳥情報を無くし、こうした汎用性の高いDBを会員に直接利用してもらうことを考えています。
Leave a comment