13: ラスターデータ

更新日

2025年7月17日

はじめに

第4回の講義で、空間データにはベクターデータとラスターデータの2種類があるという話をしました。しかし、これまでの演習では、ベクターデータだけを扱い、ラスターデータに触れる機会がありませんでした。そこで今回の演習では、Rでラスターデータを分析する例をいくつか紹介します。

演習

準備

パッケージのロード

# パッケージのロード
library(tidyverse)
library(readxl)
library(sf)
library(terra)
library(units)
library(ggrepel)

基盤地図情報数値標高モデル

データ

国土地理院が、航空測量などによって作成した数値標高モデル（DEM：Digital Elevation Mode）を公開しています。数値標高モデルは、標高のメッシュデータで、1mメッシュと5mメッシュ（航空レーザ測量）および10mメッシュ（地形図の等高線）の3種類のデータがあります。

数値標高モデルは基盤地図情報ダウンロードサービスからダウンロードできますが、ダウンロードするには利用者登録が必要です。

また、ダウンロードされるファイルはJPGIS（GML）形式で、そのままでは一般のGISソフトウェアで利用することができませんので、以下に示す変換ソフトなどを利用して、GeoTIFFなどのファイル形式に変換する必要があります。

種類	名称
Windowsアプリ	基盤地図情報標高DEMデータ変換ツール（株式会社エコリス）
MacOSアプリ	基盤地図標高変換（DemConv）（片柳由明さん）
QGISプラグイン	QuickDEM4JP（MIERUNE Inc. ）

今回は、GeoTIFFに変換済みのファイルを用意しましたので、以下のリンクからダウンロードしてください（学内専用）。

elevation.tiff ※数値標高モデル（国土地理院）の10mメッシュ（標高）データを加工して作成）

ダウンロードしたファイルをプロジェクトのデータフォルダに移動してください。

次に、国土数値情報から、全国の国・都道府県の機関データ（2022年）をダウンロードしてください。

国・都道府県の機関データ

ダウンロードしたP28-22.zipを展開し、P28-22フォルダをdataフォルダに移動してください。

標高データの地図表示

ラスターデータを読み込むには、terra::rast関数を使います。

# 標高データを読み込む
elevation <- rast("data/elevation.tiff")

標高データを描画します。

# 標高データの描画
plot(elevation, range = c(0,1053), col = terrain.colors(100))

※数値標高モデル（国土地理院）の10mメッシュ（標高）データを加工して作成）

佐賀市から福岡市にかけての範囲の標高データですね。

経路上の標高データの抽出

ここでは、ラスタの抽出を行なってみましょう。福岡県庁と佐賀県庁を結ぶ経路上の標高を抽出します。具体的には、経路を多くの点に分割し、その点の標高を抽出します。

まず、国土数値情報の国・都道府県の機関データから、福岡県庁と佐賀県庁のポイントデータを抽出します。

# 県庁ポイントデータの読み込み
fs <- 
  read_sf("data/P28-22/P28-22.geojson") |>
  filter(P28_003 %in% c("福岡県庁", "佐賀県庁")) |> 
  st_transform(crs(elevation))

osrm::osrmRoute関数を使って、福岡県庁と佐賀県庁の間の最短経路を探索します。

# 県庁間の最短経路探索
route <- osrm::osrmRoute(loc = fs)

最短経路を地図に重ねてみましょう。

# 県庁間の最短経路（地図）
plot(elevation, range = c(0,1053), col = terrain.colors(100))
plot(vect(route), add = TRUE, col = 'red')

terra::extract関数を使って、福岡県庁と佐賀県庁の標高データを抽出します。

# 県庁の標高
terra::extract(elevation, fs)

  ID elevation
1  1       4.1
2  2       4.5

福岡県庁と佐賀県庁いずれも、その標高はおよそ4 mであることが分かります。

次に、福岡県庁と佐賀県庁を結ぶ経路上の標高を調べましょう。まず、福岡県庁と佐賀県庁の最短経路上に、sf::st_segmentize関数を使って、線分上に点を追加します。さらにsf::st_cast関数を使って、POINTデータに変換します。

# 経路のセグメント化
route_seg <- route |> 
  st_segmentize(dfMaxLength = 50) |> 
  st_cast('POINT')

Warning in st_cast.sf(st_segmentize(route, dfMaxLength = 50), "POINT"):
repeating attributes for all sub-geometries for which they may not be constant

glimpse(route_seg)

Rows: 1,537
Columns: 5
$ src      <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "…
$ dst      <chr> "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "…
$ duration <dbl> 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, 4…
$ distance <dbl> 56.9749, 56.9749, 56.9749, 56.9749, 56.9749, 56.9749, 56.9749…
$ geometry <POINT [°]> POINT (130.4177 33.60706), POINT (130.4175 33.60686), P…

それぞれのポイントについて、福岡県庁からの距離を計算します。そのためにまず、ポイントから次のポイントまでの距離をsf::st_distance関数で計算しています。そして、cumsum関数を使って、その累積距離を計算しています。

# 福岡県庁からの距離
route_seg <- route_seg |> 
  mutate(
    geometry2 = geometry[row_number() + 1],
    length = st_distance(geometry, geometry2, by_element = TRUE),
    distance = cumsum(length))
glimpse(route_seg)

Rows: 1,537
Columns: 7
$ src       <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", …
$ dst       <chr> "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", …
$ duration  <dbl> 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, …
$ distance  [m] 26.57625 [m], 53.15250 [m], 79.72875 [m], 106.30500 [m], 132.8…
$ geometry  <POINT [°]> POINT (130.4177 33.60706), POINT (130.4175 33.60686), …
$ geometry2 <POINT [°]> POINT (130.4175 33.60686), POINT (130.4173 33.60667), …
$ length    [m] 26.57625 [m], 26.57625 [m], 26.57625 [m], 26.57625 [m], 26.576…

terra::extract関数を使って、経路上のポイントの標高を抽出します。

# 経路上の標高を抽出
route_elev <- terra::extract(elevation, route_seg)
route_seg <- route_seg |> 
  mutate(elevation = route_elev$elevation) |> 
  st_transform(6670)
glimpse(route_seg)

Rows: 1,537
Columns: 8
$ src       <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", …
$ dst       <chr> "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", …
$ duration  <dbl> 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, 47.83667, …
$ distance  [m] 26.57625 [m], 53.15250 [m], 79.72875 [m], 106.30500 [m], 132.8…
$ geometry  <POINT [m]> POINT (-54040.43 67474.16), POINT (-54055.65 67452.41)…
$ geometry2 <POINT [°]> POINT (130.4175 33.60686), POINT (130.4173 33.60667), …
$ length    [m] 26.57625 [m], 26.57625 [m], 26.57625 [m], 26.57625 [m], 26.576…
$ elevation <dbl> 3.9, 4.0, 4.0, 4.0, 4.0, 4.0, 3.9, 3.9, 3.8, 3.8, 3.8, 3.8, …

結果を地図に表示します。

# 標高地図の描画
ggplot() + 
  geom_sf(aes(color = elevation), data = route_seg) +
  scale_color_gradientn(colors = terrain.colors(100)) +
  theme_bw()

横軸に福岡県庁からの距離を、縦軸に標高をとったグラフを描いてみましょう。ここでは、ggplot2::geom_area関数でエリアプロットにしています。結構高低差がありますね。

# 標高グラフの描画
route_seg |> 
  mutate(
    distance = set_units(distance, 'km'),
    elevation = set_units(elevation, 'm')) |> 
ggplot(aes(x = distance, y = elevation)) + 
  geom_area(fill = gray(0.5)) +
  theme_minimal()

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_align()`).

経路上の地点に対応する標高を抽出することは、例えばハイキングの計画を立てることを想像してみると、役に立つことがイメージできるのではないでしょうか。ルートのどこに長い上り坂があるのか、休憩地点を計画したり、所要時間を見積もるのに使えると思います。また、道路建設などの土木計画にも適用可能でしょう。

高解像度土地利用土地被覆図

データ

宇宙航空研究開発機構（JAXA）地球観測研究センター（EORC）のALOS/ALOS-2解析研究プロジェクトおよび「課題分野型研究：生態系研究グループ」で作成した、高解像度土地利用土地被覆図が公開されています）。

これは、だいち2号（ALOS-2）が観測したデータなど基に作成した日本の土地利用、土地被覆のデータです。

JAXAのALOS利用推進プロジェクトのウェブサイトからデータをダウンロードすることが可能です。利用にはユーザ登録が必要ですが、利用条件等への同意のもと、だれでも無償で利用できます。

今回は、JAXAが提供している「日本域10m解像度土地利用土地被覆図【2024年】（バージョン25.04）」について、佐賀県周辺のデータを切り出したデータを用意したので、これをダウンロードして使用してください（学内専用）。

lulc_saga.tiff ※提供：高解像度土地利用土地被覆図（JAXA）

また、佐賀県の行政界データも使用します。国土数値情報から佐賀県の行政界データ（2025年）をダウンロードしてください。

行政区域データ

ダウンロードフォルダでN03-20250101_41_GML.zipを展開し、 N03-20250101_41_GMLフォルダをプロジェクトのデータフォルダに移動してください。

土地利用土地被覆図

ダウンロードしたデータを読み込んで表示してみましょう。

土地利用土地被覆データを、terra::rast関数を使って読み込みます。

# 土地利用土地被覆データ
lulc <- rast('data/lulc_saga.tiff')

佐賀県の行政区域データを読み込みます。

# 佐賀県の行政区域データ
saga <- read_sf('data/N03-20240101_41_GML/N03-20240101_41.geojson') |>
  rename(city = N03_004) |> group_by(city) |> summarise() |> 
  st_transform(crs(lulc))

土地利用土地被覆図に行政区域を重ねて描画します。

# 土地利用土地被覆図の描画
plot(lulc)
plot(vect(saga), add = TRUE)

今回利用するデータは、土地利用土地被覆が以下の15のカテゴリに分類されていて、カテゴリによって塗り分けられた地図が描画されます。

日本域高解像度土地利用土地被覆図（2024年）のデータ分類
no	categories	カテゴリ	色
#1	Water bodies	水域	■
#2	Built-up	人工構造物	■
#3	Paddy field	水田	■
#4	Cropland	畑地	■
#5	Grassland	草地	■
#6	DBF (deciduous broad-leaf forest)	落葉広葉樹	■
#7	DNF (deciduous needle-leaf forest)	落葉針葉樹	■
#8	EBF (evergreen broad-leaf forest)	常緑広葉樹	■
#9	ENF (evergreen needle-leaf forest)	常緑針葉樹	■
#10	Bare	裸地	■
#11	Bamboo forest	竹林	■
#12	Solar panel	ソーラーパネル	■
#13	Wetland	湿地	■
#14	Greenhouse	農業用温室	■
#15	Rock reef and Tidal flat	岩礁・干潟	■

ソーラーパネル比率

ソーラーパネルだけを取り出して地図にしてみます。佐賀県内にはそれほど大きな太陽光発電所はないようです（大牟田のメガソーラーが目立ちますね）。

# ソーラーパネル地図
solar <- lulc
values(solar) <- values(lulc) == 12
plot(solar)
plot(vect(saga), add = TRUE, border = "white")

では、佐賀県市町別のソーラーパネル面積比率を計算してみましょう。

行政区域ポリゴンを使って、市町ごとにカテゴリ数をカウントします。このような目的には、terra::freq関数が使えます。

# 市町別カテゴリ数のカウント
count <- freq(lulc, zones = vect(saga))
count <- count |> 
  mutate(name = str_c('cat', str_pad(value, 2, 'left', '0'))) |> 
  select(-layer, -value) |>  
  pivot_wider(values_from = count)
glimpse(count)

Rows: 20
Columns: 16
$ zone  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
$ cat01 <dbl> 19846, 425, 29598, 146271, 5367, 61459, 406, 3315, 1630, 2747, 8…
$ cat02 <dbl> 161535, 55909, 309096, 1075067, 100511, 535453, 83296, 127769, 3…
$ cat03 <dbl> 99214, 22636, 237713, 604037, 42272, 348937, 17264, 79404, 11570…
$ cat04 <dbl> 184936, 45585, 201303, 946243, 80403, 383952, 19195, 98179, 3042…
$ cat05 <dbl> 58462, 9514, 221688, 302478, 33439, 399256, 20588, 102126, 12227…
$ cat06 <dbl> 8635, 3510, 84808, 206293, 33305, 265933, 7547, 68978, 2865, 443…
$ cat07 <dbl> 80, 24, 255, 1809, 231, 937, 51, 107, 1, 58, 153, 163, 62, 84, 2…
$ cat08 <dbl> 56752, 27225, 952142, 612091, 111854, 1651044, 43680, 389308, 28…
$ cat09 <dbl> 37844, 24791, 1068850, 1701869, 147008, 2012160, 84043, 263437, …
$ cat10 <dbl> 7496, 2084, 41324, 60767, 4697, 97977, 2130, 10744, 3948, 9641, …
$ cat11 <dbl> 39310, 12395, 332850, 230223, 47895, 835812, 28527, 177036, 1729…
$ cat12 <dbl> 1890, 969, 4747, 5688, 630, 11689, 588, 3276, 697, 1199, 3700, 1…
$ cat13 <dbl> 9801, 1289, 60759, 72966, 3222, 132726, 693, 20354, 4047, 6737, …
$ cat14 <dbl> 4993, 1297, 10767, 42309, 2396, 64642, 1153, 1653, 736, 2755, 25…
$ cat15 <dbl> 3473, 466, 2594, 22368, 438, 4343, NA, 93, 676, 335, 62, 3906, N…

カウントに占めるソーラーパネルの比率を計算します。 dplyr::select関数で、cat01（水域）を取り除き、それ以外のcatで始まる列（cat02〜cat15）の合計をtotal列として追加しています。そして、cat12（ソーラーパネル）がtotalに占める比率（％）をsolar列として追加しています。

solar_ratio <- 
  bind_cols(saga, count) |> select(-cat01) |> 
  mutate(
    total = rowSums(across(starts_with("cat")), na.rm = TRUE),
    solar = cat12 / total * 100
  ) |> 
  st_transform(6670)

市町ごとのソーラーパネル比率を塗り分け地図にしてみます。最もソーラーパネル比率が高かったのは、上峰町の0.47％でした。

ggplot() + 
  geom_sf(aes(fill = solar), data = solar_ratio) +
  scale_fill_distiller(direction = 1) +
  geom_text_repel(aes(geometry = geometry, label = glue::glue("{city}\n{format(solar, digit = 2)}％")), 
                  stat = "sf_coordinates", data = solar_ratio) + 
  theme_minimal() + theme(axis.title = element_blank())

VIIRS Nighttime Lights（VNL）

データ

NASAとNOAA（アメリカ海洋大気庁）の共同による極軌道衛星システム（Joint Polar-orbiting Satellite System：JPSS）の一環として2011年に打ち上げられた人工衛星が夜間光データを提供しています。この衛星に搭載されているセンサ群（Visible and Infrared Imaging Suite：VIIRS）のDay Night Band（DNB）センサが夜間に観測された光の明るさを記録したものです。

現在このデータは、EOG（Earth Observation Group）によってVIIRS Nighttime Lights（VNL）として公開されており、自由に利用することが可能です。

今回は、Annual VNL V2から2021年の日本周辺のデータ切り出したものを準備したので、以下のリンクからダウンロードしてください。

nighttime_light.tiff

ダウンロードしたファイルを、プロジェクトのdataフォルダに移動してください。

次に、内閣府の県民経済計算のサイトから、県内総生産のデータをダウンロードしてください。

県民経済計算（平成23年度 - 令和3年度）（2008SNA、平成27年基準計数）

一番上の、1.県内総生産（生産側、名目）のExcelファイル（soukatu1.xlsx）をダウンロードし、プロジェクトのdataフォルダに移動してください。

そして、都道府県ポリゴンデータをダウンロードしてください。

都道府県ポリゴンデータ

夜間光データの可視化

夜間光データをterra::rast関数で読み込みます。

# 夜間光データの読み込み
vnl <- rast("data/nighttime_light.tiff")

まずはそのままプロットしてみます。

# 夜間光データのプロット
plot(vnl, range = c(0, 50), col = hcl.colors(100, "Cividis"))

夜間光データの平方根を取ると見た目が良くなることが知られています。

# 夜間光データ（平方根）のプロット
vnl2 <- vnl
values(vnl2) <- sqrt(values(vnl2))

Warning in sqrt(values(vnl2)): NaNs produced

plot(vnl2, range = c(0, sqrt(50)), col = hcl.colors(100, "Cividis"))

夜間光データと県内総生産

夜間光の光量データと地域経済指標との間に相関が見られることが知られています。ここでは、2021年の県内総生産（名目値）と2020年の夜間光データの相関関係を確認します。

まず、県内総生産のデータをreadxl::read_excel関数を使って読み込みます。

# 県民総生産データの読み込み
grp <- read_excel("data/soukatu1.xlsx", skip = 6, n_max = 47,
                  col_types = c("text", "text", rep("skip", 11), "numeric", "skip"),
                  col_names = c("ID", "name", "grp"))

都道府県のポリゴンを使って、都道府県内の夜間光データの合計値を計算します。ここではterra::zonal関数を使います。 fun引数に”sum”を指定すると、合計値を抽出してくれます。

# 夜間光データの都道府県集計
japan <- read_sf("data/prefectures.geojson")
vnl_pref <- zonal(vnl, vect(japan), fun = 'sum') |> 
  rename(vnl = tmp)

夜間光データと県内総生産のデータからなるデータフレームを作成します。

# データフレームの作成
dat <- grp |> bind_cols(vnl_pref)
glimpse(dat)

Rows: 47
Columns: 4
$ ID   <chr> "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11",…
$ name <chr> "北海道", "青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県", "茨城県", "栃木県", "群…
$ grp  <dbl> 20540923, 4464610, 4701411, 9649597, 3545316, 4282525, 7844733, 1…
$ vnl  <dbl> 270707.56, 40543.29, 29884.65, 64071.10, 24471.00, 29114.16, 4919…

横軸に夜間光データ、縦軸に県内総生産をとってプロットしてみます。縦軸・横軸ともに常用対数軸にしていることに注意してください。

東京都と北海道をのぞいて、かなりよくフィットしていることがわかります。東京都は夜間光に比べて県内総生産が大きく、逆に北海道は夜間光に比べて県内総生産が小さいことがわかります。その他、沖縄県と青森県も、やや回帰直線から外れていますが、その他の府県は、かなりよく当てはまっており、全体で見ると夜間光と県内総生産に強い相関が見られることが分かります。

# 夜間光と県内総生産の散布図＋回帰直線
ggplot(dat, aes(x = vnl, y = grp)) +
  geom_point() + geom_smooth(method = "lm") +
  geom_text_repel(aes(label = name), min.segment.length = 0) +
  scale_x_log10() + scale_y_log10() + 
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

Warning: ggrepel: 19 unlabeled data points (too many overlaps). Consider
increasing max.overlaps