4 空間データ操作

必須パッケージ

  • この章では、Chapter 3 で使用したものと同じパッケージが必要である。
  • さらに、Section 4.3 では、以下の二つのデータセットも読み込む必要がある。
elev = rast(system.file("raster/elev.tif", package = "spData"))
grain = rast(system.file("raster/grain.tif", package = "spData"))

4.1 イントロダクション

ジオコンピュテーションにおいて、ベクタデータセット間の空間結合やラスタデータセットのローカルおよびフォーカル演算などの空間演算は重要な要素である。 この章では、空間オブジェクトがその位置と形状に基づいて、さまざまな方法で変更できることを紹介する。 空間的な操作の多くは、非空間的(属性的)な操作に相当するため、前の章で示したデータセットの部分集合や結合といった概念がここでも適用できる。 これは特にベクタ操作に当てはまる。ベクタ属性の操作に関する Section 3.2 は、空間的な対応である空間部分集合(Section 4.2.1 で取り上げている)を理解するための基礎となるものである。 空間結合(Section 4.2.4)と属性集計(Section 4.2.6)には、前の章で説明した非空間的な対応関係がある。

しかし、空間演算は非空間演算と異なる点がいくつもある。 例えば空間結合は、ターゲットデータセットと交差する、または一定の距離内にあるモノのマッチングなど、結合方法は多数ある。一方、前章の Section 3.2.4 で説明した属性結合は、1つの方法でしかできない(ただし fuzzyjoin パッケージのドキュメントで説明した、ファジー結合を使う場合は別)。 オブジェクト間の空間的関係のさまざまなタイプ(intersect と disjoint を含む)については、Section 4.2.2 で説明する。 空間オブジェクトのもう一つのユニークな側面は距離である。すべての空間オブジェクトは空間を通じて関連しており、距離計算を行うことでこの関係の強さを調べることができる。これは、Section 4.2.8 のベクタデータで説明する。

ラスタの空間演算には、部分集合(Section 4.3.1)、および複数のラスタ「タイル」を単一のオブジェクトに統合する(Section 4.3.8)方法がある。 マップ代数は、ラスタのセルの値を、周囲のセル値を参照して、あるいは参照せずに変更する操作を対象とする。 多くのアプリケーションに不可欠なマップ代数の概念は、Section 4.3.2 で紹介している。ローカル、フォーカル、ゾーンのマップ代数演算については、それぞれ Section 4.3.3、Section 4.3.4、Section 4.3.5 のセクションで解説している。ラスタデータセット全体を表す要約統計量を生成するグローバルマップ代数操作と、ラスタの距離計算については、Section 4.3.6 で説明する。 演習の前の最後の章では、2つのラスタデータセットを合成する方法について説明し、再現可能な例を挙げて実演している。

2つの空間オブジェクトを使用する空間操作は、両方のオブジェクトが同じ座標参照系を持つことに依存することに留意する。この点については、Section 2.4 で紹介し、Chapter 7 でさらに詳しく解説する。

4.2 ベクタデータに対する空間演算

ここでは、sf パッケージのシンプルフィーチャとして表現されたベクタ地理データに対する空間演算の概要を説明する。 Section 4.3 は、terra パッケージのクラスと関数を使用したラスタデータセットの空間演算を紹介する。

4.2.1 空間部分集合

空間部分集合とは、空間オブジェクトを取り出し、別のオブジェクトと空間的に関連するフィーチャだけを含む新しいオブジェクトを返す処理である。 属性部分集合(Section 3.2.1 で説明)と同様に、sf データフレームの部分集合は、角括弧を使用して作成できる ([) を使い、x[y, , op = st_intersects] という文法を使う。ここで x は行の部分集合が返される sf オブジェクト、y は「部分集合・オブジェクト」、op = st_intersects は部分集合を行うために使われる位相関係(二項述語 binary predicate としても知られている)を指定するオプションの引数である。 op 引数がないときに使用されるデフォルトの位相関係は、st_intersects()である。つまり、コマンド x[y, ]x[y, , op = st_intersects] と同じであり x[y, , op = st_disjoint] とは異なる(これらのトポロジカル関係の意味と他のトポロジカル関係は次のセクションで説明する)。 tidyversefilter() 関数も使用できるが、この方法は、以下の例で見るように、より冗長である。

空間部分集合を示すために、spData パッケージの nznz_height データセットを使用する。これは、ニュージーランドの 16 の主要地域と 101 の最高地点に関する地理データをそれぞれ含み(Figure 4.1)、投影座標系で表示されるものである。 次のコードでは、Canterbury を表すオブジェクトを作成し、空間部分集合を使用して、対象地域のすべての High Point を返す。

canterbury = nz |> filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
赤い三角形はニュージーランドの 101 の High Point を表し、中央のカンタベリー地域付近に集まっている(左)。カンタベリーのポイントは、`[`部分集合演算子で作成された(グレーでハイライト、右)。

FIGURE 4.1: 赤い三角形はニュージーランドの 101 の High Point を表し、中央のカンタベリー地域付近に集まっている(左)。カンタベリーのポイントは、[部分集合演算子で作成された(グレーでハイライト、右)。

属性による部分集合と同様に、コマンド x[y, ]nz_height[canterbury, ] と同等)は、ソース オブジェクトの内容を使用して、ターゲットのフィーチャの部分集合を返す。 しかし、y がクラス logical または integer のベクトルである代わりに、空間部分集合では xy の両方が地理オブジェクトでなければならない。 具体的には、この方法で空間部分集合に使用されるオブジェクトは、クラス sf または sfc を持つ必要がある。nznz_height は共に地理ベクタデータフレームで、クラス sf を持つ。操作の結果、ターゲット nz_height オブジェクトの中で canterbury 地域と交差する(この場合はその中にある点)フィーチャを表す別の sf オブジェクトが返される。

空間部分集合には、ターゲットオブジェクトのフィーチャが選択される部分集合オブジェクトとどのような空間的関係を持たなければならないかを決める、様々な位相関係を用いることができる。 これには、Section 4.2.2 で見るように、touchescrosseswithin が含まれる。 デフォルトの設定 st_intersects は、ソース「部分集合」オブジェクトに touchescrosseswithin するターゲット内のフィーチャを返す「全て」トポロジー関係である。 上記のように、op = 引数で別の空間演算子を指定することができる。次のコマンドでは、st_intersects() の逆で、カンタベリーと交差しない点を返す( Section 4.2.2 を参照)。

nz_height[canterbury, , op = st_disjoint]
空白の引数(, , で示される)は、[ が、sf オブジェクトの3番目の引数である op を強調するために含まれていることに注意。 これを使うと、部分集合操作をいろいろな方法で変更することができる。 例えば、nz_height[canterbury, 2, op = st_disjoint] は同じ行を返すが、2番目の属性列のみを含む (詳細は sf:::`[.sf` and the ?sf)。

ベクタデータの空間部分集合についてこれだけ知っておけば多くの応用例を使うことができる。 st_intersects()st_disjoint() 以外のトポロジカル・リレーションをすぐに学びたいのであれば、残りを飛ばして次の章(Section 4.2.2)に飛んでも構わない。 ここからは、その他の部分集合化の方法などの詳細について説明する。

空間部分集合を行うもう一つの方法は、位相演算子によって返されるオブジェクトを使用することである。 これらのオブジェクトは、それ自体、例えば、連続する領域間の関係のグラフネットワークを探索する際に有用であるが、以下のコードチャンクで示されるように、部分集合にも使用できる。

sel_sgbp = st_intersects(x = nz_height, y = canterbury)
class(sel_sgbp)
#> [1] "sgbp" "list"
sel_sgbp
#> Sparse geometry binary predicate list of length 101, where the
#> predicate was `intersects'
#> first 10 elements:
#>  1: (empty)
#>  2: (empty)
#>  3: (empty)
#>  4: (empty)
#>  5: 1
#>  6: 1
#>  7: 1
#>  8: 1
#>  9: 1
#>  10: 1
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]

上記のコードチャンクは、クラス sgbp のオブジェクト(疎な幾何学二項述語、空間演算における長さ x のリスト)を作成し、それを論理ベクタ sel_logicalTRUEFALSE の値のみを含み、dplyr のフィルタ関数でも使用できるもの)に変換している。 関数 lengths() は、nz_height のどのフィーチャが y任意の物体と交差しているかを特定する。 この場合、1が最も大きな値であるが、より複雑な操作を行う場合には、例えば、ソースオブジェクトの2つ以上のフィーチャと交差するフィーチャのみを部分集合するような方法を用いることができる。

注意: 論理的な出力を返す別の方法として、st_intersects() などの演算子で sparse = FALSE (「疎行列ではなく密行列を返す」という意味) をセットする方法がある。例えば、st_intersects(x = nz_height, y = canterbury, sparse = FALSE)[, 1] というコマンドは、sel_logical と同じ出力を返す。 注意: sgbp オブジェクトを含むソリューションは、多対多のオペレーションに対応し、より低いメモリ要件で動作するため、より一般的な方法である。

sf オブジェクトと dplyr データ操作コードの互換性を高めるために作成された sf 関数 st_filter() でも同じ結果を得ることができる。

canterbury_height3 = nz_height |>
  st_filter(y = canterbury, .predicate = st_intersects)

この時点で、(行名以外は)同じバージョンの canterbury_height が3つある。1つは [ 演算子を用いて作成し、もう1つは中間選択オブジェクトを介して作成し、最後は sf の便利な関数 st_filter() を用いて作成した。 次のセクションでは、二つのフィーチャが空間的に関連しているかどうかを識別するために使用できる、二項述語としても知られている、さまざまなタイプの空間的関係性を探る。

4.2.2 トポロジー関係

トポロジー関係は、オブジェクト間の空間的な関係を表す。 「二項位相関係」(binary topological relationships)とは、2次元以上の点(一般的には点、線、ポリゴン)の順序集合で定義される2つの物体間の空間関係について論理的に記述したもの(答えは TRUEFALSE しかない)である (Egenhofer and Herring 1990)。 このように言うと、かなり抽象的に聞こえるだろうが、実際、位相関係の定義と分類は、1966年に初めて書籍として出版された数学的基礎に基づいている (Spanier 1995)。 代数的位相幾何学の分野は21世紀まで続いている (Dieck 2008)

トポロジー関係は数学的な起源を持つが、一般的な空間的関係をテストするためによく使われる関数を視覚化することで、直感的に理解することが可能である。 Figure 4.2 は、様々なジオメトリペアとその関連性を示している。 Figure 4.2 の3番目と4番目のペア(左から右、そして下)は、ある関係では順序が重要であることを示している。関係 equalsintersectscrossestouchesoverlaps は対称であり、function(x, y) が真なら function(y, x) も真となるが、 containswithin など幾何学の順序が重要である関係は、そうではない。 各ジオメトリペアには、次節で説明する FF2F11212 のような「DE-9IM」文字列があることを確認しておこう。

Egenhofer and Herring (1990)の Figure 1 と Figure 2 を参考にした、ベクトル幾何学間のトポロジー関係。関数(x, y)が真となる関係が、各ジオメトリのペアについて印刷されており、xはピンク、yは青で表されている。各ペアの空間的関係の性質は、Dimensionally Extended 9-Intersection Model 文字列で記述されている。

FIGURE 4.2: Egenhofer and Herring (1990)の Figure 1 と Figure 2 を参考にした、ベクトル幾何学間のトポロジー関係。関数(x, y)が真となる関係が、各ジオメトリのペアについて印刷されており、xはピンク、yは青で表されている。各ペアの空間的関係の性質は、Dimensionally Extended 9-Intersection Model 文字列で記述されている。

sf では、異なる種類の位相関係をテストする関数を「二項述語」と呼ぶ。これは、コマンド vignette("sf3")(訳注:日本語版)、およびヘルプページ ?geos_binary_pred で見ることができる。 位相関係が実際にどのように機能するかを見るために、Figure 4.2 で説明した関係を基に、前の章(Section 2.2.4)で学んだベクタギオメトリの表現方法の知識を統合して、簡単な再現性のある例を作ってみよう。 なお、ポリゴンの頂点の座標(x、y)を表す表形式のデータを作成するために、R の基本関数 cbind() を使って、座標点を表す行列、POLYGON、そして最後に sfc オブジェクトを作成する(Chapter 2 で説明)。

polygon_matrix = cbind(
  x = c(0, 0, 1, 1,   0),
  y = c(0, 1, 1, 0.5, 0)
)
polygon_sfc = st_sfc(st_polygon(list(polygon_matrix)))

次のコマンドで、空間的な関係を示す追加の形状を作成する。これらの形状は、上で作成したポリゴンの上にプロットすると、Figure 4.3 に示すように、互いに空間的に関連するようになる。 関数 st_as_sf() と引数 coords を使って、座標を表す列を含むデータフレームから、点を含む sf オブジェクトに効率的に変換していることに注目してみよう。

line_sfc = st_sfc(st_linestring(cbind(
  x = c(0.4, 1),
  y = c(0.2, 0.5)
)))
# create points
point_df = data.frame(
  x = c(0.2, 0.7, 0.4),
  y = c(0.1, 0.2, 0.8)
)
point_sf = st_as_sf(point_df, coords = c("x", "y"))
点、線、ポリゴンのオブジェクトを配置し、トポロジー関係を表現。

FIGURE 4.3: 点、線、ポリゴンのオブジェクトを配置し、トポロジー関係を表現。

簡単なクエリを作ってみよう。point_sf の点のうち、ポリゴン polygon_sfc と何らかの形で交差しているものはどれか? この問題は、見れば答えることができる(点1は接していて、点3は中にある)。 この質問には、空間についての関数 st_intersects() を用いて、次のように答えることができる。

st_intersects(point_sf, polygon_sfc)
#> Sparse geometry binary predicate... `intersects'
#>  1: 1
#>  2: (empty)
#>  3: 1

その結果は、直感と一致するはずである。 1点目と3点目は真の値(1)、2点目は偽の値(空のベクタで表される)が返され、ポリゴンの境界の外にある。 予想外のは、その結果がベクタのリストという形になっていることだ。 この疎行列 (sparse matrix) 出力は、関係が存在する場合にのみ登録され、複数フィーチャに対する位相幾何学的操作のメモリ要件を軽減する。 前節で見たように、TRUE または FALSE の値からなる密な行列 (dense matrix) が返されるのは、sparse = FALSE の時である。

st_intersects(point_sf, polygon_sfc, sparse = FALSE)
#>       [,1]
#> [1,]  TRUE
#> [2,] FALSE
#> [3,]  TRUE

上記の出力では、各行がターゲット(引数 x)オブジェクトの特徴量、各列が選択オブジェクト(y)の特徴量を表している。 この場合、y オブジェクト polygon_sfc にはフィーチャが1つしかないので、Section 4.2.1 で見たように部分集合に使える結果は1列だけである。

st_intersects() は、フィーチャが接触しているだけの場合でも TRUE を返す。intersects は、Figure 4.2 に示されているように、多くのタイプの空間的関係を識別する「全て捕まえる」トポロジー操作である。 Figure 4.2 より限定的な質問としては、どの点がポリゴン内にあるか、どのフィーチャが y と共有の境界線上にあるか、またはそれを含んでいるか、などがある。 こうした問いには、次のように答えることができる(結果は示していない)。

st_within(point_sf, polygon_sfc)
st_touches(point_sf, polygon_sfc)

点 1 は境界ポリゴンに接触しているが、境界ポリゴン内にはない。 st_intersects() の反対は st_disjoint() で、これは選択したオブジェクトと空間的に何ら関係のないオブジェクトだけを返す(注:[, 1] は結果をベクトルに変換する)。

st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1]
#> [1] FALSE  TRUE FALSE

関数 st_is_within_distance() は、選択オブジェクトにほぼ接触しているフィーチャを検出する。この関数には、さらに dist という引数がある。 ターゲットオブジェクトが選択されるまでに必要な距離を設定することができる。 点 2 は polygon_sfc の最も近い頂点から 0.2 単位以上離れているが、距離を 0.2 に設定すると、まだ選択されていることに注意されたい。 これは、距離が最も近い辺まで測定されるためで、この場合、Figure 4.3 の点2の真上にあるポリゴンの部分である。 (点 2 とポリゴンの実際の距離が 0.13 であることは、コマンド st_distance(point_sf, polygon_sfc) で確認できる。) ‘is within distance’ という二値空間述語は以下のコードチャンクで示され、その結果、すべての点がポリゴンから 0.2 単位以内にあることが示される。

st_is_within_distance(point_sf, polygon_sfc, dist = 0.2, sparse = FALSE)[, 1]
#> [1] TRUE TRUE TRUE
位相関係を計算する関数では、空間インデックスを用いることで、空間クエリの性能を大幅に向上させることができる。 これは、Sort-Tile-Recursive (STR) アルゴリズムにより実現される。 次のセクションで紹介する st_join 関数も、空間インデックスを利用している。 詳しくは、以下のサイトを参照。 https://www.r-spatial.org/r/2017/06/22/spatial-index.html

4.2.3 DE-9IM 文字列

前節で示した二項述語の根底には、DE-9IM(Dimensionally Extended 9-Intersection Model)というものがある。 名前からして暗号のようであり、簡単なテーマではない。 しかし、空間的な関係をよりよく理解するために、学習する価値があるだろう。 さらに、DE-9IM の高度な使い方として、カスタム空間述語を作成することも可能である。 このモデルは当初、発明者によって「2つのフィーチャの境界、内部、外部の交点の次元」を意味する「DE + 9IM」と表示されていたが (Clementini and Di Felice 1995)、現在は「DE-9IM」と表記されている (Shen, Chen, and Liu 2018)

DE-9IM 文字列がどのように機能するかを示すために、Figure 4.2 の最初のジオメトリペアの様々な関連性を見てみよう。 Figure 4.4 は、各オブジェクトの内部、境界、外部のあらゆる組み合わせの交点を示す9交差点モデル(9IM)を示している。最初のオブジェクト x の各コンポーネントを列とし、y の各コンポーネントを行として配置すると、各要素間の交点が強調されたファセット図形が作成される。

Dimensionally Extended 9 Intersection Model (DE-9IM) の仕組みを説明する図。凡例にない色は、異なる構成要素間の重なりを表している。太い線は2次元の交わりを強調する。例えば、オブジェクト x の境界とオブジェクト y の内部の交わりは、中央上部のファセットで示されている。

FIGURE 4.4: Dimensionally Extended 9 Intersection Model (DE-9IM) の仕組みを説明する図。凡例にない色は、異なる構成要素間の重なりを表している。太い線は2次元の交わりを強調する。例えば、オブジェクト x の境界とオブジェクト y の内部の交わりは、中央上部のファセットで示されている。

DE-9IM 文字列は、各タイプの関係の次元から導き出される。 この場合、Figure 4.4 の赤い交点は、 Table 4.1 に示すように、0(点)、1(線)、2(ポリゴン)の次元を持つ。

TABLE 4.1: ジオメトリ x、y の内部、境界、外部の関係を示す表。
Interior (x) Boundary (x) Exterior (x)
Interior (y) 2 1 2
Boundary (y) 1 1 1
Exterior (y) 2 1 2

この行列を「行単位」で一列にすると(つまり、1行目、2行目、3行目の順に連結する)、文字列 212111212 が得られる。 もうひとつの例で、このシステムを紹介する。 Figure 4.2 に示す関係(3列目1行目のポリゴンペア)は、DE-9IM システムでは以下のように定義できる。

  • 大きなオブジェクト x内部y の内部、境界、外部との交点は、それぞれ 2、1、2 の寸法を持つ
  • 大きなオブジェクト x境界y の内部、境界、外部との交点はそれぞれ F, F, 1 の次元を持ち、ここで ‘F’ は ‘false’ を意味し、オブジェクトは不連続である
  • x外部y の内部、境界、外部との交点はそれぞれ F、F、2 の寸法を持つ。大きなオブジェクトの外部は y の内部や境界に接触しないが、小さなオブジェクトと大きなオブジェクトの外部は同じ面積をカバーする

これら3つの構成要素を連結すると、文字列 212 , FF1 , FF2 が作成される。 これは、関数 st_relate() で得られた結果と同じである( Figure 4.2 の他の形状がどのように作成されたかは、この章のソースコードを参照されたい)。

xy2sfc = function(x, y) st_sfc(st_polygon(list(cbind(x, y))))
x = xy2sfc(x = c(0, 0, 1, 1,   0), y = c(0, 1, 1, 0.5, 0))
y = xy2sfc(x = c(0.7, 0.7, 0.9, 0.7), y = c(0.8, 0.5, 0.5, 0.8))
st_relate(x, y)
#>      [,1]       
#> [1,] "212FF1FF2"

DE-9IM 文字列を理解することで、新しい二値空間述語を開発することができる。 ヘルプページ ?st_relate では、チェスの駒を利用して、ポリゴンが境界を共有する「クイーン」(queen)と点のみを共有する「ルーク」(rook、将棋でいう飛車と同じ動き)関係に対する関数定義がそれぞれ記載されている。 「クイーン」の関係は、「境界-境界」の関係(Table 4.1 の2列目と2行目のセル、または DE-9IM 文字列の5番目の要素)が空であってはならないという意味で、パターン F***T**** に対応し、「ルーク」の関係では同じ要素が1でなければならない(線形交点を意味する)ことを意味している。 これらは以下のように実装されている。

st_queen = function(x, y) st_relate(x, y, pattern = "F***T****")
st_rook = function(x, y) st_relate(x, y, pattern = "F***1****")

先に作成したオブジェクト x をベースに、新たに作成した関数を用いて、グリッドの中央のマスに対して、どの要素が「クイーン」「ルーク」であるかを以下のように調べることができる。

grid = st_make_grid(x, n = 3)
grid_sf = st_sf(grid)
grid_sf$queens = lengths(st_queen(grid, grid[5])) > 0
plot(grid, col = grid_sf$queens)
grid_sf$rooks = lengths(st_rook(grid, grid[5])) > 0
plot(grid, col = grid_sf$rooks)
9つの形状を持つグリッドの中央の正方形に対する「クイーン」(左)と「ルーク」(右)の関係を見つけるためのカスタムバイナリ空間述語のデモ。

FIGURE 4.5: 9つの形状を持つグリッドの中央の正方形に対する「クイーン」(左)と「ルーク」(右)の関係を見つけるためのカスタムバイナリ空間述語のデモ。

4.2.4 空間結合

2つの非空間データセットを結合するには、Section 3.2.4 で説明されているように、共有の「キー」変数に依存する。 空間データ結合は同じ概念を適用するが、その代わりに前節で説明した空間関係に依存する。 属性データの場合と同様に、結合では、ソースオブジェクト( y )からターゲットオブジェクト( 結合関数の引数 x )に新しい列を追加する。

例えば、地球上にランダムに分布する10個の点があり、そのうちの陸地にある点はどの国のものか? このアイデアを reproducible example に実装することで、地理データを扱うスキルが身に付き、空間結合がどのように機能するかを知ることができる。 出発点は、地表にランダムに散らばる点を作ることである。

set.seed(2018) # 再現できるように seed を設定
(bb = st_bbox(world)) # 世界の境界
#>   xmin   ymin   xmax   ymax 
#> -180.0  -89.9  180.0   83.6
random_df = data.frame(
  x = runif(n = 10, min = bb[1], max = bb[3]),
  y = runif(n = 10, min = bb[2], max = bb[4])
)
random_points = random_df |> 
  st_as_sf(coords = c("x", "y")) |> # 座標を設定
  st_set_crs("EPSG:4326") # 地理 CRC を設定

Figure 4.6 で示したシナリオでは、random_points オブジェクト(左上)には属性データがないのに対し、world (右上)には凡例で示した国名のサンプルを含む属性があることがわかる。 空間結合は、以下のコードチャンクに示すように、st_join() で実装されている。 出力は、random_joined のオブジェクトで、Figure 4.6 (左下)に図示されている。 結合データセットを作成する前に、空間部分集合を用いて、ランダムな点を含む国だけを含む world_random を作成し、結合データセットで返される国名の数が4であることを検証している( Figure 4.6 の右上のパネル参照)。

world_random = world[random_points, ]
nrow(world_random)
#> [1] 4
random_joined = st_join(random_points, world["name_long"])
空間結合の図解。ソースワールドオブジェクト(右上)からランダムポイント(左上)に新しい属性変数が追加され、最後のパネルで表されるデータになる。

FIGURE 4.6: 空間結合の図解。ソースワールドオブジェクト(右上)からランダムポイント(左上)に新しい属性変数が追加され、最後のパネルで表されるデータになる。

デフォルトでは、st_join() は左結合を行う。つまり、結果は y にマッチしない行を含む x の全ての行を含むオブジェクトとなる( Section 3.2.4 を参照)。Inner Join をする場合には、left = FALSE とする。 空間部分集合と同様に、st_join() で使用されるデフォルトの位相演算子は st_intersects() である。これは join 引数を設定することで変更できる(詳細は ?st_join を参照)。 上の例では、ポリゴンレイヤからポイントレイヤへの列の追加を示しているが、ジオメトリの種類に関係なく同じ方法で行える。 このような場合、例えば x にポリゴンが含まれ、それぞれが y の複数のオブジェクトと一致する場合、空間結合では重複するフィーチャが発生するため、y の一致するオブジェクトごとに新しい行が作成される。

4.2.5 非重複結合

2つの地理データセットが接触していなくても、地理的に強い関係がある場合がある。 すでに spData パッケージに含まれているデータセット cycle_hirecycle_hire_osm が良い例となる。 これらをプロットすると、Figure 4.7 に示すように、しばしば密接に関連しているが、接触していないことがわかる。このベースバージョンは、以下のコードで作成される。

plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")

以下のように、同じ st_intersects() な点があるかどうかを確認することができる。

any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
#> [1] FALSE

FIGURE 4.7: 公式データ(青)とOpenStreetMapのデータ(赤)に基づく、ロンドンにおける自転車レンタルポイントの空間分布。

cycle_hire_osm の変数 capacitycycle_hire に含まれる公式の「ターゲット」データに結合する必要があるとする。 このとき、オーバーラップしない結合が必要である。 最も簡単な方法は、トポロジカル演算子 st_is_within_distance() を使用することである。

sel = st_is_within_distance(cycle_hire, cycle_hire_osm, dist = 20)
summary(lengths(sel) > 0)
#>    Mode   FALSE    TRUE 
#> logical     304     438

これは、ターゲットオブジェクト cycle_hire の中に、閾値距離 cycle_hire_osm 内に 438 個の点があることを示している。 それぞれの cycle_hire_osm ポイントに関連するを取得する方法は? 解答は再び st_join()だが、引数 dist を追加する(20 m 以下に設定)。

z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, dist = 20)
nrow(cycle_hire)
#> [1] 742
nrow(z)
#> [1] 762

結合結果の行数がターゲットより多いことに注意。 これは、cycle_hire の一部のサイクルレンタル・ステーションが、cycle_hire_osm で複数のマッチングを行っているためである。 重なった点の値を集約して平均値を返すには、Chapter 3 で学習した集約方法を用いると、対象と同じ行数を持つオブジェクトが得られる。

z = z |> 
  group_by(id) |> 
  summarize(capacity = mean(capacity))
nrow(z) == nrow(cycle_hire)
#> [1] TRUE

近くのステーションの収容台数は、ソース cycle_hire_osm のデータの収容台数のプロットとこの新しいオブジェクトの結果を比較することで検証できる(プロットは表示していない)。

plot(cycle_hire_osm["capacity"])
plot(z["capacity"])

この結合の結果、単純なフィーチャの属性データは空間演算で変更されたが、各フィーチャに関連するジオメトリは変更されていない。

4.2.6 空間的な集計

属性データの集約と同様に、空間データの集約では、集約された出力は集約されていない入力よりも少ない行数で済む。 平均値や合計値などの統計的な集約関数は、変数の複数の値を要約し、グループ化変数ごとに単一の値を返す。 Section 3.2.3 は、aggregate()group_by() |> summarize() が属性変数に基づいてデータを集約する方法を示したが、このセクションでは、同じ関数が空間オブジェクトでどのように機能するかを示す。 aggregate() group_by() |> summarize()

ニュージーランドの例に戻って、各地域の高所の平均的な高さを求めるとする。ターゲットオブジェクト(x または nz_height)の値がどのようにグループ化されるかを定義するのは、ソース(この場合は y または nz)のジオメトリである。 これは、base R の aggregate() メソッドを使えば、1行のコードで可能である。

nz_agg = aggregate(x = nz_height, by = nz, FUN = mean)

前のコマンドの結果は、(空間)集約オブジェクト(nz)と同じジオメトリを持つ sf オブジェクトである。これは、コマンド identical(st_geometry(nz), st_geometry(nz_agg)) で確認することができる。 先の操作の結果を Figure 4.8 に示す。これは、ニュージーランドの16の地域それぞれにおける nz_height のフィーチャの平均値を示している。 また、次のように st_join() の出力を「tidy」関数 group_by()summarize() にパイプすることでも、同じ結果を生成することができる。

ニュージーランドの各地域の上位101の高さの平均値。

FIGURE 4.8: ニュージーランドの各地域の上位101の高さの平均値。

nz_agg2 = st_join(x = nz, y = nz_height) |>
  group_by(Name) |>
  summarize(elevation = mean(elevation, na.rm = TRUE))

出来上がった nz_agg オブジェクトは、集計オブジェクト nz と同じジオメトリを持つが、関数 mean() を用いて、各地域の x の値をまとめた列が新たに追加されている。 ここでは、mean() の代わりに、median()sd() など、グループごとに単一の値を返す他の関数を使用することも可能である。 注: aggregate()group_by() |> summarize() のアプローチの違いの一つは、前者は一致しない地域名に対して NA の値を返すのに対し、後者は地域名を保持することである。 したがって、「tidy」アプローチは、集計関数や結果の列名の点でより柔軟である。 新しいジオメトリも作成する集計操作については、Section 5.2.7 で説明している。

4.2.7 不一致レイヤを結合

空間一致 (Spatial congruence) は、 空間的集計に関連する重要な概念である。 集合体(ここでは y と呼ぶ)は、2つのオブジェクトが境界を共有している場合、ターゲットオブジェクト(x)と一致している。 行政区域のデータでは、大きな単位、例えばイギリスのMiddle Layer Super Output Area(MSOAs)や他の多くのヨーロッパ諸国の地区が、多くの小さな単位で構成されていることがよくあることである。

対照的に、不一致 (incongruent) 集約オブジェクトは、ターゲットと共通の境界を共有しない (Qiu, Zhang, and Zhou 2012)。 これは、Figure 4.9 で説明されている空間集約(およびその他の空間操作)において問題となる。各サブゾーンの重心を集約すると、正確な結果を得ることができない。 面積補間は、単純な面積加重法や「ピクノフィラティック」 (pycnophylactic) 法などのより洗練されたアプローチを含む様々なアルゴリズムを使用して、1セットの面積単位から別の単位に値を転送することによってこの問題を克服している (Tobler 1979)

大きな凝集帯(半透明の青い枠)に対して、一致する面単位(左)と不一致する面単位(右)を示した図。

FIGURE 4.9: 大きな凝集帯(半透明の青い枠)に対して、一致する面単位(左)と不一致する面単位(右)を示した図。

spData パッケージには、incongruent ( Figure 4.9 の右側のパネルにある黒い縁取りのある色のついたポリゴン) と aggregating_zones ( Figure 4.9 の右側のパネルにある半透明の青い縁取りのある二つのポリゴン) という名前のデータセットが含まれている。 ここで、incongruentvalue 列が、百万ユーロ単位の地域総所得を指すと仮定しよう。 基礎となる9つの空間ポリゴンの値を、aggregating_zones の2つのポリゴンにどのように移せばいいのだろうか?

このための最も簡単で有用な方法は、面積加重空間補間で、incongruent オブジェクトから aggregating_zones の新しい列に、重なり合う面積に比例して値を転送する。入力と出力のフィーチャの空間交差が大きければ大きいほど、対応する値も大きくなる。 これは、以下のコードチャンクに示すように、st_interpolate_aw() で実装されている。

iv = incongruent["value"] # 転送する値だけを残す
agg_aw = st_interpolate_aw(iv, aggregating_zones, extensive = TRUE)
#> Warning in st_interpolate_aw.sf(iv, aggregating_zones, extensive = TRUE):
#> st_interpolate_aw assumes attributes are constant or uniform over areas of x
agg_aw$value
#> [1] 19.6 25.7

この場合、所得が小さなゾーンに均等に分布していると仮定すると、総所得はいわゆる空間的に広範な変数(面積とともに増加する)なので、集計ゾーンに入る交差点の値を合計することは意味がある(したがって、上記の警告メッセージが表示されるの)。 これは、平均所得やパーセンテージのような空間的に集中しがちな変数では異なるだろうが、面積が大きくなればなるほど増加するわけではない。 st_interpolate_aw() は、空間的に集約された変数でも同様に動作する。extensive パラメータを FALSE に設定すると、集約の際に合計関数ではなく平均を使用する。

4.2.8 距離の関係

位相関係が二値(あるフィーチャが別のフィーチャと交差するかしないか)であるのに対して、距離関係は連続的である。 2つのオブジェクト間の距離は、st_distance() 関数で計算される。 これは、ニュージーランドの最高地点と、Section 4.2.1 で作成したカンタベリー地域の地理的重心との間の距離を求めるコードである。

nz_highest = nz_height |> slice_max(n = 1, order_by = elevation)
canterbury_centroid = st_centroid(canterbury)
st_distance(nz_highest, canterbury_centroid)
#> Units: [m]
#>        [,1]
#> [1,] 115540

この結果には、2つの驚くべきことがある。

  • units、距離は10万インチではなく、10万メートルであることを教えてくれる。
  • 結果には1つの値しか含まれないが、行列として返される

この2番目の特徴は、st_distance() のもう一つの便利な特徴である、オブジェクト xy の特徴のすべての組み合わせの間の距離行列を返す機能を示唆するものである。 このコマンドは、nz_height の最初の 3 つのフィーチャと、オブジェクト co で表されるニュージーランドのオタゴおよびカンタベリー地域との間の距離を求めるものである。

co = filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)
#> Units: [m]
#>        [,1]  [,2]
#> [1,] 123537 15498
#> [2,]  94283     0
#> [3,]  93019     0

なお、nz_height の2番目と3番目のフィーチャと、co の2番目のフィーチャの距離は0である。 これは、点とポリゴンの間の距離は、ポリゴンの任意の部分への距離を指すという事実を示している。 nz_height の2番目と3番目の点は、オタゴ内にある。これは、プロットすることで確認できる(結果は示していない)。

plot(st_geometry(co)[2])
plot(st_geometry(nz_height)[2:3], add = TRUE)

4.3 ラスタデータに対する空間演算

このセクションでは、ラスタデータセットを操作するためのさまざまな基本メソッドを紹介した Section 3.3 をベースに、より高度で明示的な空間ラスタ操作を実演する。また、Section 3.3 で手動で作成したオブジェクト elevgrain を使用する。 これらのデータセットは、読者の便宜を図るため spData パッケージにも含まれている。

4.3.1 空間部分集合

前の章(Section 3.3)では、特定のセル ID や行と列の組み合わせに関連する値を取得する方法を紹介した。 ラスタオブジェクトは、位置(座標)などの空間オブジェクトを抽出することも可能である。 部分集合に座標を使用するには、terra の関数 cellFromXY() で座標をセル ID に「変換」することができる。 別の方法として、terra::extract()(注意: tidyverse の中にも extract() という関数がある)を使って値を抽出することができる。 以下に、座標 0.1, 0.1 に位置する点を覆うセルの値を求める方法を示す。

id = cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
# the same as
terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))

ラスタオブジェクトは、以下のコードのように、別のラスタオブジェクトの内部に部分集合することもできる。

clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip]
# we can also use extract
# terra::extract(elev, ext(clip))

これは、Figure 4.10 に示すように、2 番目のラスタ(ここでは clip)の範囲内にある最初のラスタオブジェクト(この場合は elev)の値を取得することになる。

元のラスタ(左)。ラスタマスク(中)。ラスタをマスクした出力(右)。

FIGURE 4.10: 元のラスタ(左)。ラスタマスク(中)。ラスタをマスクした出力(右)。

上記の例では、特定のセルの値を返したが、多くの場合、ラスタデータセットの部分集合操作による空間出力が必要である。 これは、[ 演算子で、drop = FALSE を使って行うことができる。 以下のコードは、elev の先頭行の2つのセルをラスタオブジェクトとして返す(出力の最初の2行のみ表示)。

elev[1:2, drop = FALSE]    # spatial subsetting with cell IDs
#> class       : SpatRaster 
#> dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#> ...

空間部分集合のもう一つの一般的な使用例は、logical(または NA)値のラスタを使用して、同じ範囲と解像度の別のラスタをマスクする場合である(Figure 4.10 に図示)。 この場合、[mask() 関数を使用することができる(結果は示していない)。

# ラスタのマスクを作成
rmask = elev
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)

上記のコードでは、rmask というマスクオブジェクトを作成し、NATRUE にランダムな値を割り当てている。 次に、elev のうち、TRUE となる値を rmask に保持したい。 つまり、elevrmask でマスクしたい。

# 空間的に部分集合を作成
elev[rmask, drop = FALSE]     # [ operator を使用
mask(elev, rmask)             # mask() を使用

上記の方法は、一部の値(例えば、間違っていると予想される値)を NA に置き換えるために使用することも可能である。

elev[elev < 20] = NA

これらの操作は、実際のところ、2つのラスタをセル単位で比較するローカルの論理操作である。 次のサブセクションでは、これらの操作と関連する操作についてより詳しく説明する。

4.3.2 マップ代数

「マップ代数」(Map algebra) という用語は、1970年代後半に、地理的なラスタデータおよび(あまり目立たないが)ベクタデータを分析するための「規則、機能、および技術のセット」を表すために作られたものである (Tomlin 1994) ここでは、マップ代数をより狭く定義し、周囲のセル、ゾーン、またはすべてのセルに適用される統計関数を参照して、ラスタセル値を修正または要約する操作としている。

ラスタデータセットは暗黙的に座標を保存しているだけなので、マップ代数演算は高速になる傾向があり、そのため古いことわざでは「ラスタは高速だがベクタは補正が効く」とされている。 ラスタデータセットのセルの位置は、その行列の位置と、データセットの解像度および原点(ヘッダに格納)を使用して計算することができる。 しかし、処理にあたっては、処理後のセル位置が変わらないことを確認すれば、セルの地理的な位置はほとんど関係ない。 さらに、2つ以上のラスタデータセットが同じ範囲、投影、解像度を共有している場合、それらを行列として処理することができる。

これは、マップ代数が terra パッケージで動作する方法である。 まず、ラスタデータセットのヘッダを照会し、(マップ代数演算が複数のデータセットに対して行われる場合)データセットの互換性を確認する。 第二に、マップ代数はいわゆる一対一の位置対応を保持しており、セルは移動できないことを意味している。 これは、行列の掛け算や割り算などで値の位置が変わる行列代数とは異なる。

マップ代数(またはラスタデータによる地図作成)では、ラスタ操作を4つのサブクラスに分け (Tomlin 1990)、それぞれが1つまたは複数のグリッドを同時に処理するようにしている。

  1. ローカル (Local)、つまりセル単位の操作
  2. フォーカル (Focal)、つまり近傍(Nighborhood)オペレーション。 多くの場合、出力セルの値は、3×3の入力セルブロックの結果
  3. ゾーン (Zonal) 演算は、フォーカル演算と似ているが、新しい値を計算する周囲の画素グリッドは不規則なサイズと形状を持つことができる。
  4. グローバル (Global) またはラスタ単位の操作。 つまり、出力セルは1つまたは複数のラスタ全体から潜在的にその値を引き出す。

このトポロジーは、マップ代数演算を、各ピクセル処理ステップに使用するセル数と出力の種類によって分類したものである。 なお、ラスタ演算は、地形、水文解析、画像分類などの分野ごとの分類方法もある。 以下では、各タイプのマップ代数演算の使用方法について、動作例を参照しながら説明する。

4.3.3 ローカルオペレーション

ローカル操作は、1つまたは複数のレイヤにおけるすべてのセル単位の操作で構成されている。 ラスタ代数とは、ローカル操作の典型的な使用例で、ラスタからの値の加算や減算、ラスタの二乗や乗算が含まれる。 ラスタ代数では、特定の値(下の例では5)より大きいラスタセルをすべて見つけるなどの論理演算も可能である。 terra パッケージは、以下のように、これらすべての操作に対応している( Figure 4.11 )。

elev + elev
elev^2
log(elev)
elev > 5
elev ラスタオブジェクトのさまざまなローカル操作の例:2つのラスタの加算、二乗、対数変換の適用、論理演算の実行。

FIGURE 4.11: elev ラスタオブジェクトのさまざまなローカル操作の例:2つのラスタの加算、二乗、対数変換の適用、論理演算の実行。

ローカル演算のもう一つの良い例は、デジタル標高モデルを低標高(クラス1)、中標高(クラス2)、高標高(クラス3)にグループ化するように、数値の間隔をグループに分類することである。 classify() コマンドを使って、まず再分類行列を作る必要がある。ここで、最初の列はクラスの下限、2番目の列は上限に対応する。 3列目は、1列目と2列目で指定した範囲の新しい値を表している。

rcl = matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
rcl
#>      [,1] [,2] [,3]
#> [1,]    0   12    1
#> [2,]   12   24    2
#> [3,]   24   36    3

ここでは、0~12、12~24、24~36の範囲のラスタ値をそれぞれ1、2、3に再分類している。

recl = classify(elev, rcl = rcl)

classify() 関数は、カテゴリー別ラスタのクラス数を減らしたい場合にも使用できる。 Chapter 14 では、いくつかの追加的な再分類を実施する予定である。

算術演算子とは別に、app()tapp()lapp() 関数も使用することができる。 より効率的であるため、大規模なラスタデータが存在する場合に適している。 さらに、出力ファイルを直接保存することも可能である。 関数 app() は、ラスタの各セルに関数を適用し、複数のレイヤの値を1つのレイヤにまとめる(合計を計算するなど)ために使用される。 tapp() は、app() の拡張で、ある操作を行いたいレイヤの部分集合(index の引数を参照)を選択することができるようになっている。 最後に、関数 lapp() は、レイヤを引数として各セルに関数を適用することができる。lapp() のアプリケーションを以下に示す。

正規化差分植生指数(normalized difference vegetation index, NDVI)の算出は、よく知られたローカル(ピクセル単位)のラスタ処理である。 正の値は生きた植物が存在することを示す(ほとんどが0.2以上)。 NDVI は、Landsat や Sentinel などの衛星システムから得られるリモートセンシング画像の赤色および近赤外(near-infrared、NIR)バンドから算出されるものである。 植物は可視光線、特に赤色光を強く吸収し、近赤外光を反射するため、NVDI の計算式が成り立っている。

\[ \begin{split} NDVI&= \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}\\ \end{split} \]

ザイオン国立公園のマルチスペクトル衛星ファイルについて、NDVI を計算してみよう。

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)

ラスタオブジェクトは、青、緑、赤、近赤外(NIR)の4つの衛星バンドを持っている。 次のステップは、NDVI の計算式を R の関数に実装することである。

ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}

この関数は、2つの数値引数(nirred)を受け取り、NDVI 値を含む数値ベクタを返す。 lapp()fun 引数として使用することができる。 この関数が必要とするのは2つのバンド(元のラスタの4つではない)であり、それらはNIR、赤の順である必要があることを覚えておく必要がある。 そのため、入力ラスタを部分集合し、計算を行う前に multi_rast [[c(4, 3)] ] で部分集合してから計算を行う。

ndvi_rast = lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)

その結果を右図( Figure 4.12 )に示すように、同じ領域のRGB画像(同図の左図)と比較することができる。 これにより、NDVI値が最も大きいのは北部の密林地帯、最も低いのは北部の湖と雪山の尾根に関連していることがわかる。

ザイオン国立公園の衛星ファイルの例で計算されたRGB画像(左)とNDVI値(右)

FIGURE 4.12: ザイオン国立公園の衛星ファイルの例で計算されたRGB画像(左)とNDVI値(右)

予測マッピングも、ローカルラスタ操作の興味深い応用例である。 応答変数は、例えば、種の豊富さ、地滑りの存在、木の病気、作物の収穫量など、空間における測定または観測された点に対応する。 その結果、様々なラスタ(標高、pH、降水量、気温、土地被覆、土壌等級など)から、宇宙や空中の予測変数を簡単に取得することができるようになる。 その後、lm()glm()gam() または機械学習技術を使用して、予測因子の関数として応答をモデル化する。 したがって、ラスタオブジェクトの空間予測は、予測ラスタ値に推定係数を適用し、出力ラスタ値を合計することで行うことができる(Chapter 15 を参照)。

4.3.4 フォーカルオペレーション

ローカルな機能とは、1つのセル(複数の層からなる可能性もある)を対象とするものであるのに対し、フォーカルな機能とは、中心(焦点)のセルとその近隣のセルを考慮したものである。 近傍領域(カーネル、フィルタ、移動窓とも呼ばれる)は、通常、3×3セル(つまり、中心セルとその周囲8個の近傍セル)のサイズであるが、ユーザーが定義する他の形状(必ずしも長方形ではない)をとることができる。 フォーカル操作は、指定された近傍領域内のすべてのセルに集約関数を適用し、対応する出力を中心セルの新しい値として使用し、次の中心セル(Figure 4.13)へと進む。 この操作は、空間フィルタリング (spatial filtering) やコンボリューション (convolution) などと呼ばれることもある (Burrough, McDonnell, and Lloyd 2015)

R では、focal() 関数で空間フィルタリングを行うことができる。 移動窓の形状を matrix で定義する。その値は重みに対応する(以下のコードチャンクの w パラメータを参照)。 次に、fun パラメータで、この近傍領域に適用したい関数を指定することができる。 ここでは、最小値を選んでいるが、sum()mean()var() など、他の要約関数も使用可能である。

r_focal = focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)

この関数は、例えば、プロセス中の NA を削除すべきか(na.rm = TRUE)、しないか(na.rm = FALSE)などの追加引数も受け付ける。

焦点演算による入力ラスタ(左)と出力ラスタ(右)-3×3の移動窓で最小値を求める。

FIGURE 4.13: 焦点演算による入力ラスタ(左)と出力ラスタ(右)-3×3の移動窓で最小値を求める。

期待通りの出力が得られるかどうか、すぐに確認することができる。 この例では、最小値は常に移動ウィンドウの左上隅でなければならない(入力ラスタは、左上隅から始まるセルの値を行単位で1つずつ増加させることによって作成したことを思い出してほしい)。 この例では、重み付け行列は1だけで構成されており、各セルの出力に対する重みが同じであることを意味しているが、これは変更可能である。

画像処理では、焦点関数やフィルターが重要な役割を担っている。 ローパスフィルタやスムージングフィルタは、平均関数を用いて極端な部分を除去する。 カテゴリーデータの場合、平均値を最頻値に置き換えることができる。 それに対して、ハイパスフィルターはフィーチャを強調する。 ここでは、ライン検出のラプラスフィルターやソーベルフィルターがその例として挙げられるだろう。 R での使い方は、focal() のヘルプページで確認してみよう(この章の最後の演習でも使用する)。

傾斜、アスペクト、流れ方向などの地形特性を計算する地形処理では、焦点関数に依存している。 terrain() は、これらの指標の計算に使用することができる。ただし、勾配を計算する Zevenbergen and Thorne 法を含むいくつかの地形アルゴリズムは、この terra 関数には実装されていない。 その他、曲率、寄与率、湿潤指数など多くのアルゴリズムが、オープンソースのデスクトップ型地理情報システム(GIS)ソフトウェアに実装されている。 Chapter 10 は、このような GIS 機能を R 内からアクセスする方法を示している。

4.3.5 ゾーンオペレーション

ゾーン演算は、フォーカル演算と同様に、複数のラスタセルに集計関数を適用する。 しかし、前節で紹介したフォーカル操作の場合の事前定義された近傍ウィンドウとは対照的に、ゾーン操作の場合は、通常カテゴリー値で構成される第二ラスタがゾーンフィルター(または「ゾーン」)を定義する。 そのため、ゾーンフィルタを定義するラスタセルは、必ずしも隣接している必要はない。 粒径ラスタはその良い例で、Figure 3.2 の右側のパネルに示されているように、異なる粒径がラスタ全体に不規則に広がっていることがわかる。 最後に、ゾーン操作の結果は、ゾーンごとにグループ化された要約表となる。このため、この操作は、GIS の世界ではゾーン統計とも呼ばれる。 これは、ラスタオブジェクトを返すフォーカルオペレーションとは対照的である。

次のコードチャンクは、zonal() 関数を使用して、例えば、各粒度クラスに関連する平均標高を計算する。

z = zonal(elev, grain, fun = "mean")
z
#>   grain elev
#> 1  clay 14.8
#> 2  silt 21.2
#> 3  sand 18.7

これは、各カテゴリーの統計値 、ここでは各粒度クラスの平均高度を返す。 注: as.raster 引数を TRUE に設定することで、各ゾーンの統計情報を計算したラスタを取得することも可能である。

4.3.6 グローバルな操作と距離

グローバル操作は、ラスタデータセット全体が1つのゾーンに相当するゾーン操作の特殊なケースである。 最も一般的なグローバル操作は、最小値や最大値など、ラスタデータセット全体の記述統計である。これらの操作については、Section 3.3.2 ですでに説明した。

それ以外にも、距離や重みのラスタの計算にもグローバル操作は有効である。 最初のケースでは、各セルから特定のターゲットセルまでの距離を計算することができる。 例えば、最も近い海岸までの距離を計算したい場合がある(terra::distance() も参照)。 また、地形も考慮したい。つまり、純粋な距離だけでなく、海岸に行くときに山脈を越えないようにしたいのである。 そのためには、標高が1メートル増えるごとにユークリッド距離が「伸びる」ように、標高で距離に重みをつければいいのだ。 Visibility と viewshed の計算もグローバル操作に属する ( Chapter 10 の演習では、viewshed ラスタを計算する)。

4.3.7 ベクタ処理における写像代数の対応

多くのマップ代数演算はベクタ処理に対応するものである (Liu and Mason 2009)。 最大距離のみを考慮した距離ラスタの計算(グローバル演算)は、ベクタバッファ演算(Section 5.2.5)と同等である。 ラスタデータの再分類(入力に応じてローカルまたはゾーン関数)は、ベクタデータの溶解(Section 4.2.4)と同等である。 2つのラスタを重ね合わせ(ローカル操作)、一方がマスクを表す NULL または NA の値を含む場合、ベクタクリッピング(Section 5.2.5)に類似している。 空間クリッピングとよく似たものに、2つのレイヤを交差させるものがある( Section 4.2.1 )。 違いは、これら2つのレイヤ(ベクタまたはラスタ)は、単に重複する領域を共有することである(例として Figure 5.8 を参照)。 ただし、表現には注意が必要である。 ラスタデータモデルとベクタデータモデルでは、同じ言葉でも微妙に意味が異なることがある。 集計とは、ベクタデータの場合はポリゴンを分解することであり、ラスタデータの場合は解像度を上げることである。 実際、ポリゴンを分解したり集約したりすると、解像度が下がるという見方もできる。 しかし、ゾーン演算は、セル解像度の変更と比較して、より優れたラスタ同等物である可能性がある。 ゾーン演算は、あるラスタのセルを、別のラスタのゾーン(カテゴリ)に合わせて、集約関数(上記参照)を使って分解することができる。

4.3.8 ラスタのマージ

NDVI を計算し(Section 4.3.3 参照)、さらに、調査地域内の観測の標高データから地形属性を計算したいとする。 このような計算には、リモートセンシングの情報が必要である。 リモートセンシング画像は、特定の空間範囲をカバーするシーンに分割されることが多く、1つの調査エリアが複数のシーンにまたがっていることがよくある。 そして、調査対象シーンをマージする必要がある。 一番簡単なのは、これらのシーンをマージする、つまり並べることである。 これは例えば、デジタル標高データ(SRTM、ASTER)で可能である。 以下のコードでは、まずオーストリアとスイスの SRTM 標高データをダウンロードする(国番号については、geodata 関数 country_codes() を参照)。 第二段階では、2つのラスタを1つに統合する。

aut = geodata::elevation_30s(country = "AUT", path = tempdir())
ch = geodata::elevation_30s(country = "CHE", path = tempdir())
aut_ch = merge(aut, ch)

terramerge() コマンドは、2つの画像を合成し、重なった場合は、最初のラスタの値を使用する。

重複する値が互いに対応しない場合、マージアプローチはほとんど意味がない。 このようなケースは、撮影日が異なるシーンの分光画像を合成する場合によくある。 merge() コマンドはそのまま使え、出来上がった画像にはっきりとした枠が表示される。 一方、mosaic() コマンドでは、オーバーラップする領域に対して関数を定義することができる。 例えば、平均値を計算することもできる。この場合、マージされた結果における明確な境界線は滑らかになるだろうが、ほとんどの場合、それを消すことはできない。 R によるリモートセンシングの詳しい紹介は、Wegmann, Leutner, and Dech (2016) を参照。

4.4 演習

E1. It was established in Section 4.2 that Canterbury was the region of New Zealand containing most of the 100 highest points in the country. How many of these high points does the Canterbury region contain?

Bonus: plot the result using the plot() function to show all of New Zealand, canterbury region highlighted in yellow, high points in Canterbury represented by red crosses (hint: pch = 7) and high points in other parts of New Zealand represented by blue circles. See the help page ?points for details with an illustration of different pch values.

E2. Which region has the second highest number of nz_height points, and how many does it have?

E3. Generalizing the question to all regions: how many of New Zealand’s 16 regions contain points which belong to the top 100 highest points in the country? Which regions?

  • Bonus: create a table listing these regions in order of the number of points and their name.

E4. Test your knowledge of spatial predicates by finding out and plotting how US states relate to each other and other spatial objects.

The starting point of this exercise is to create an object representing Colorado state in the USA. Do this with the command colorado = us_states[us_states$NAME == "Colorado",] (base R) or with with the filter() function (tidyverse) and plot the resulting object in the context of US states.

  • Create a new object representing all the states that geographically intersect with Colorado and plot the result (hint: the most concise way to do this is with the subsetting method [).
  • Create another object representing all the objects that touch (have a shared boundary with) Colorado and plot the result (hint: remember you can use the argument op = st_intersects and other spatial relations during spatial subsetting operations in base R).
  • Bonus: create a straight line from the centroid of the District of Columbia near the East coast to the centroid of California near the West coast of the USA (hint: functions st_centroid(), st_union() and st_cast() described in Chapter 5 may help) and identify which states this long East-West line crosses.

E5. Use dem = rast(system.file("raster/dem.tif", package = "spDataLarge")), and reclassify the elevation in three classes: low (<300), medium and high (>500). Secondly, read the NDVI raster (ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))) and compute the mean NDVI and the mean elevation for each altitudinal class.

E6. Apply a line detection filter to rast(system.file("ex/logo.tif", package = "terra")). Plot the result. Hint: Read ?terra::focal().

E7. Calculate the Normalized Difference Water Index (NDWI; (green - nir)/(green + nir)) of a Landsat image. Use the Landsat image provided by the spDataLarge package (system.file("raster/landsat.tif", package = "spDataLarge")). Also, calculate a correlation between NDVI and NDWI for this area (hint: you can use the layerCor() function).

E8. A StackOverflow post shows how to compute distances to the nearest coastline using raster::distance(). Try to do something similar but with terra::distance(): retrieve a digital elevation model of Spain, and compute a raster which represents distances to the coast across the country (hint: use geodata::elevation_30s()). Convert the resulting distances from meters to kilometers. Note: it may be wise to increase the cell size of the input raster to reduce compute time during this operation (aggregate()).

E9. Try to modify the approach used in the above exercise by weighting the distance raster with the elevation raster; every 100 altitudinal meters should increase the distance to the coast by 10 km. Next, compute and visualize the difference between the raster created using the Euclidean distance (E7) and the raster weighted by elevation.