R言語 Advent Calendar 2023の16日目の記事です。
背景
去年くらいから札幌の都市部でも羆害事件や目撃報告が増えており、ヒグマがことさら身近な印象です。公開されている情報をRで収集・整理・可視化してみます。
webスクレイピング
rvestパッケージを使って札幌市ヒグマ出没情報のページから出没情報データをスクレイピングしていきます。表形式のデータを抜き出します。
library(tidyverse)
library(rvest)
higuma_url <- "https://www.city.sapporo.jp/kurashi/animal/choju/kuma/syutsubotsu/"
page_content <- rvest::read_html(higuma_url)
higuma_tables <-
page_content |>
rvest::html_elements("table")
higuma_tables
このページには表がふたつあり、一方が今年度の出没情報で他方が昨年度1–3月の出没情報です。
{xml_nodeset (2)}
[1] <table width="100%"><tbody>\n<tr>\n<th>No</th>\n<th>日時</th ...
[2] <table width="100%"><tbody>\n<tr>\n<th>No</th>\n<th>日時</th ...
今年度の情報のみを取り出します。中身を確認すると、以下のような形式になっています。
higuma_table <-
higuma_tables[[1]] |>
html_table()
head(higuma_table)
| No | 日時 | 場所 | 地図 | 内容 |
|---|---|---|---|---|
| 227 | 2023年12月6日(水曜日) | 南区定山渓無番地(定山渓湯の沢林道入口から約750mの地点) | 地図\n(グーグルマップへのリンク） | 足跡を確認 |
| 226 | 2023年12月4日(月曜日) | 南区藤野902番1付近(藤野聖山園敷地内) | 地図\n(グーグルマップへのリンク） | 足跡を確認 |
| 225 | 2023年12月3日(日曜日)\n7時57分 | 南区簾舞3条1丁目3番付近 | 地図\n(グーグルマップへのリンク） | ヒグマらしき動物を目撃 |
| 224 | 2023年12月2日(土曜日)\n18時45分 | 南区石山610番5付近(石山穴の沢線の道路上) | 地図\n(グーグルマップへのリンク） | ヒグマを目撃 |
| 223 | 2023年11月30日(木曜日)\n6時0分 | 南区簾舞1条4丁目2番付近(簾舞小北側道路上) | 地図\n(グーグルマップへのリンク） | ヒグマらしき動物を目撃 |
| 222 | 2023年11月28日(火曜日)\n21時55分 | 中央区南28条西13丁目1番付近(藻岩山斜面) | 地図\n(グーグルマップへのリンク） | ヒグマを目撃 |
場所列から地図を作ることもできそうですが、山中の地点が多数あることと表記に揺らぎがあることから厄介そうです。代わりに地図列のgoogle mapへのリンクに埋め込まれた緯度経度を取得することにします。以下のようにhtmlタグのクラスを指定し、urlリンクを取り出します。
map_links <-
higuma_tables[[1]] |>
rvest::html_elements("td") |>
rvest::html_elements("a") |>
as.character()
head(map_links)
#> [1] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=42.980178%2C141.071612&z=17\">地図</a>"
#> [2] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=42.930313%2C141.300129&z=17\">地図</a>"
#> [3] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=42.960721%2C141.27308&z=17\">地図</a>"
#> [4] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=42.9491537%2C141.3249449&z=17\">地図</a>"
#> [5] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=42.964626%2C141.265745&z=17\">地図</a>"
#> [6] "<a href=\"https://www.google.com/maps/d/viewer?mid=z1tD77I5wNwc.kjf5OMKDCpLY&t=p&ll=43.022111%2C141.342557&z=17\">地図</a>"
正規表現で緯度・経度に相当する部分のみを取り出します。文字列ベクトルxから緯度・経度のdata.frameを作成する関数を用意し、一括処理しました。
parse_geo <-
function(x){
strngir::str_extract(x, regex("(?<=ll=)(.*)(?=&)")) |>
strngir::str_split("%2C", simplify = TRUE) |>
data.frame() |>
dplyr::transmute(lat = as.numeric(X1), lon = as.numeric(X2))
}
higuma_location <-
parse_geo(map_links)
head(higuma_location)
#> lat lon
#> 1 42.98018 141.0716
#> 2 42.93031 141.3001
#> 3 42.96072 141.2731
#> 4 42.94915 141.3249
#> 5 42.96463 141.2657
#> 6 43.02211 141.3426
可視化に向けてデータを整えておきます。時間列を日付データに変換し、目撃情報に応じて危険度を場合分けし、マップ可視化時のポップアップ用の列を作りました。子グマと一緒なら極めて危険、ヒグマが視認されたら危険、痕跡あるいはヒグマらしきものが確認された場合はやや危険、の３水準としています。
higuma_data <-
higuma_table |>
dplyr::transmute(No,
date = lubridate::ymd(stringr::str_sub(日時, 1, 12)),
risk = dplyr::case_when(
stringr::str_detect(`内容`, "親子") ~ "extreme",
stringr::str_detect(`内容`, "目撃|カメラ") & stringr::str_detect(`内容`, "らしき", negate = TRUE) ~ "high",
TRUE ~ "medium"),
popup = paste0(`日時`, "<br>", `内容`)) |>
dplyr::bind_cols(higuma_location)
head(higuma_data)
#> # A tibble: 6 × 6
#> No date risk popup lat lon
#> <int> <date> <chr> <chr> <dbl> <dbl>
#> 1 227 2023-12-06 medium "2023年12月6日(水曜日)<br>足跡を… 43.0 141.
#> 2 226 2023-12-04 medium "2023年12月4日(月曜日)<br>足跡を… 42.9 141.
#> 3 225 2023-12-03 medium "2023年12月3日(日曜日)\n7時57分<… 43.0 141.
#> 4 224 2023-12-02 high "2023年12月2日(土曜日)\n18時45分… 42.9 141.
#> 5 223 2023-11-30 medium "2023年11月30日(木曜日)\n6時0分<… 43.0 141.
#> 6 222 2023-11-28 high "2023年11月28日(火曜日)\n21時55… 43.0 141.
過去５年分についても月別に報告数を整理したPDFファイルがアップロードされていたので、比較用に打ち込みました。12月から翌3月までは合計報告数のみなので省略しています。
higuma_past_data <-
tibble::tribble(~ year, ~ month, ~ count,
2022, 4, 18,
2022, 5, 17,
2022, 6, 15,
2022, 7, 24,
2022, 8, 22,
2022, 9, 14,
2022, 10, 30,
2022, 11, 17,
2021, 4, 5,
2021, 5, 21,
2021, 6, 30,
2021, 7, 46,
2021, 8, 31,
2021, 9, 19,
2021, 10, 15,
2021, 11, 8,
2020, 4, 5,
2020, 5, 12,
2020, 6, 13,
2020, 7, 31,
2020, 8, 13,
2020, 9, 6,
2020, 10, 7,
2020, 11, 4,
2019, 4, 21,
2019, 5, 35,
2019, 6, 42,
2019, 7, 30,
2019, 8, 47,
2019, 9, 7,
2019, 10, 11,
2019, 11, 1,
2018, 4, 5,
2018, 5, 21,
2018, 6, 35,
2018, 7, 37,
2018, 8, 17,
2018, 9, 5,
2018, 10, 12,
2018, 11, 5)
データ概観
データを概観します。今年度分を積み上げ棒グラフ、過去５年分を折れ線グラフとして表示しました。
higuma_data |>
dplyr::mutate(month = lubridate::month(date)) |>
dplyr::summarise(count = dplyr::n(), .by = c(risk, month)) |>
ggplot2::ggplot(ggplot2::aes(month, count, fill = risk)) +
ggplot2::theme_light(base_size = 20) +
ggplot2::geom_bar(stat = "identity") +
ggplot2::geom_line(data = higuma_past_data, inherit.aes = FALSE,
ggplot2::aes(month, count, col = as.character(year), group = year)) +
ggplot2::scale_fill_manual(values = c("grey20", "grey50", "grey80")) +
ggplot2::scale_x_continuous(breaks = 4:12, labels = paste0(4:12, "月")) +
ggplot2::labs(x = NULL, y = "報告数 [件]", col = "年", fill = "危険度")
今年６月は報告数が過去５年と比較して多く、１日あたり２回程度の報告があったようです。10・11月も比較的報告数が多めです。ただし７・８月の報告数はむしろ少なめで、通年で見ると2023年は極端にヒグマの目撃が多かった年だったというわけではないようです。今年の７・８月は道民には暑すぎる天候経過だったので、もしかしたら人間側のアウトドアな活動が少なく報告も減ったなどの影響もあるかもしれません。とはいえ、このデータを見る限りでは東北などの劇的な変化とは状況が異なる印象です。
地図化
leafletを使って地図を作成していきます。危険度ごとにマッピングするアイコンの色を変えてみました。Icon rainbowさんのアイコンを使わせていただきました。事前にqiitaにアップロードした画像をiconUrlで指定して読み込んでいます。
higuma_icon <-
list(extreme = leaflet::makeIcon(iconWidth = 24, iconHeight = 24, iconAnchorX = 12, iconAnchorY = 12,
iconUrl = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/106460/dc6a4d52-91a0-e4b6-ed5f-1d8bae858c98.png"),
high = leaflet::makeIcon(iconWidth = 24, iconHeight = 24, iconAnchorX = 12, iconAnchorY = 12,
iconUrl = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/106460/24995871-942e-e42f-ebc6-40f969af18c7.png"),
medium = leaflet::makeIcon(iconWidth = 24, iconHeight = 24, iconAnchorX = 12, iconAnchorY = 12,
iconUrl = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/106460/fdd9319e-a633-b63d-deb5-60fed4412dee.png"))
higuma_data |>
leaflet::leaflet() |>
leaflet::addProviderTiles(provider = providers$Esri.WorldImagery) |>
leaflet::addMarkers(data = dplyr::filter(higuma_data, risk == "extreme"),
lng = ~ lon, lat = ~ lat, popup = ~ popup, icon = higuma_icon$extreme) |>
leaflet::addMarkers(data = dplyr::filter(higuma_data, risk == "high"),
lng = ~ lon, lat = ~ lat, popup = ~ popup, icon = higuma_icon$high) |>
leaflet::addMarkers(data = dplyr::filter(higuma_data, risk == "medium"),
lng = ~ lon, lat = ~ lat, popup = ~ popup, icon = higuma_icon$medium) |>
leaflet::addMarkers(lat = 43.015, lng = 141.41) |> # 札幌ドーム
leaflet::addMarkers(lat = 43.0625, lng = 141.3535) # 札幌時計台
画像中心あたりのマーカーが札幌駅近くにある時計台、その右下のマーカーが札幌ドームで直線距離は7 kmほどです。札幌駅から一時間ほど歩いて山に入ればヒグマゾーンで、濃赤マーカーで示した子連れの危険な個体も都市部辺縁で散見されていることがよくわかります。
自分はインドア派ですが屋外系アクティビティでは改めてよく注意しようと思います。