いい感じの作図|北海道の河川集水域地図

こんにちは。博士後期課程3年の植村です。

 

今年はなんとブログの投稿がここまでゼロ!

こんな年はかつてあったのでしょうか・・・!!

そこで、それを阻止すべく、How to 系の投稿をします。

 

研究、特に論文書き(描き)のプロセスで欠かせないものの一つに作図があると思います。中でも、生態学であれば、調査地図や分布図がFig. 1にくることも多いです(Fig. 1に持ってくる論文が好きです)。

そんな地図ですが、凝ったものを作成しようとするとなかなか時間がかかります。数十年前までは手描き、一昔前まではデータ整理をExcel→作図をGISソフトウェア(例えば、QGIS)でカチカチとクリックしながら行った方も多いと思います。しかし、現在ではRを用いて、データ整理から作図までほとんど全てが完結します。

今回はRを用いて、北海道の水系ごとの集水域地図(+アルファ)を作成したいと思います。シンプルな作図なのですが、Googleで日本語検索してもこの類がほとんど見つかりませんでした。集水域は河川や河川に生息する生物にとって重要な環境要素なので、作図の”需要”があると勝手に思っています。

そんな中、数少ない日本語ソースとして参考にさせていただいたサイトを先に挙げます。

 

手順ですが、以下の3ステップで整理、作図をします。では、早速コードを示します(注:下記のコードの使用は自己責任でお願いいたします)。

# PC: Apple Macbook Pro 14 (Apple M3 Pro, macOS 15.1.1)
# R version 4.4.1 (2024-06-14)

# 必要なライブラリ(パッケージ)を読み込む
pacman::p_load(sf,
               dplyr,
               ggplot2,
               mapview,
               ggspatial,
               tidyverse)

# 1. 白地図を作成するためのデータ整理する-----------------------------------
# -------------------------------------------------------------------------
# データをダウンロードする
## Retrieve from https://www1.gsi.go.jp/geowww/globalmap-gsi/download/data/gm-japan/gm-jpn-all_u_2_2.zip
# -------------------------------------------------------------------------
# 行政区域データを読み込む(データフレームには”gm_jpn”と名づけた)
gm_jpn <- read_sf("gm-jpn-all_u_2_2/polbnda_jpn.shp") |>
  ## CRSを変換する (= WGS84 "地理座標系" = LngLat object)
  st_transform(crs = 4326) |>
  ## 行政区域レベルのデータを都道府県レベルのデータに集約する
  group_by(nam) |>   # nam = 都道府県名 
  summarize()

# データフレームを確認する
gm_jpn |> head()

# 北海道の”白”地図を描画する
ggplot() +
  geom_sf(data = gm_jpn[gm_jpn$nam == "Hokkai Do",]) +
  theme_void()




# 2. 集水域地図を作成するためのデータ整理する-------------------------------
# -------------------------------------------------------------------------
# データをダウンロードする
## Retrieve from https://nlftp.mlit.go.jp/ksj/gmlold/datalist/gmlold_KsjTmplt-W12.html
# -------------------------------------------------------------------------
# 北海道の流域界・非集水域データを読み込む(データフレームには”basin_hkd”と名づけた)
basin_hkd <- st_read("W12-52A_01_WatershedBoundary_hkd.shp") |>
  ## CRSを変換する (= WGS84 "地理座標系" = LngLat object)
  st_transform(crs = 4326)

# データフレームを確認する
basin_hkd |> head()

# データフレームの列名を日本語から英語に変換する(日本語のままでも操作可能)
colnames(basin_hkd)[colnames(basin_hkd) == "管轄地建"] <- "kankatsuchiken"
colnames(basin_hkd)[colnames(basin_hkd) == "水系NO"] <- "suikei_NO"
colnames(basin_hkd)[colnames(basin_hkd) == "河川CD"] <- "kasen_CD"
colnames(basin_hkd)[colnames(basin_hkd) == "河川名"] <- "kasenmei"
colnames(basin_hkd)[colnames(basin_hkd) == "水系名"] <- "suikeimei"
colnames(basin_hkd)[colnames(basin_hkd) == "水系CD"] <- "suikei_CD"
colnames(basin_hkd)[colnames(basin_hkd) == "都道府県CD"] <- "todofuken_CD"
colnames(basin_hkd)[colnames(basin_hkd) == "都道府県名"] <- "todofukenmei"
colnames(basin_hkd)[colnames(basin_hkd) == "都道府県NO"] <- "todofuken_NO"
colnames(basin_hkd)[colnames(basin_hkd) == "備考"] <- "bikou"
# -------------------------------------------------------------------------
# 河川名ごとに集水域面積(”area”と名付けた)を計算する
basin_hkd$area <- basin_hkd |>
  st_area() |>
  as.numeric()

# 具体例で確認する(河川名ごとに異なっていればOK)
unique(basin_hkd[basin_hkd$suikeimei == "斜里川",]$area)

# 水系名ごとに集水面積(”area_suikei”と名付けた)を計算する
basin_hkd <- basin_hkd |>
  group_by(suikeimei) |>
  mutate(area_suikei = sum(area)) |>
  ungroup()  # Remove grouping

# 具体例で確認する(河川名ごとに異なっていなければOK)
unique(basin_hkd[basin_hkd$suikeimei == "斜里川",]$area_suikei)



# 3. 北海道の水域名ごとに集水域を描画し集水面積で色付けした地図を作成する

# 動的なmapviewで河川名ごとの集水域を確認する(結構便利)
basin_hkd |> mapview()

Screenshot

# 北海道の河川名ごとに集水域(集水域面積で色付け)を描画する
ggplot() +
  geom_sf(data = basin_hkd,
          aes(fill = area)) +
  theme_void()


# 完成地図を描画する
ggplot() +
  ## 北海道の白地図を描画する
  geom_sf(data = gm_jpn[gm_jpn$nam == "Hokkai Do",],
          color = "#000000", fill = "#ffffff",
          linewidth = .15,
          show.legend = F) +
  ## 北海道の水域名ごとに集水域(集水域面積で色付け)を描画する
  geom_sf(data = basin_hkd, aes(fill = log(area_suikei),
                                color = log(area_suikei))) +
  ## 距離スケールバーを入れる
  annotation_scale(location = "tr", width_hint = .2) +
  scale_x_continuous(limits = c(139.6, 145.7), breaks = seq(141, 145, 2)) +
  scale_y_continuous(limits = c(41.5, 45.5), breaks = seq(42, 45, 1)) +
  scale_fill_distiller(type = "seq", palette = "Greys") +
  scale_color_distiller(type = "seq", palette = "Greys") +
  labs(title = NULL,
       x = "Longitude",
       y = "Latitude",
       fill = "log(basin area)",
       color = "log(basin area)"
       ) +
  theme_bw() +
  ## 凡例の位置を微調整する
  theme(legend.position = c(.85, .175),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.key.height = unit(.25, 'cm'),
        legend.key.width = unit(.5, 'cm'),
        axis.ticks.length = unit(-0.075, "cm"))



 

いかがでしたか?こんなシンプルなコードで”いい感じ”の作図ができてしまうR、素晴らしいですね。

Rなど作図や統計の手段としてのソフトウェアの発達は、研究の発展とも大きく関係していると思います。実際、最近の研究では、生態学や進化学の分野において必要なコードを共有している論文がより速く引用数を増やしているという結果も出ています(Maitner et al. 2024 Ecol. Evol.)。

このブログの読者層は修士課程の受験生〜釣り人まで幅広いですが、今回のポストが特にRの初学者の方々にとって”Rで作図”のきっかけとなりましたら幸いです。

 

 

では、D論に戻ります。

 

(植村)