小松市のクマ出没予測マップの作成

「クマ共生ハッカソン」の一環として、石川県小松市のクマ出没予測マップを作成しました。クマそのものの生息の有無ではなく、出没情報の有無をモデル化して、出没確率を予測したものです。RのsdmTMBパッケージを使用して、時空間自己相関を組み込んだ種分布モデル（Species Distribution Model: SDM）を作成し、予測マップを作成しました。

前回、QGISで作成した土地利用とクマ出没情報の位置を重ねた地図から、小松市近辺を拡大したものを貼っておきます。

データ

クマ出没情報は、いしかわオープンデータカタログの「石川県クマ出没データ」からダウンロードしたものに前処理を加えたものを使用します。

土地利用データとして、国土数値情報の土地利用3次メッシュデータを利用します。石川県にかさなる部分をダウンロードしておきます。

また、国土数値情報の行政区域データから石川県のデータをダウンロードしておきます。

以下のRコードでは、tidyverse, sf, nanoparquet, sdmTMBの各パッケージを読み込み、データファイルの場所を指定しています。

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(nanoparquet)
library(sdmTMB)

# データディレクトリ
data_dir <- "data"

# クマ出没データ
data_file <- file.path(data_dir, "kuma_data.parquet")

# 土地利用データ
landuse_files <- 
  c("L03-a-21_5436.geojson",
    "L03-a-21_5536.geojson",
    "L03-a-21_5537.geojson",
    "L03-a-21_5636.geojson",
    "L03-a-21_5637.geojson") |>
  purrr::map_chr(~ file.path(data_dir, .x))

# 石川県行政区域データ
ishikawa_pref_file <- file.path(data_dir, "N03-20250101_17.geojson")

クマ出没データ

Parquet形式で保存しておいたクマ出没データを読み込みます。すべての出現タイプを組み合わせたものになります。土地利用データにあわせて、500mメッシュコードから3次（1km）メッシュコードを作成しておきます。

kuma_data <- read_parquet(data_file) |>
  # 3次メッシュコード
  # 500mメッシュコードの先頭から8文字
  dplyr::mutate(mesh_code = str_sub(`500mメッシュコード`, 1, 8))

土地利用データ

土地利用データを読み込みます。データ中の各土地利用面積はm²単位ですが、km²単位に変換します。また、3次メッシュコードを文字型に変換しておきます。

landuse <- landuse_files %>%
  purrr::map(\(f) {
    sf::st_read(f, quiet = TRUE)
  }) |>
  dplyr::bind_rows() |>
  # m^2単位をkm^2単位にする
  dplyr::mutate_at(2:14, `/`, 10^6) |>
  # 3次メッシュコードを文字型にする
  dplyr::mutate(mesh_code = as.character(`メッシュ`))

行政区域データ

国土数値情報の石川県の行政区域データを読み込みます。

ishikawa <- sf::st_read(ishikawa_pref_file, quiet = TRUE)

小松市の行政区域を抽出します。

komatsu <- ishikawa |>
  dplyr::filter(N03_004 == "小松市") |>
  sf::st_union()

土地利用3次メッシュのうち小松市の行政区域と重なるものを抽出します。

int <- sf::st_intersects(landuse, komatsu, sparse = FALSE)
landuse_komatsu <- landuse[int, ]

クマ出没情報の地図

小松市のクマ出没データについて、3次メッシュごと年ごとの出没情報のありなしを集計して、p_sumに格納します。

p_sum <- kuma_data |>
  dplyr::filter(`市町名` == "小松市") |>
  dplyr::rename(year = `出没年`) |>
  dplyr::group_by(mesh_code, year) |>
  dplyr::summarise(present = as.numeric(n() > 0),
                   .groups = "drop")

6年間（2019〜2024年）で、3次メッシュごとにクマの出没情報があった年の割合を地図でしめします。

geometry <- landuse_komatsu |>
  dplyr::select(mesh_code, geometry)

prop <- p_sum |>
  dplyr::group_by(mesh_code) |>
  dplyr::summarise(prop = n() / 6, .groups = "drop")

geometry |>
  dplyr::left_join(prop, by = "mesh_code") |>
  dplyr::mutate(prop = replace_na(prop, 0)) |>
  ggplot(aes(fill = prop)) +
  geom_sf() +
  scale_fill_viridis_c(name = "割合") +
  theme_minimal(base_family = "Noto Sans JP")

図 2: 小松市の３次メッシュごとの、クマの出没情報があった年の割合（2019〜2024年）

統計モデル

全出没タイプを統合したデータについて、出没情報のありなしをモデル化します。

説明変数として用いる可能性のある環境データとして以下を用意しました。森林からの出没は、森林と農地・市外との境界付近で多発することが分かっていますので、森林面積の2次関数をあてはめることを想定して2乗の項もつくっています。同様に、河川からの出没に対応することを想定して、河川地及び湖沼面積も用意しています。とはいえ、これらすべてを使用するとは限りません。

forest: 3次メッシュ中の森林面積(km²)
forest2: 3次メッシュ中の森林面積(km²)の2乗
water: 3次メッシュ中の河川地及び湖沼面積(km²)
water2: 3次メッシュ中の河川地及び湖沼面積(km²)の2乗

env <- landuse_komatsu |>
  dplyr::mutate(forest = `森林`,
                forest2 = `森林`^2,
                water = `河川地及び湖沼`,
                water2 = `河川地及び湖沼`^2) |>
  as.data.frame() |>
  dplyr::select(mesh_code, forest, forest2, water, water2)

3次メッシュコードと年のすべての組み合わせをつくっておき、これに出没のありなしデータと環境データを結合させます。

同時に、ブナ豊凶をしめすダミー変数buna_poorを加えます。大凶作の年（2020年）は1、それ以外の年は0をとります。

p_year <- tidyr::expand_grid(
  mesh_code = unique(landuse_komatsu$mesh_code),
  year = 2019:2024
) |>
  dplyr::left_join(p_sum, by = c("mesh_code", "year")) |>
  dplyr::mutate(present = replace_na(present, 0)) |>
  dplyr::left_join(env, by = "mesh_code") |>
  dplyr::mutate(buna_poor = if_else(year == 2020, 1, 0))

sdmTMBで使用するメッシュ作成のため、UTM座標系の座標を求めます。

まず、landuse_komatsuの座標系をWGS84 (EPSG: 4326)に変換します。各メッシュの重心の座標を求め、これからさらに対応するUTM座標系(53N)を求めています。

coord <- landuse_komatsu |>
  sf::st_transform(4326) |>
  sf::st_make_valid() |>
  sf::st_centroid() |>
  sf::st_coordinates() |>
  as.data.frame() |>
  dplyr::rename(longitude = X, latitude = Y) |>
  add_utm_columns(ll_names = c("longitude", "latitude"))

出没あり・なしデータに経度・緯度データを結合します。

coord2 <- coord |>
  bind_cols(
    landuse_komatsu |>
      dplyr::select(mesh_code)
    )
p_year_coord <- p_year |>
  dplyr::left_join(coord2, by = "mesh_code")

メッシュの作成

sdmTMBで使用するメッシュを作成します。辺の最小の長さ(cutoff)は2kmとしました。この値は大きくするほど、ガウスランダム場がよりなめらかになります。

mesh <- make_mesh(p_year_coord, xy_cols = c("X", "Y"), cutoff = 2)
plot(mesh)

モデルの作成とあてはめ

2種類のモデルを用意しました。

森林面積(F)とその2乗、河川・湖沼面積(W)、ブナの豊凶(B, 大凶作: 1, それ以外0)を説明変数とするモデル
森林面積(F)とその2乗、ブナの豊凶を説明変数とするモデル

ともに、これら説明変数のほかに、時空間自己相関をモデルに組み込んでいます。

目的変数Yはクマ出没情報のあり(Y=1)、なし(Y=0)とします。Yは、確率pのベルヌーイ分布にしたがうとします。

\[ Y \sim \mathrm{Bernoulli}(p) \]

確率pは、そのロジットlogit(p)が、線形予測子（以下の式はモデル1の線形予測子）と、空間および時間の変量効果（\(\epsilon\)）の和に等しいとします。

\[ \mathrm{logit}(p) = \beta + \beta_\mathrm{F} F + \beta_\mathrm{F2}F^2 + \beta_\mathrm{W}W + \beta_\mathrm{B}B + \epsilon \]

sdmTMB関数を使用して、それぞれのモデルをあてはめます。空間の変量効果はガウスランダム場、時間の変量効果は1次の自己回帰にしたがうとしています。

これらの環境要因以外の要因（誘引物など）は、ガウスランダム場としてメッシュに割り当てることを想定しています。

以下のRコードのモデル式では、forestは森林面積（モデル式ではF）、waterは河川・湖沼面積（モデル式ではW）、buna_poor（モデル式ではB）はブナ大凶作の年かどうかを示すダミー変数（1: 大凶作、0: 大凶作以外）です。

fit1 <- sdmTMB(present ~ forest + forest2 + water + buna_poor,
              data = p_year_coord,
              mesh = mesh,
              family = binomial(link = "logit"),
              spatial = "on",
              time = "year",
              time_varying = ~ 1,
              time_varying_type = "ar1",
              extra_time = 2025)
fit2 <- sdmTMB(present ~ forest + forest2 + buna_poor,
              data = p_year_coord,
              mesh = mesh,
              family = binomial(link = "logit"),
              spatial = "on",
              time = "year",
              time_varying = ~ 1,
              time_varying_type = "ar1",
              extra_time = 2025)

結果

あてはめ結果をチェックします。モデル1です。

sanity(fit1)
## ✔ Non-linear minimizer suggests successful convergence
## ✔ Hessian matrix is positive definite
## ✔ No extreme or very small eigenvalues detected
## ✔ No gradients with respect to fixed effects are >= 0.001
## ✔ No fixed-effect standard errors are NA
## ✔ No standard errors look unreasonably large
## ✔ No sigma parameters are < 0.01
## ✔ No sigma parameters are > 100
## ✔ Range parameter doesn't look unreasonably large

モデル2です。

sanity(fit2)
## ✔ Non-linear minimizer suggests successful convergence
## ✔ Hessian matrix is positive definite
## ✔ No extreme or very small eigenvalues detected
## ✔ No gradients with respect to fixed effects are >= 0.001
## ✔ No fixed-effect standard errors are NA
## ✔ No standard errors look unreasonably large
## ✔ No sigma parameters are < 0.01
## ✔ No sigma parameters are > 100
## ✔ Range parameter doesn't look unreasonably large

どちらも問題は見つかりませんでした。

AICを比較します。

c(AIC(fit1), AIC(fit2))
## [1] 861.7945 860.8953

AICの値はモデル2の方が多少小さく、またモデル2の方がより単純なモデルですので、より予測能力が高いモデルとしてモデル2を採用します。

なお、河川・湖沼面積の2乗も含めたモデルもいちおうあてはめてみましたが、AICは小さくなりませんでした。

2つのモデルのあてはめの結果の要約を表示します。

summary(fit1)
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: present ~ forest + forest2 + water + buna_poor
## Mesh: mesh (isotropic covariance)
## Time column: year
## Data: p_year_coord
## Family: binomial(link = 'logit')
##  
## Conditional model:
##             coef.est coef.se
## (Intercept)    -6.84    1.18
## forest         14.02    2.19
## forest2       -13.50    1.74
## water           2.07    1.94
## buna_poor       2.59    1.17
## 
## Time-varying parameters:
##                    coef.est   coef.se
## (Intercept)-2019  0.4800000 0.5100000
## (Intercept)-2020 -0.3300000 0.8400000
## (Intercept)-2021  0.3600000 0.4600000
## (Intercept)-2022 -0.2600000 0.4800000
## (Intercept)-2023 -0.6900000 0.4000000
## (Intercept)-2024  1.1400000 0.5000000
## (Intercept)-2025 -0.5500000 1.2600000
## rho-(Intercept)  -0.4827451 0.8322732
## 
## Matérn range: 7.59
## Spatial SD: 2.07
## Spatiotemporal IID SD: 0.54
## ML criterion at convergence: 420.897
## 
## See ?tidy.sdmTMB to extract these values as a data frame.

summary(fit2)
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: present ~ forest + forest2 + buna_poor
## Mesh: mesh (isotropic covariance)
## Time column: year
## Data: p_year_coord
## Family: binomial(link = 'logit')
##  
## Conditional model:
##             coef.est coef.se
## (Intercept)    -6.57    1.13
## forest         13.71    2.16
## forest2       -13.42    1.73
## buna_poor       2.59    1.16
## 
## Time-varying parameters:
##                    coef.est   coef.se
## (Intercept)-2019  0.4700000 0.5100000
## (Intercept)-2020 -0.3200000 0.8400000
## (Intercept)-2021  0.3600000 0.4600000
## (Intercept)-2022 -0.2600000 0.4800000
## (Intercept)-2023 -0.6900000 0.3900000
## (Intercept)-2024  1.1400000 0.5000000
## (Intercept)-2025 -0.5500000 1.2500000
## rho-(Intercept)  -0.4837248 0.8280261
## 
## Matérn range: 7.60
## Spatial SD: 2.03
## Spatiotemporal IID SD: 0.54
## ML criterion at convergence: 421.448
## 
## See ?tidy.sdmTMB to extract these values as a data frame.

モデル2のbuna_poor（モデル式ではB）の係数の最尤推定値は2.59でした。この値は、たとえばブナが大凶作でないときの出現確率が0.1だったとした場合、ブナが大凶作になると出現確率が0.6にまで上昇するという値になります。

予測

モデル2を用いて、2025年の3次メッシュごとのクマ出没確率を予測します。

ブナ大凶作の場合

ブナ大凶作の場合の予測値を求めます。なお、今年2025年は大凶作の予報です。

newdata <- env |>
  dplyr::left_join(coord2, by = "mesh_code") |>
  dplyr::mutate(year = 2025, buna_poor = 1)
pred_buna_poor <- predict(fit2, newdata, type = "response")

地図で表示します。

pred_buna_poor |>
  sf::st_as_sf() |>
  ggplot(aes(fill = est)) +
  geom_sf() +
  scale_fill_viridis_c(name = "出没確率", limits = c(0, 1)) +
  theme_minimal(base_family = "Noto Sans JP")

森林と市街地・農地との境界付近に1に近い出没確率が予測されたメッシュがかなりあります。

なお、これは人里近くにクマが出没する確率です。南側の森林の広がる地域は、出没情報があがってこないだけで、クマが生息している地域です。出没確率の予測値が低くてもクマはいますので、立ち入る際はじゅうぶんな注意が必要です。

ブナ大凶作でない場合

ブナ大凶作でない場合の予測値を求めます。

newdata <- env |>
  dplyr::left_join(coord2, by = "mesh_code") |>
  dplyr::mutate(year = 2025, buna_poor = 0)
pred_buna <- predict(fit2, newdata, type = "response")
summary(pred_buna$est)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001443 0.0006987 0.0024193 0.0338008 0.0162508 0.6176153

地図で表示します。

pred_buna |>
  sf::st_as_sf() |>
  ggplot(aes(fill = est)) +
  geom_sf() +
  scale_fill_viridis_c(name = "出没確率", limits = c(0, 1)) +
  theme_minimal(base_family = "Noto Sans JP")

ブナ大凶作でない場合には、ブナ大凶作の場合と比較して、出没確率は低く予測されました。

実際の観測データとの比較

2025年の出没メッシュを集計します。

pre_2025 <- file.path(data_dir, "bear_sightings_r7.csv") |>
  readr::read_csv() |>
  dplyr::filter(市町名 == "小松市") |>
  dplyr::mutate(mesh_code = jpmesh::coords_to_mesh(`経度`, `緯度`, 1) |>
                  as.character()) |>
  dplyr::group_by(mesh_code) |>
  dplyr::summarise(n = n(), .groups = "drop")
## Rows: 103 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): WKT, 出没日, 時間, 市町名, 場所, 備考, 頭数, 月
## dbl (2): 緯度, 経度
## lgl (1): 番号
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

予測結果のデータフレームと結合します。予測は、ブナ大凶作としておこなっています。

bind_2025 <- pred_buna_poor |>
  dplyr::left_join(pre_2025, by = "mesh_code") |>
  dplyr::mutate(pre_2025 = if_else(is.na(n), 0, 1))

横軸を出没確率、縦軸を出没の有無としてプロットします。

ggplot(bind_2025, aes(x = est, y = pre_2025)) +
  geom_point(color = "red", size = 3) +
  scale_x_continuous(breaks = seq(0, 1, 0.2)) +
  scale_y_continuous(breaks = c(0, 1), minor_breaks = NULL) +
  labs(x = "2025年の出没予測確率", y = "実際の出没のあり(1)なし(0)") +
  theme_bw(base_family = "Noto Sans JP")

図 6: 2025年ブナ大凶作の場合のクマ出没確率の予測値と、実際の出没状況の比較

予測確率が0.36程度のメッシュですでに出没が見られています。予測マップで確率が小さくても油断しないことが必要です。

おわりに

繰り返しますが、予測マップで出没確率が低いと予測されていても、クマに対する注意は必要です。森林が広がる場所は、出没情報は少なくともクマの生息域ですから、立ち入るときは音を立てて人の存在を知らせるといった対策が必要になります。クマの生息しない市街地でも、生ゴミを放置しないといった基本的な対策は日ごろからとっておくべきでしょう。

また、モデルの予測精度についても、今後のデータでさらに検証が必要と考えられます。