【Google Colab×R】津波浸水想定域データを可視化する方法

はじめに

QGISとRstudio両方をつかって、可視化や分析をするのが大変…という悩みを解決するために、Google ColabでRコードを使って図示をする方法を紹介します。

参考書籍の「Rではじめる地理空間データの統計解析入門」ではRコードを使って地図上に可視化する方法がたくさん載っているため、ぜひ参考にしてみてください!

Google Colabを使うメリットデメリット

Google Colabを使うメリット

  • いろいろなデバイスで編集確認ができる
  • 1回閉じても実行しないで出力の確認ができる
  • RとPythonを併用できる
  • QGISなどと比べてエラー内容が分かりやすい

iPhoneやiPadで編集できる、電源が落ちてしまっても1から実行する必要がない、Rで難しいGoogleドライブとのマウントやGoogleドライブのファイル読み込みをPythonで対応できる、エラー内容を検索するとすぐに直し方が分かるなどのメリットがあります。

Google Colabを使うデメリット

  • Rパッケージのダウンロードを毎回しないといけない

Rstudioと同じようにパッケージをインストールすると毎回時間がかかってしまいます。

ステップ1国土数値情報ダウンロードサイトからシェープファイルをダウンロードする

国土数値情報ダウンロードサイトリンク↓

https://nlftp.mlit.go.jp

今回は津波浸水想定データをダウンロードしていきます。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A40-2024.html

下の方から岩手県を選択してダウンロードします。

ファイルはGoogleドライブのマイドライブの「+新規」から新しいフォルダを作成(今回は、A40-22_03_GMLという名前など)からダウンロードファイルを保存することで読み込みのときのエラーを減らすことができます。そのままZIPファイルの中身をGoogleドライブに移すと読み込めないことがあります。

これでデータをダウンロードして保存する下準備が終わります。

ステップ2 GoogleColabにアクセスする

GoogleColabにアクセスします。

ノートブックを新規作成します。

ファイルの移動から先ほどシェープファイルを保存したフォルダに移動します

ステップ3 実際にGoogleドライブをGoogleColabで読み込む

Googleドライブと繋げておきます。基本的には、ドライブと繋げる、読み込みなどはPythonが適しているため、右下のランタイムオプションは最初のままPythonにしておきます。

入力コードです。

# 1. マウント
from google.colab import drive
drive.mount('/content/drive')

# 2. 地図パッケージの土台となるシステム部品を入れておく(約30秒)
!apt-get -qq update
!apt-get -qq install libudunits2-dev libgdal-dev libgeos-dev libproj-dev

# 3. Rの呼び出し準備
%load_ext rpy2.ipython

以下がセルの出力です。

Mounted at /content/drive
W: Skipping acquire of configured file 'main/source/Sources' as repository 'https://r2u.stat.illinois.edu/ubuntu jammy InRelease' does not seem to provide it (sources.list entry misspelt?)
The rpy2.ipython extension is already loaded. To reload it, use:
  %reload_ext rpy2.ipython

ステップ4 ライブラリなどの読み込み

以下からはRコードを使うため、1行目に「%%R」と入力を忘れないようにします。Rstudioのようにライブラリをインストールしたところ20分程かかったため、以下のようにサーバーを指定するのがおすすめです。

入力コード

%%R
# 4. 高速な完成品サーバーを指定
options(repos = c(CRAN = "https://packagemanager.posit.co/cran/__linux__/jammy/latest"))

# 5. すべて一緒にインストール(これで相性問題が起きません!約1〜2分で終わります)
install.packages(c("sf", "NipponMap", "RColorBrewer", "showtext"))

# 6. パッケージの読み込み
library(sf)
library(NipponMap)
library(RColorBrewer)
library(showtext)

# 7. 日本語フォントの設定
font_add_google("Noto Sans JP", "notojp")
showtext_auto()

print("すべてのパッケージの準備が完了しました!")

出力コード

 [1] "すべてのパッケージの準備が完了しました!"
Installing packages into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
trying URL 'https://packagemanager.posit.co/cran/__linux__/jammy/latest/src/contrib/sf_1.1-1.tar.gz'
trying URL 'https://packagemanager.posit.co/cran/__linux__/jammy/latest/src/contrib/NipponMap_0.2.tar.gz'
trying URL 'https://packagemanager.posit.co/cran/__linux__/jammy/latest/src/contrib/RColorBrewer_1.1-3.tar.gz'
trying URL 'https://packagemanager.posit.co/cran/__linux__/jammy/latest/src/contrib/showtext_0.9-8.tar.gz'

The downloaded source packages are in
	‘/tmp/RtmpSbWtF0/downloaded_packages’
Loading required package: sysfonts
Loading required package: showtextdb

ステップ5 Googleドライブのファイル読み込み

Googleドライブのファイルを読み込んでいきます。これはPythonコードのため、「%%R」はいらないです。

入力コード

#作業フォルダへの移動
import os
os.chdir('/content/drive/MyDrive/A40-22_03_GML')

ステップ6 実際に過地図で可視化していく

色分けの確認

実際にどのような色分けがあるか確認していきます。

入力コード

%%R
display.brewer.all()

出力コード

見にくいですが、RdYlGnが浸水ダメージが大きそうな赤から、ダメージが小さそうな緑にグラデーションとなっているため、今回はRdYlGnを使います。

shp(シェープファイル)を読み込んで、sf形式としてRで操作

国土数値情報ダウンロードサイトのようなサイトでは.shpとしてダウンロードしますが、Rではsf形式で操作していきます。

入力コード

%%R
# sfパッケージの呼び出し
library(sf)

# Pythonで移動したフォルダ内にいるので、ファイル名だけで読み込めます
# (※ファイル名はダウンロードしたものに合わせてください。ここでは岩手県の例です)
tsunami_data <- read_sf("A40-22_03.shp")

# 座標系が最初からどう設定されているか確認する
print(st_crs(tsunami_data))

# 試しに形だけプロットしてみる
plot(st_geometry(tsunami_data), col="blue")

出力コード

Coordinate Reference System:
  User input: JGD2011 
  wkt:
GEOGCRS["JGD2011",
    DATUM["Japanese Geodetic Datum 2011",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Japan - onshore and offshore."],
        BBOX[17.09,122.38,46.05,157.65]],
    ID["EPSG",6668]]

岩手県の海外線が出てきます!

ファイルの列名の確認

列名が分からないとplotするときに列を選ぶことができないため確認します。

入力コード

%%R
tsunami_data

出力コード

Simple feature collection with 233714 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 141.5845 ymin: 38.93819 xmax: 142.0724 ymax: 40.45024
Geodetic CRS:  JGD2011
# A tibble: 233,714 × 4
   A40_001 A40_002 A40_003                                              geometry
   <chr>   <chr>   <chr>                                           <POLYGON [°]>
 1 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7035 38.9386, 141.7035 38.93855,…
 2 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7034 38.9386, 141.7034 38.93855,…
 3 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7023 38.93856, 141.7023 38.93851…
 4 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7019 38.93865, 141.7019 38.93861…
 5 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7035 38.9386, 141.7035 38.93864,…
 6 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7025 38.93874, 141.7025 38.93869…
 7 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7016 38.93879, 141.7016 38.93875…
 8 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7034 38.93882, 141.7034 38.93878…
 9 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.704 38.93886, 141.704 38.93882, …
10 岩手県  03      0.01m以上 ~ 0.3m未満 ((141.7035 38.93891, 141.7035 38.93887…
# ℹ 233,704 more rows
# ℹ Use `print(n = ...)` to see more rows

浸水カテゴリの確認

浸水データが入っている列がA40_003と分かったため、その列のデータのカテゴリ分けがどのようになっているか確認します。

入力コード

%%R
# A40_003列に含まれるパターンと、その件数を一覧表示する
table(tsunami_data$A40_003)

# (もし種類だけをシンプルに見たい場合はこちら)
# unique(tsunami_data$A40_003)

出力コード

 0.01m以上 ~ 0.3m未満   0.3m以上 ~ 0.5m未満   0.5m以上 ~ 1.0m未満 
                 25936                  23280                  37319 
  1.0m以上 ~ 3.0m未満 10.0m以上 ~ 20.0m未満              20.0m以上 
                 51930                  14870                   1783 
  3.0m以上 ~ 5.0m未満  5.0m以上 ~ 10.0m未満 
                 49406                  29190 

8つのカテゴリ分けがあることが分かります。

浸水カテゴリの順序設定と可視化

浸水の深さが順序順になるように設定して地図上に浸水域を確認していきます。

入力コード

%%R
library(RColorBrewer)

# 1. 正しい「深さの順番」を小さい順に定義します
depth_levels <- c(
  "0.01m以上 ~ 0.3m未満",
  "0.3m以上 ~ 0.5m未満",
  "0.5m以上 ~ 1.0m未満",
  "1.0m以上 ~ 3.0m未満",
  "3.0m以上 ~ 5.0m未満",
  "5.0m以上 ~ 10.0m未満",
  "10.0m以上 ~ 20.0m未満",
  "20.0m以上"
)

# 2. A40_003列を「順序付きファクター型」に変換します
# (これにより、見た目は文字のまま、裏側で上記の正しい順番がセットされます)
tsunami_data$A40_003 <- factor(tsunami_data$A40_003, levels = depth_levels, ordered = TRUE)

# 3. 8段階のグラデーションカラーを作成
pal <- brewer.pal(8, "Blues")

# 4. 地図を描画
plot(tsunami_data["A40_003"], 
     pal = pal, 
     main = "津波浸水想定域", 
     key.pos = 4,         # 凡例を右側に配置
     key.width = lcm(4))  # 凡例の幅を広げて文字を見やすくする

出力コード

少し海外線が青いですが、ほぼ海岸線データで見えなくなっています。

海岸線を消す

海岸線で見にくいため海岸線を消していきます

入力コード

%%R
# border = NA を追加して、黒い枠線を消します
plot(tsunami_data["A40_003"], 
     pal = pal, 
     border = NA,         # ← ここを追加!
     main = "津波浸水想定域", 
     key.pos = 4, 
     key.width = lcm(4))

出力コード

カテゴリ分けされていることは分かりましたが、画質の問題や、全部青のため、見にくいです。

色の変更

最後に見やすい色にしてPDFに保存していきます。

入力コード

%%R
library(RColorBrewer)

# 1. 色を「RdYlGn」にし、浅い(緑) → 深い(赤)になるよう反転させます
pal <- rev(brewer.pal(8, "RdYlGn"))

# 2. 高画質なPDFファイルとして保存を開始(widthとheightで用紙サイズを大きく設定)
# ※ファイル名は自由に変更可能です
pdf("iwate_tsunami_highres.pdf", width = 10, height = 15)

# 3. 描画処理(画面には出ず、PDFに書き込まれます)
plot(pref3["A40_003"], 
     pal = pal, 
     border = NA,         # 枠線は「なし」が一番綺麗です
     main = "岩手県 津波浸水想定域", 
     key.pos = 4, 
     key.width = lcm(4),
     key.length = 0.8)    # 凡例の長さもバランスよく調整

# 4. 保存を完了する(※この行を忘れないでください!)
dev.off()

出力コード

拡大すると少し見えるようになりました。(岩手県が大きいことがよくわかります)

まとめ

今回はGoogleColabを使ってRで地図に浸水想定データを図示しました。ぜひご自身の気になる県や、気になるテーマを図示してみてください。

コメント

タイトルとURLをコピーしました