6 ラスタとベクタの相互作用

必須パッケージ

この章では、以下のパッケージが必要である。

6.1 イントロダクション

この章では、Chapter 2 で紹介したラスタとベクタの地理データモデル間の相互作用に焦点を当てる。 主に4つの技法を紹介する。 最初は、ベクタオブジェクトを使用したラスタの切り落とし (crop) とマスク (mask) から始める(Section 6.2)。 次に、さまざまな種類のベクタデータを使ってラスタ値を抽出する(Section 6.3)。 最後は、ラスタベクタ変換である(Section 6.4 と Section 6.5)。 以上の概念を、実世界でどのように応用できるかを理解できるよう、これまでの章で使用したデータを用いて実際に操作していく。

6.2 ラスタの切り落とし(crop)

多くの地理データプロジェクトでは、リモートセンシング画像(ラスタ)や行政境界線(ベクタ)など、さまざまなソースからのデータを統合している。 入力されたラスタデータセットの範囲は、対象地域よりも大きいことがよくある。 この場合、入力データの空間的な広がりを統一するために、ラスタの切り落とし (crop) やマスク (mask) が有効である。 この2つの処理により、オブジェクトのメモリ使用量と、その後の解析に必要な計算資源が削減される。ラスタデータを含む魅力的なマップを作成する前に必要な前処理工程となる場合がある。

ここでは、2つのオブジェクトを使ってラスタ切り落としを説明する。

  • SpatRaster のオブジェクト srtm は、米国ユタ州南西部の標高(海抜メートル)を表す。
  • ベクタ (sf) オブジェクト zion は、ザイオン国立公園 (Zion National Park) を表す。

ターゲットと切り取りオブジェクトの両方が同じ投影である必要がある。 したがって、以下のコードは Chapter 2 でインストールされた spDataLarge パッケージからデータセットを読み込むだけでなく、zion を「再投影」している(この話題は Chapter 7 で取り上げている)。

srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))

srtm のラスタを切り出すために terra パッケージの crop() を使用する。 この関数は、その第 2 引数に渡されたオブジェクトの範囲に基づいて、その第 1 引数に渡されたオブジェクトの長方形の範囲に縮小する。 この機能を、Figure 6.1 (B) を生成する以下のコマンドで示している

srtm_cropped = crop(srtm, zion)

crop() に関連するものとして、terra 関数 mask() があり、これは第2引数に渡されたオブジェクトの境界外の値を NA に設定するものである。 したがって、次のコマンドは、ザイオン国立公園の境界の外側のすべてのセルをマスクする(Figure 6.1 (C))。

srtm_masked = mask(srtm, zion)

crop()mask() は、ほとんどの場合で一緒に使う。 (a)ラスタの範囲を目的の領域に限定し、(b)領域外の値をすべて NA に置き換える。

srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)

mask() の設定を変更すると、異なる結果が得られる。 例えば、updatevalue = 0 を設定すると、国立公園外のすべてのピクセルが0に設定される。 inverse = TRUE を設定すると、公園の境界の内側をすべてマスクする (詳細は ?mask を参照)(Figure 6.1 (D))。

srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
ラスタクロップ、ラスタマスクの説明図。

FIGURE 6.1: ラスタクロップ、ラスタマスクの説明図。

6.3 ラスタ抽出

ラスタ抽出は、地理的(通常はベクタ)な「範囲選択」オブジェクトに基づいて、特定の位置の「ターゲット」ラスタに関連する値を識別して返す処理である。 結果は、使用する範囲選択の種類(点、線、ポリゴン)と、terra::extract() 関数に渡される引数に依存する。この関数を使って、ラスタ抽出のデモを行う。 ラスタ抽出の逆、つまりベクタオブジェクトに基づいてラスタセル値を割り当てるのがラスタ化で、Section 6.4 で説明する。

基本的な例として、ラスタセルの特定のの値を抽出してみよう。 そのために、ザイオン国立公園内の30カ所のサンプルを収録した zion_points を使用する(Figure 6.2)。 次のコマンドは、srtm から標高値を抽出し、各ポイントの ID(ベクタの行ごとに 1 つの値)と関連する srtm の値を含むデータフレームを作成する。 さて、出来上がったオブジェクトを cbind() 関数で zion_points データセットに追加してみよう。

data("zion_points", package = "spDataLarge")
elevation = terra::extract(srtm, zion_points)
zion_points = cbind(zion_points, elevation)
ラスタ抽出に使用した点の位置。

FIGURE 6.2: ラスタ抽出に使用した点の位置。

ラスタ抽出は、範囲選択でも機能する。 そして、線に接するラスタセルごとに1つの値を抽出する。 しかし、抽出されたラスタ値の各ペア間の距離を正しく取得することが難しいため、線の断片に沿った値を得るための線抽出アプローチは推奨されていない。

この場合、線を多くの点に分割し、その点の値を抽出するのが良い方法である。 これを示すために、以下のコードでは、Figure 6.3 (A) に示したザイオン国立公園の北西から南東に向かう直線、zion_transect を作成する(ベクタデータモデルについての復習は Section 2.2 を参照)。

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) |>
  st_linestring() |> 
  st_sfc(crs = crs(srtm)) |>
  st_sf(geometry = _)

線の範囲選択から高さを抽出することの有用性は、ハイキングの計画を立てることを想像してみるとよくわかる。 以下に示す方法は、ルートの「標高プロファイル」を提供し(線は直線である必要はない)、長い上り坂による所要時間を見積もるのに便利である。

まず、各断片に固有の id を追加する。 次に、st_segmentize() 関数を使って、与えられた密度(dfMaxLength)で線に沿って点を追加し、st_cast() で点に変換することができる。

zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")

これで大きな点の集合ができたので、断片の最初の点と、それ以降の各点との距離を導き出してみたい。 このケースでは、1つの断片しかないが、原理的には、このコードはいくつの断片でも動作するはずである。

zion_transect = zion_transect |> 
  group_by(id) |> 
  mutate(dist = st_distance(geometry)[, 1]) 

最後に、断片の各ポイントの標高値を抽出し、この情報をメインオブジェクトに結合することができる。

zion_elev = terra::extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)

その結果、zion_transect、Figure 6.3 (B) に示すように、標高プロファイルを作成することができる。

ラスタ抽出に使用した線の位置(左)と、その線に沿った標高(右)。

FIGURE 6.3: ラスタ抽出に使用した線の位置(左)と、その線に沿った標高(右)。

ラスタ抽出のための地理ベクタオブジェクトの最後のタイプは、ポリゴンである。 線と同様に、ポリゴンも1ポリゴンあたり多くのラスタ値を返す傾向がある。 これは以下のコマンドで示され、ID (ポリゴンの行番号)と srtm (関連する標高値)の列名を持つデータフレームが生成される。

zion_srtm_values = terra::extract(x = srtm, y = zion)

このような結果を利用して、ポリゴンごとのラスタ値の要約統計量を生成することで、例えば、単一の地域を特徴付けることや、多くの地域を比較することができる。 要約統計の生成は以下のコードで示されている。このコードは、ザイオン国立公園の標高値の要約統計を含むオブジェクト zion_srtm_df を作成する(Figure 6.4 (A) を参照)。

group_by(zion_srtm_values, ID) |> 
  summarize(across(srtm, list(min = min, mean = mean, max = max)))
#> # A tibble: 1 × 4
#>      ID srtm_min srtm_mean srtm_max
#>   <dbl>    <int>     <dbl>    <int>
#> 1     1     1122     1818.     2661

前述のコードチャンクは、Chapter 3 で説明されているように、ポリゴン ID ごとのセル値の要約統計を提供するために tidyverse を使用した。 その結果、例えば公園の最高標高が海抜約 2,661 m であることなど、有用な要約が得られる(標準偏差など、他の要約統計もこの方法で計算できる)。 この例ではポリゴンが1つしかないので、1行のデータフレームが返されるが、このメソッドは複数の範囲選択ポリゴンが使用されている場合に動作する。

同様のアプローチは、ポリゴン内のカテゴリ的なラスタ値の出現をカウントする場合にも有効である。 これは、Figure 6.4 (B) の spDataLarge パッケージの土地被覆データセット ( nlcd ) を使って説明され、以下のコードで実証されている。

nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
zion2 = st_transform(zion, st_crs(nlcd))
zion_nlcd = terra::extract(nlcd, zion2)
zion_nlcd |> 
  group_by(ID, levels) |>
  count()
#> # A tibble: 7 × 3
#> # Groups:   ID, levels [7]
#>      ID levels         n
#>   <dbl> <fct>      <int>
#> 1     1 Developed   4205
#> 2     1 Barren     98285
#> 3     1 Forest    298299
#> 4     1 Shrubland 203701
#> # … with 3 more rows
連続(左)とカテゴリ(右)のラスタ抽出に使用した範囲。

FIGURE 6.4: 連続(左)とカテゴリ(右)のラスタ抽出に使用した範囲。

terra パッケージはポリゴン内のラスタ値を高速に抽出するが、 extract() は大きなポリゴンデータセットを処理する際のボトルネックになることがある。 exactextractr パッケージは、 exact_extract() 関数を通してピクセル値を抽出するための 超高速の代替手段 を提供する。 また、exact_extract() 関数は、デフォルトで、ポリゴンによってオーバーラップされた各ラスタセルの割合を計算し、より正確である(詳細については、以下の注を参照)。

ポリゴンは通常不規則な形状をしているため、ラスタのセルの一部にしか重ならないことがある。 より詳細な結果を得るために、extract() 関数には exact と呼ばれる引数がある。 exact = TRUE とすると、出力データフレームに fraction という列が追加され、ポリゴンによってカバーされる各セルの割合が格納される。 これは、連続ラスタの加重平均や、カテゴリラスタのより正確なカバレッジを計算するのに便利である。 デフォルトでは、この操作はより多くの計算を必要とするため、FALSE に設定されている。 exactextractr::exact_extract() 関数は、常に各セルにおけるポリゴンの被覆率を計算する。

6.4 ラスタ化

ラスタ化とは、ベクタオブジェクトをラスタオブジェクトに変換して表現することである。 通常、出力されたラスタは定量的な解析(地形の解析など)やモデリングに利用される。 Chapter 2 で見たように、ラスタデータモデルには、ある種の手法を助長するような特徴がある。 さらに、ラスタライゼーションは地理的なデータ集計の一種と考えることができ、結果として得られる値はすべて同じ空間分解能を持つため、データセットを簡素化することができる。

terra パッケージには、この作業を行うための関数 rasterize() が含まれている。 その最初の2つの引数は、x、ラスタ化されるベクタオブジェクトと、y、出力の範囲、解像度、CRSを定義する「テンプレートラスタ」オブジェクトである。 入力ラスタの地理的解像度が低すぎる(セルサイズが大きすぎる)と、ベクタデータの地理的変動を完全に見逃す可能性があり、高すぎる場合は計算時間がかかりすぎる可能性がある。 適切な地理的解像度を決定する際に従うべき単純なルールはなく、結果の使用目的によって大きく左右される。 例えば、ラスタ化の出力を既存のラスタに合わせる必要がある場合など、ターゲット解像度がユーザーに課されることがよくある。

ラスタ化を実演するために、入力ベクタデータ cycle_hire_osm_projected (ロンドンの自転車レンタルポイントに関するデータセットを Figure 6.5 (A)に図示)と同じ範囲とCRS、空間解像度1000メートルのテンプレートラスタを使用することにする。

cycle_hire_osm = spData::cycle_hire_osm
cycle_hire_osm_projected = st_transform(cycle_hire_osm, "EPSG:27700")
raster_template = rast(ext(cycle_hire_osm_projected), resolution = 1000,
                       crs = st_crs(cycle_hire_osm_projected)$wkt)

ラスタ化は非常に柔軟な操作で、結果はテンプレートとなるラスタの性質だけでなく、入力ベクタの種類(ポイント、ポリゴンなど)や、rasterize() 関数が取るさまざまな引数に依存する。

この柔軟性を説明するために、3つの異なるアプローチでラスタ化を試みる。 まず、サイクルハイヤーポイントの有無を表すラスタ(有無ラスタと呼ばれる)を作成する。 この場合、rasterize() は、xy (前述のベクタとラスタのオブジェクト)に加えて、field で指定された空でないすべてのセルに転送する値(図示の結果 Figure 6.5 (B))を1つだけ引数として要求することになる。

ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template,
                       field = 1)

fun 引数は、近接した複数の観測値をラスタオブジェクトの関連セルに変換するために使用される要約統計量を指定する。 デフォルトでは、fun = "last" が使用されるが、fun = "length" などの他のオプションも使用できる。この場合、各グリッドセル内のサイクルレンタルポイントの数をカウントする(この操作の結果は、Figure 6.5 (C) に示される)。

ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template, 
                       fun = "length")

新しい出力である ch_raster2 は、各グリッドセル内のサイクルハイヤーポイントの数を示している。 レンタルサイクルの拠点は、capacity 変数で記述される自転車の数が異なるため、各グリッドセルの収容台数 (capacity) はどの程度なのかという疑問が生じる。 これを計算するためには、フィールド("capacity") を sum することが必要で、その結果、Figure 6.5 (D) のような出力が得られる。以下のコマンドで計算する(mean など他の要約関数も使用できる)。

ch_raster3 = rasterize(vect(cycle_hire_osm_projected), raster_template, 
                       field = "capacity", fun = sum)
点のラスタ化例。

FIGURE 6.5: 点のラスタ化例。

また、カリフォルニア州のポリゴンと境界線をベースにしたデータセット(下記作成)は、線のラスタ化を表している。 ポリゴンオブジェクトをマルチリンストリングにキャストした後、0.5度の分解能を持つテンプレートラスタを作成する。

california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = rast(ext(california), resolution = 0.5,
                        crs = st_crs(california)$wkt)

線またはポリゴンのラスタ化を考慮する場合、有用な追加引数の1つは touches である。 デフォルトでは FALSE であるが、TRUE に変更すると、線またはポリゴンの境界で接触しているすべてのセルが値を得る。 touches = TRUE による線分ラスタ化を以下のコードで実行する ( Figure 6.6 (A))。

california_raster1 = rasterize(california_borders, raster_template2,
                               touches = TRUE)

ポリゴンのラスタ化と比較すると、touches = FALSE がデフォルトで、Figure 6.6 (B) に示すように、範囲選択ポリゴン内に中心点があるセルのみが選択されることになる。

california_raster2 = rasterize(california, raster_template2) 
線とポリゴンのラスタ化の例。

FIGURE 6.6: 線とポリゴンのラスタ化の例。

6.5 空間ベクタ化

空間ベクタ化は、ラスタ化(Section 6.4)と対で、方向が逆になる。 空間的に連続したラスタデータを、点、線、ポリゴンなどの空間的に離散したベクタデータに変換する。

表現に注意。 R では、単にベクトル化と言った場合、for ループなどを 1:10 / 2 のように置き換えることができるようにすることを指す(Wickham (2019) 参照)。

最も単純なベクタ化は、ラスタセルの中心を点に変換することである。 as.points() は、NA 以外のすべてのラスタグリッドセルに対して実行する(Figure 6.7)。 なお、ここでは、st_as_sf() を使って、結果のオブジェクトを sf クラスに変換することもしている。

elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_point = as.points(elev) |> 
  st_as_sf()
elev オブジェクトのラスタ表現と点表現。

FIGURE 6.7: elev オブジェクトのラスタ表現と点表現。

空間ベクタ化のもう一つの一般的なタイプは、例えば連続した高さや温度の線(等温線)を表す等高線の作成である。 人工ラスタ elev は平行線を生成するため、実世界のデジタル標高モデル(DEM)を使用しよう(読者への課題:これを検証し、なぜこうなるのかを説明しなさい)。 等高線は terra 関数 as.contour() で作成することができる。この関数は、 filled.contour() のラッパーである(図示せず)。

dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
cl = as.contour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)

contour() , rasterVis::contourplot() , tmap::tm_iso() などの関数で、既存のプロットに等高線を追加することもできる。 Figure 6.8 に示すように、等高線 (isoline) にはラベルを付けることができる。

モンゴル山南麓の等高線を重ねた陰影付き DEM。

FIGURE 6.8: モンゴル山南麓の等高線を重ねた陰影付き DEM。

ベクタ化の最後のタイプとして、ラスタをポリゴンに変換しよう。 これには、terra::as.polygons() を使う。各ラスタセルを5つの座標からなるポリゴンに変換し、そのすべてがメモリに保存される(ラスタがベクタと比較して高速である理由を説明する!)。

以下では、grain オブジェクトをポリゴンに変換し、その後、同じ属性値を持つポリゴン間の境界を解消することで説明している(as.polygons()dissolve の引数も参照)。

grain = rast(system.file("raster/grain.tif", package = "spData"))
grain_poly = as.polygons(grain) |> 
  st_as_sf()

grain 粒状データセットのポリゴンは、長方形のピクセルをつなぐことで定義される直方体の境界を持つ。 ポリゴンを滑らか (smooth) にするために、Chapter 5 のパッケージ smoothr を使用することができる。 平滑化処理はポリゴン境界の鋭いエッジを除去するため、平滑化ポリゴンは元のピクセルとは空間的範囲が正確に同じにはならない(例は、smoothr website を参照)。 そのため、平滑化されたポリゴンをさらに解析に使用する場合は注意が必要です。

ラスタ(左)をポリゴン(dissolve = FALSE、中央)と集約ポリゴン(dissolve = TRUE、右)にベクトル化した図。

FIGURE 6.9: ラスタ(左)をポリゴン(dissolve = FALSE、中央)と集約ポリゴン(dissolve = TRUE、右)にベクトル化した図。

6.6 演習

Some of the exercises use a vector (zion_points) and raster dataset (srtm) from the spDataLarge package. They also use a polygonal ‘convex hull’ derived from the vector dataset (ch) to represent the area of interest:

library(sf)
library(terra)
library(spData)
zion_points_path = system.file("vector/zion_points.gpkg", package = "spDataLarge")
zion_points = read_sf(zion_points_path)
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
ch = st_combine(zion_points) |>
  st_convex_hull() |> 
  st_as_sf()

E1. Crop the srtm raster using (1) the zion_points dataset and (2) the ch dataset. Are there any differences in the output maps? Next, mask srtm using these two datasets. Can you see any difference now? How can you explain that?

E2. Firstly, extract values from srtm at the points represented in zion_points. Next, extract average values of srtm using a 90 buffer around each point from zion_points and compare these two sets of values. When would extracting values by buffers be more suitable than by points alone?

  • Bonus: Implement extraction using the exactextractr package and compare the results.

E3. Subset points higher than 3100 meters in New Zealand (the nz_height object) and create a template raster with a resolution of 3 km for the extent of the new point dataset. Using these two new objects:

  • Count numbers of the highest points in each grid cell.
  • Find the maximum elevation in each grid cell.

E4. Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 by 6 km) and plot the result.

  • Resample the lower resolution raster back to the original resolution of 3 km. How have the results changed?
  • Name two advantages and disadvantages of reducing raster resolution.

E5. Polygonize the grain dataset and filter all squares representing clay.

  • Name two advantages and disadvantages of vector data over raster data.
  • When would it be useful to convert rasters to vectors in your work?