いきものログ提供のデータをRに読み込んで地図上に可視化
Shinya Uryu
October 07, 2015
library(magrittr)
library(ggplot2)
library(viridis)
library(XML)
library(RCurl)
library(gridExtra)
library(sp)
library(readr)
library(dplyr)環境省が運営している「いきものログ」というwebサイトがある。身近な生物の目撃情報について共有する投稿型のシステムであり、気に入っている点の一つに、誰かが集めてくれた生物の情報を自分も利用できることがある。いきものログでは、投稿された生物データを.csvや.kml、.shpのファイルとしてダウンロードすることができる。いますぐ使い道が思いつくわけではないが、気になるサービスなのでそのデータをRで利用できるようにする場合のメモ。
今回は夏真っ盛りいつの間にか秋になっていた！ということでカブトムシについての情報を得てみたい。生物名「カブトムシ」で検索すると、747件のデータが該当した（2015年10月6日 現在）。データのダウンロードは手動で行う（Javascriptで制御しているみたいなのでRでダウンロードする方法がわからん）。
dir.create(path = "1510ikilog")
unzip(zipfile = "/Users/uri/Downloads/ikilog_csv20151006012449.zip",
exdir = "1510ikilog/")
list.files(path = "1510ikilog/")データの読み込み
いきものログのページでは、CSV形式でのダウンロードと書かれているが、実際にダウンロードされるのはTXT形式のファイルである。それはそれとして、ダウンロードしてきたファイルをRに読み込む手順。データフレームに変換し、必要な行および列に整形する。また、データの可視化のために位置情報を含んだデータ行のみを抽出する。
# 空のdata.frameオブジェクトを用意する
df_beetle <- data.frame()
read_delim(file = "1510ikilog/ikilog_20151006012449.txt",
locale = locale(encoding = "CP932"),
delim = "\t") %>%
# readr パッケージを使用しない場合...
# read.table("1510ikilog/ikilog_20151006012449.txt", fileEncoding = "cp932",
# sep = "\t", stringsAsFactors = FALSE,
# na.strings = "",
# header = TRUE, fill = TRUE) %>%
.[-1:-4, ] %>%
dplyr::select(., c(1:32)) %>%
dplyr::filter(., !is.na(life_darwincore_id), 英語項目名 == "日本語") %>% {
nrow(.) %>% print(.) # 全データが取得できているか確認
# 位置座標の情報が含まれているデータのみを抽出
df_beetle <<- dplyr::rename(.,
longitude = decimalLongitude,
latitude = decimalLatitude) %>%
dplyr::filter(., !is.na(longitude), !is.na(latitude)) %>%
dplyr::mutate_each(funs(as.numeric), longitude, latitude, individualCount) %>%
.[-ncol(.)] %T>% {
dplyr::glimpse(.)
}
}## [1] 747
## Observations: 727
## Variables: 31
## $ 英語項目名 (chr) "日本語", "日本語", "日本語", "日本語", "日本語", "日本語", "日本語",...
## $ life_darwincore_id (chr) "8342145", "8245591", "8214856", "8213250",...
## $ registrantID (chr) "クリス", "花虫山海", "chasuke", "クマゲラ", "yasu", "...
## $ year (chr) "2015", "2015", "2015", "2015", "2015", "20...
## $ month (chr) "7", "8", "8", "8", "7", "7", "7", "7", "7"...
## $ day (chr) "28", "23", "19", "4", "13", "15", "12", "3...
## $ eventRemarks (chr) "ふくおか生きもの見つけ隊調査報告", "ふくおか生きもの見つけ隊調査報告", NA,...
## $ scientificName (chr) "Allomyrina dichotoma dichotoma", "Allomyri...
## $ Japanese_Name (chr) "カブトムシ", "カブトムシ", "カブトムシ", "カブトムシ", "カブトムシ"...
## $ open_status (chr) "0", "0", "0", "0", "0", "0", "0", "0", "0"...
## $ comment (chr) "立派なオスですが、上翅が傷ついていました。", "朝、庭で見つけました。", NA,...
## $ picture1 (chr) "8342145-1.JPG", "8245591-1.JPG", "8214856-...
## $ caption1 (chr) "かぶと", "元気なカブトムシ", "カブトムシ1", "カブトムシ雌", "カブト...
## $ picture2 (chr) NA, NA, "8214856-2.JPG", NA, "7920159-2.JPG...
## $ caption2 (chr) NA, NA, "カブトムシ2", NA, "カブトムシ", NA, "コガネムシ科"...
## $ picture3 (chr) NA, NA, "8214856-3.JPG", NA, "7920159-3.JPG...
## $ caption3 (chr) NA, NA, "カブトムシ3", NA, "カブトムシ", NA, NA, NA, ...
## $ movie_url (chr) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ latitude (dbl) 33.51783, 33.61917, 35.45622, 35.94424, 36....
## $ longitude (dbl) 130.4902, 131.1378, 138.7633, 139.7243, 140...
## $ mesh1 (chr) "5030", "5031", "5338", "5339", "5440", "53...
## $ mesh20km (chr) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ mesh2 (chr) "503023", "503131", "533816", "533975", "54...
## $ mesh3 (chr) "50302329", "50313141", "53381641", "533975...
## $ mesh5k (chr) "5030232", "5031311", "5338161", "5339752",...
## $ stateProvince (chr) "福岡県", "福岡県", "山梨県", "埼玉県", "茨城県", "東京都", "...
## $ county (chr) "大野城市", "豊前市", "富士吉田市", "さいたま市岩槻区", "水戸市", ...
## $ municipality (chr) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ locality (chr) NA, "福岡県豊前市大字赤熊", NA, NA, NA, NA, NA, NA, N...
## $ verbatimLocality (chr) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ individualCount (dbl) 1, 1, 1, 1, 1, 1, 1, 1, 10, 3, 1, 2, 1, 6, ...最終的な件数は727件になった
地図上にプロット
いきものログのページでできる可視化をやってみる。
日本地図の用意
まずは元になる日本地図を用意する。今回はGitHubにてdataofjapanのlandリポジトリにて公開されている地球地図日本( http://www.gsi.go.jp/kankyochiri/gm_jpn.html )提供のShapefileを変換したJSON形式のファイルを使用させてもらう。GeoJsonだとちょっと重いのでTopoJsonにするのが良い。
jp <- rgdal::readOGR(dsn = "https://raw.githubusercontent.com/dataofjapan/land/master/japan.geojson",
layer = "OGRGeoJSON",
stringsAsFactors = FALSE)
# jp <- geojsonio::topojson_read(x = "https://raw.githubusercontent.com/dataofjapan/land/master/japan.topojson")
jp %<>% fortify()
jp_map <- ggplot() +
geom_map(data = jp,
map = jp,
aes(x = long, y = lat, map_id = id),
fill = "white", color = "black",
size = 0.5) +
coord_map() +
ggthemes::theme_map()class(jp_map)## [1] "gg" "ggplot"こちらのggplotオブジェクトにデータを付加していく。
緯度経度情報に基づくポイントデータの表示
jp_map +
geom_point(data = df_beetle,
aes(x = longitude, y = latitude),
color = "#ffff33", fill = "#ff7f00",
shape = 21, alpha = 1/2)
都道府県ごとの報告個体数集計
# 都道府県ごとの報告個体数を集計する
df_beetle %>% group_by(stateProvince) %>%
summarise(count = sum(individualCount)) %>%
rename(Kanji = stateProvince) -> df_pref
# 都道府県別に図示するための処理
"https://en.wikipedia.org/wiki/Prefectures_of_Japan" %>%
getURL() %>%
readHTMLTable(which = 4) %>%
select(Kanji, ISO) %>%
dplyr::mutate(id = as.numeric(gsub(pattern = "JP-", x = ISO, replacement = ""))) %>%
left_join(., df_pref, by = "Kanji") -> df_pref## Warning in left_join_impl(x, y, by$x, by$y): joining character vector and
## factor, coercing into character vectorjp_map +
geom_map(data = df_pref,
map = jp,
aes(fill = count, map_id = id, stat = "identity")) +
scale_fill_viridis(alpha = 0.8)
市区町村別の図示は細かくなりすぎるのでちょっとパス。
メッシュ
80km四方、20km四方、10km四方、5km四方、1km四方と表示できるが、80km四方（地域一次メッシュ）のみ。国土交通省 国土地理院 ２０万分１地勢図図を元にして、新たに、地域一次メッシュのGeoJsonを用意した。
jp80 <- geojsonio::topojson_read(x = "https://gist.githubusercontent.com/uribo/b09d642351c03dd975aa/raw/6cd2f7e70922b24f89b442c96c001dc3ab794317/japan_one_twenty_map.topojson") %>%
fortify() %>%
group_by(group) %>%
mutate(new.id = paste0(round(as.numeric(min(lat) * 1.5), 2), as.numeric(substr(min(long), 2, 3)))) %>% # 一次メッシュの算出
rename(old = id, id = new.id)
jp_map80 <- ggplot() +
geom_map(data = jp80,
map = jp80,
aes(x = long, y = lat, map_id = id),
fill = "white", color = "black",
size = 0.5) +
coord_map() +
ggthemes::theme_map()
df_beetle %>% group_by(mesh1) %>%
summarise(count = sum(individualCount)) %>%
rename(id = mesh1) -> df_mesh1jp_map80 +
geom_map(data = df_mesh1,
map = jp80,
aes(fill = count, map_id = id, stat = "identity")) +
scale_fill_viridis()
報告個体数の経年変化
df_beetle %>% dplyr::group_by(year) %>%
summarise(n = sum(individualCount)) %$% {
ggplot(, aes(year, n)) +
geom_bar(stat = "identity", fill = "#71D91080") +
ggthemes::theme_hc()
}
おまけ
こんな感じで図示するのもいかがか。
# ポイントを適当なまとまりで表示
p1 <- jp_map +
stat_binhex(data = df_beetle,
bins = 40,
alpha = 0.80,
aes(x = longitude, y = latitude,
fill = cut(..count.., c(1, 4, 9, 49, 99, 499, 999, 4999, 5000, Inf)))) +
scale_fill_brewer(
palette = "OrRd",
labels = c("1", "2-4", "5-9", "10-49", "50-99",
"100-499", "500-999", "1000-4999", "5000+"))
# 報告個体数でポイントの大きさを重み付け
p2 <- jp_map +
geom_point(data = df_beetle,
aes(x = longitude, y = latitude),
color = "#ffff33", fill = "#ff7f00",
shape = 21, alpha = 1/2, size = df_beetle$individualCount * 3)
gridExtra::grid.arrange(p1, p2, ncol = 2)
出典
これらの地図は、環境省いきものログ（調査名：いきものみっけ）2015年10月6日時点の情報を使用し、瓜生真也(@u_ribo)作成（複製）したものである。