10 GISソフトウェアへのブリッジ

必須パッケージ

  • 本章では、QGIS 、SAGA 、GRASS がインストールされていること、および以下のパッケージが添付されていることを条件とする。
# remotes::install_github("paleolimbot/qgisprocess")
library(qgisprocess)
library(Rsagacmd)
library(rgrass)
library(rstac)
library(gdalcubes)

10.1 イントロダクション

R の特徴は、その使い方にある。 コマンドを入力し、Enter(RStudio のソースエディタでコードを書いている場合は Ctrl+Enter )を押すと、対話的にコマンドが実行される。 このようなコンピュータとの対話方法は、コマンドラインインターフェイス(CLI) と呼ばれる(下記の注の定義参照)。 CLI は R だけのものではない。46 一方、専用の GIS パッケージでは、グラフィカルユーザーインターフェイス(Grafical User Interface, GUI\index{graphical user interface)に重点が置かれる傾向がある。 QGIS、SAGA、GRASS、gvSIG をシステムターミナルや組み込み CLI から操作することは可能であるが、「マウス操作」が一般的である。 QGIS (Sherman 2008) の作成者である Gary Sherman によれば、これは多くの GIS ユーザーがコマンドラインの利点を見逃していることを意味するとのことである。

「近代的」GIS ソフトウェアの発展に伴い、マウス操作を好む人がほとんどである。それはいいことであるが、とてつもない量のコマンドラインには、柔軟性とパワーが待っている。繰り返しコマンドラインで作業をする場合、同じことを GUI で行う場合よりも短い時間で行うことができる。

「CLI vs GUI」 の議論は逆境になることもあるが、その必要はない。どちらの選択肢も、目の前のタスクとユーザーのスキルセットに応じて、互換的に使用することができる47 。 R が提供するような優れた CLI(そして RStudio のような IDE によって強化されている)の利点は数多くある。 良い CLI とは、

  • 繰り返の作業の自動化を促進する
  • 優れた科学的実践とデータサイエンスのバックボーンである透明性と再現性を可能にする
  • 既存の機能を修正したり、新しい機能を実装するためのツールを提供することで、ソフトウェア開発を促進する
  • 多くの分野や産業で求められている、将来性のあるプログラミングスキルを身につけることができる
  • ユーザーフレンドリーで高速、効率的なワークフローを実現する

一方、GUIベースのGIS システム(特にQGIS )も有利である。 良い GIS GUI とは、

  • 学習曲線が「浅い」ので、新しい言語を何時間も学ぶことなく、地理データを探索し、可視化することができる
  • トレース、スナップ、トポロジーツールなど、「デジタイジング」(新しいベクタデータセットの作成)のための優れたサポートを提供する48
  • 地上基準点によるジオリファレンス(ラスタ画像と既存地図とのマッチング)、オルソレクションが可能である
  • 立体視マッピング(LiDAR、Structure from Motionなど)に対応する
  • オブジェクト指向のリレーショナルデータモデル、トポロジー、高速な(空間)クエリを備えた空間データベース管理システムへのアクセスを提供する

専用 GIS のもう一つの利点は、何百もの「ジオアルゴリズム」(地理的問題を解決するための計算レシピ、Chapter 11 を参照)を利用できることである。 これらの多くは、本章のテーマであり動機でもある「GIS ブリッジ」を経由する以外、R のコマンドラインから利用することはできない。49

コマンドラインインタフェースとは、ユーザーが連続したテキスト行(コマンドライン)を介してコマンドを発行することで、コンピュータプログラムと対話するための手段である。 Linuxの bash や Windows の PowerShell が一般的な例である。 CLI は、RStudio for R のような IDE で補強することができ、コードの自動補完やユーザーエクスペリエンスを向上させる機能を提供する。

R はインターフェース言語として誕生した。 その前身である S は、他の言語(特に FORTRAN)の統計アルゴリズムへのアクセスを提供していたが、直感的な 入力・評価・出力(read-evaluate-print loop, REPL)でアクセスできるようになっていた (Chambers 2016)。 Chapter 1 で述べたように、R はこの伝統を受け継ぎ、特に C++ などの多くの言語へのインタフェースを提供している。 R は GIS として設計されたものではない。 しかし、専用の GIS とのインターフェイスが可能なため、驚異的な地理空間能力を発揮する。 R は統計プログラミング言語としてよく知られているが、GIS のワークフローを再現することができ、さらに(比較的)一貫した CLI という利点があることを知らない人が多いようである。 さらに、R は、インタラクティブ/アニメーションマップ作成(Chapter 9 参照)や空間統計モデリング( Chapter 12 参照)など、ジオコンピューティングのいくつかの分野では GIS を凌駕する性能を有している。 この章では、3つの成熟したオープンソース GIS 製品への「ブリッジ」(Table 10.1 参照)に焦点を当てる。QGIS (qgisprocess; Section 10.2 )、SAGA (Rsagacmd; Section 10.3)、GRASS (rgrass; Section 10.4)である。50

TABLE 10.1: Comparison between three open-source GIS. Hybrid refers to the support of vector and raster operations.
GIS First release No. functions Support
QGIS 2002 >1000 hybrid
SAGA 2004 >600 hybrid
GRASS 1982 >500 hybrid

また、R-GIS ブリッジを補完するために、第 2 部では空間ライブラリへのインタフェース(Section 10.6.1、 空間データベース (Section 10.6.2)、 地球観測データのクラウド処理 (Section 10.7) について簡単に紹介する。)

10.2 qgisprocess から QGIS

QGIS は、最も人気のあるオープンソース GIS の 1 つである [Table 10.1 ; Graser and Olaya (2015)]。 その最大の特徴は、他のオープンソース GIS と統一されたインターフェースを提供する点にある。 つまり、QGIS から、GDAL、GRASS、SAGA を利用することができるのである (Graser and Olaya 2015)。 バージョン 3.14 以降、QGIS はコマンドライン API qgis_process を提供しており、QGIS GUI の外でこれらのすべてのジオアルゴリズム(セットアップによっては 1000 以上)を実行することが可能である。

qgisprocess パッケージは、QGIS コマンドラインユーティリティをラップ(訳注:プログラミング用語で、機能を使いやすく覆い隠す wrap こと)しているため、R セッションから QGIS、GDAL、GRASS、SAGA アルゴリズムを呼び出すことが可能である。 qgisprocess を実行する前に、QGIS と、SAGA や GRASS などの(サードパーティの)依存パッケージがすべてインストールされていることを確認してみよう。

library(qgisprocess)
#> Using 'qgis_process' in the system PATH.
#> QGIS version: 3.26.1-Buenos Aires
#> ...

このパッケージは、QGISのインストールを自動的に検出しようとし、検出できない場合は警告を発する。51 設定に失敗した場合の解決策としては、options(qgisprocess.path = "path/to/your_qgis_process")、環境変数 R_ QGISPROCESS_PATH を設定する方法が考えられる。 上記の方法は、複数のQGISがインストールされており、どれを使うかを決めたい場合にも使える。

次に、どのプロバイダ(異なるソフトウェアを意味する)が自分のコンピュータで利用できるかを調べる。

qgis_providers()
#> # A tibble: 6 × 2
#>   provider provider_title   
#>   <chr>    <chr>            
#> 1 3d       QGIS (3D)       
#> 2 gdal     GDAL             
#> 3 grass7   GRASS            
#> 4 native   QGIS (native c++)
#> 5 qgis     QGIS             
#> 6 saga     SAGA

出力表から、QGISのジオアルゴリズム(native, qgis, 3d)と、サードパーティプロバイダの GDAL、SAGA、GRASS の外部アルゴリズムを QGIS インターフェースを通して使用できることが確認できた。

これで、R から QGIS の地理計算をする準備ができた。 それでは、2つの事例を試してみよう。 最初のものは、異なる境界線を持つ2つのポリゴンデータセットを統合する方法を示している (Section 10.2.1 ) 。 もう一つは、ラスタ(Section 10.2.2 )で表現された数値標高モデルから新しい情報を導き出すことに重点を置いている。

10.2.1 ベクタデータ

異なる空間単位(地域、行政単位など)を持つ2つのポリゴンオブジェクトがある場合を考えてみよう。 この2つのオブジェクトを統合して、すべての境界線と関連する属性を含む1つのオブジェクトにすることが目標である。 Section 4.2.7 ( Figure 10.1 ) ですでに出会った不整合なポリゴンを再び利用する。 どちらのポリゴンデータセットも spData パッケージで提供されており、その両方に地理的な CRS を使用したい(Chapter 7 も参照)。

data("incongruent", "aggregating_zones", package = "spData")
incongr_wgs = st_transform(incongruent, "EPSG:4326")
aggzone_wgs = st_transform(aggregating_zones, "EPSG:4326")
Illustration of two areal units: incongruent (black lines) and aggregating zones (red borders).

FIGURE 10.1: Illustration of two areal units: incongruent (black lines) and aggregating zones (red borders).

この作業を行うアルゴリズムを見つけるには、qgis_algorithms() 関数の出力を検索すればよい。 この関数は、利用可能なすべてのプロバイダと、それらが含むアルゴリズムを含むデータフレームを返す。52

qgis_algo = qgis_algorithms()

qgis_algo オブジェクトには多くの列があるが、通常、私たちが関心を持つのは、プロバイダとアルゴリズム名の情報を組み合わせた algorithm 列だけである。 関数の短い説明文に「union」という単語が含まれていると仮定すると 、以下のコードを実行して、興味のあるアルゴリズムを見つけることができる。

grep("union", qgis_algo$algorithm, value = TRUE)
#> [1] "native:multiunion" "native:union"      "saga:fuzzyunionor" "saga:polygonunion"

上記のリストにあるアルゴリズムの1つ、"native:union" は、探している機能の可能性が高そうである。 次のステップは、このアルゴリズムが何をするのか、どう使えばいいのかを知ることである。 qgis_show_help() がそのためにある。これは、アルゴリズムが何をし、引数や出力についての要約を返す。53 これによって、出力が長くなる。 ここでは、"native:union" という引数で得られる結果を示す。

alg = "native:union"
qgis_arguments(alg) |> 
  dplyr::mutate(acceptable_values = unlist(acceptable_values)) |>
  dplyr::select(name, description, acceptable_values)

INPUTOVERLAYOVERLAY_FIELDS_PREFIXOUTPUTがあることがすぐに分かる。 これらの引数は、ベクタレイヤへのパスを期待しているようである(列名: acceptaple_values)。 ただし、qgisprocess パッケージでは、このような場合にも sf オブジェクトを提供することが可能となる。54 この機能は便利だが、もしすでに R セッション中に空間データがある場合、qgisprocess にそのパスを送る時のみ使用することを勧める。なぜなら、qgisprocess がジオアルゴリズムを実行する際に最初にすることは、R セッション中の空間データを QGIS が分かる .gpkg または .tif に変換することだからである。

最後に、QGIS に作業を任せることができる。 qgisprocess の主な機能は qgis_run_algorithm() である。 使用するアルゴリズム名とヘルプに示される名前付き引数のセットを受け取り、期待される計算を実行する。 今回のケースでは、INPUTOVERLAYOUTPUT の3つの引数が重要だと思われる。 最初の INPUT は、主なベクタオブジェクト incongr_wgs であり、2番目の OVERLAY は、aggzone_wgs である。 最後の引数、OUTPUT は新しいベクタファイルへのパスを期待する。しかし、パスを提供しない場合、qgisprocessは自動的に一時ファイルを作成する。 引数 .quiet は、qgisprocess のメッセージを減らす。

union = qgis_run_algorithm(alg, INPUT = incongr_wgs, OVERLAY = aggzone_wgs, 
                           .quiet = TRUE)
union

上記のコードを実行すると、2つの入力オブジェクトが一時的な .gpkg ファイルに保存され、選択されたアルゴリズムがそれらに実行され、一時的な .gpkg ファイルが出力として返される。 qgisprocess パッケージは、qgis_run_algorithm() の結果を、この場合は出力ファイルへのパスを含むリストとして保存する。 このファイルを R に読み戻すには、read_sf() (例, union_sf = read_sf(union [[1] ]) を使うか、st_as_sf() を使って直接読み込むことができる。

union_sf = st_as_sf(union)

QGIS の union の操作は、交差 と2つの入力レイヤーの対称差分を用いて、2つの入力レイヤーを1つのレイヤーにマージすることに注意(ちなみに、これは GRASS と SAGA でユニオン操作をするときのデフォルトでもある)。 これは st_union(incongr_wgs, aggzone_wgs) (演習参照)とは違う

その結果、union_sf、2つの入力オブジェクトよりも多くのフィーチャを持つポリゴンとなる。 しかし、これらのポリゴンの多くは小さく、実際の領域を表しているわけではなく、2つのデータセットの詳細度が異なるために生じたものであることに注意しておこう。 これらの誤差はスライバーポリゴンと呼ばれる( Figure 10.2 の左側のパネルにある赤い色のポリゴンを参照)。 スライバーを識別する一つの方法として、面積が比較的非常に小さいポリゴン、ここでは例えば25000m2を見つけ、次にそれを削除する。 適切なアルゴリズムを探そう。

grep("clean", qgis_algo$algorithm, value = TRUE)

今回見つかったアルゴリズム( v.clean )は、QGIS ではなく、GRASS GIS に含まれている。 GRASS GIS の v.clean は、空間ベクタデータのトポロジーをクリーニングする強力なツールである。 重要なのは、qgisprocessを通して使用できることである。

前のステップと同様に、このアルゴリズムのヘルプを見るところから始めよう。

qgis_show_help("grass7:v.clean")

ここでは出力を省略した。実際のヘルプテキストはかなり長く、多くの引数を含んでいる。55 これは、v.clean がマルチツールであり、さまざまな種類のジオメトリをクリーニングし、さまざまな種類のトポロジー問題を解決することができることがある。 この例では、いくつかの引数に絞って説明するが、v.clean の機能については、 このアルゴリズムのドキュメントを勧める。

qgis_arguments("grass7:v.clean") |>
  dplyr::select(name, description) |>
  dplyr::slice_head(n = 4)
#> # A tibble: 4 × 2
#>   name      description                              
#>   <chr>     <chr>                                    
#> 1 input     Layer to clean                           
#> 2 type      Input feature type                       
#> 3 tool      Cleaning tool                            
#> 4 threshold Threshold (comma separated for each tool)

このアルゴリズムの主な引数は input で、これは私たちのベクタオブジェクトである。 次に tool の選択であるが、これはクリーニングの方法である。56 v.clean には、重複した形状の削除、線間の微小角度の削除、微小領域の削除など、12 種類のツールが存在する。 今回は、後者のツール、rmarea を解説する。 いくつかのツール(rmarea を含む)は、追加の引数 threshold を必要とし、その動作は選択されたツールに依存する。 この場合、rmarea ツールは、threshold で与えられた値より小さいか等しい領域をすべて削除する。

このアルゴリズムを実行し、その出力を新しい sf オブジェクト clean_sf に変換してみよう。

clean = qgis_run_algorithm("grass7:v.clean", input = union_sf,
                           tool = "rmarea", threshold = 25000, .quiet = TRUE)
clean_sf = st_as_sf(clean)

その結果、(fig:sliver) の右側のパネルでは、予想通り、灰色のポリゴンが削除されているように見える。

灰色部分を赤で強調(左)と灰色部分を除去(右)。

FIGURE 10.2: 灰色部分を赤で強調(左)と灰色部分を除去(右)。

10.2.2 ラスタデータ

デジタル標高モデル(DEM)には、ラスタセルごとの標高情報が含まれている。 衛星航法、水流モデル、表面分析、可視化など、さまざまな用途で使用されている。 ここでは、DEM ラスタから統計学習における予測因子として利用可能な新しい情報を導き出してみたい。 例えば、様々な地形パラメータは、地滑りの予測に役立つ(Chapter 12 参照)。

このセクションでは、dem.tif を使用することにする。これは、 Mongón 調査地域のデジタル標高モデルである(Land Process Distributed Active Archive Center からダウンロード、?dem.tif も参照)。 解像度は約 30 m × 30 mで、投影型 CRS を使用している。

library(qgisprocess)
library(terra)
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))

terra パッケージの terrain() では、傾斜、アスペクト、TPI (Topographic Position Index) 、TRI (Topographic Ruggedness Index) 、粗さ、流れ方向など、地形の基本特性を算出することができる。 とはいえ、GIS は、地形の特性に関する機能が他にもたくさんあり、文脈によってはより適しているものもある。 例えば、地形湿潤指数(TWI)は、水文・生物学的プロセスの研究に有用であることがわかっている (Sørensen, Zinko, and Seibert 2006)。 このインデックスのアルゴリズムリストを、"wetness" というキーワードで検索してみよう。

qgis_algo = qgis_algorithms()
grep("wetness", qgis_algo$algorithm, value = TRUE)
#> [1] "saga:sagawetnessindex"           "saga:topographicwetnessindextwi"

上記のコードの出力は、目的のアルゴリズムが SAGA GIS ソフトウェアに存在することを示唆するものである。57 SAGA はハイブリッド GIS であるが、主にラスタ処理、ここでは特にデジタル標高モデル (土壌特性、地形属性、気候パラメータ)に重点を置いている。 したがって、SAGA は大規模な(高解像度の)ラスタデータセットの高速処理に特に優れている (Conrad et al. 2015)

"saga:sagawetnessindex" アルゴリズムは、実際には修正された TWI であり、谷底に位置するセルに対してより現実的な土壌水分ポテンシャルをもたらすものである (Böhner and Selige 2006)

qgis_show_help("saga:sagawetnessindex")
# 出力は非表示

ここでは、デフォルトの引数を使用する。 与える引数は、入力となる DEM だけである。 このアルゴリズムを使う際は、デフォルト値が研究の目的にあっているか確認する必要がある。58 "saga:sagawetnessindex" は、集水域、集水勾配、修正集水域、地形湿潤指数という4つのラスタを返す。

dem_wetness = qgis_run_algorithm("saga:sagawetnessindex", DEM = dem, 
                                 .quiet = TRUE)

その結果、dem_wetness、4つの出力へのファイルパスが記載されたリストが得られる。 qgis_as_terra() 関数で出力名を指定することで、選択した出力を読み出すことができる。

dem_wetness_twi = qgis_as_terra(dem_wetness$TWI)

Figure 10.3 の左パネルに出力された TWI マップを見ることができる。 地形湿潤指数は、単位なし。 数値が小さいほど水がたまらず、数値が大きいほど水がたまるエリアであることを示す。

また、デジタル標高モデルからの情報は、例えば、地形に分類することができる。地形は、斜面、尾根、谷などの地形を表す10のクラスからなる地形学的フォノタイプである (Jasiewicz and Stepinski 2013) 。 これらのクラスは、地滑りしやすさ、生態系サービス、人間の移動性、デジタル土壌マッピングなど、多くの研究で利用されている。

geomorphons のアルゴリズムのオリジナルの実装は GRASS GIS で作成され、qgisprocessのリストで "grass7:r.geomorphon" として見つけることができる。

grep("geomorphon", qgis_algo$algorithm, value = TRUE)
#> [1] "grass7:r.geomorphon"
qgis_show_help("grass7:r.geomorphon")
# 出力は非表示

geomorphons の計算には、入力 DEM (elevation) が必要で、オプションの引数でカスタマイズすることができる。 search – 視線距離を計算する長さ、および -m – 検索値を(セル数ではなく)メートル単位で提供することを指定するフラグが含まれている。 追加論点の詳細は、原著論文と GRASS GIS documentationに記載されている。

dem_geomorph = qgis_run_algorithm("grass7:r.geomorphon", elevation = dem, 
                                    `-m` = TRUE, search = 120, .quiet = TRUE)

出力される dem_geomorph$forms は、10個のカテゴリーからなるラスタファイルで、それぞれが地形形状を表している。 これを qgis_as_terra() でRに読み込んで可視化したり( Figure 10.3 の右パネル)、その後の計算で使ったりすることができる。

dem_geomorph_terra = qgis_as_terra(dem_geomorph$forms)

興味深いことに、Figure 10.3 に示すように、いくつかの地形と TWI 値の間にはつながりがある。 TWI値が最も大きいのは谷や窪地であり、最も小さいのは予想通り尾根であった。

Topographic wetness index (TWI, left panel) and geomorphons (right panel) derived for the Mongón study area.

FIGURE 10.3: Topographic wetness index (TWI, left panel) and geomorphons (right panel) derived for the Mongón study area.

10.3 SAGA GIS

System for Automated Geoscientific Analyses (SAGA ; Table 10.1 ) は、コマンドラインインタフェース ( Windows では saga_cmd.exe、Linux では単に saga_cmd ) を介して SAGA モジュールを実行する可能性を提供する(SAGA wiki on modules を参照)。 また、Pythonインターフェース(SAGA Python API )も用意されている。 Rsagacmd は、前者を使って R 内で SAGA を実行している。

この Section では、Rsagacmd を使用して、SAGA GIS の seeded region growing アルゴリズムを使って、2000年9月22日のペルーの Mongón 調査地域の正規化差分植生指数(normalized difference vegetation index, NDVI)の値が類似した地域を抽出する(Figure 10.4 左)。59

ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))

Rsagacmd を始めるには、saga_gis() 関数を実行する必要がある。 この関数は主に2つの目的がある。

  • 有効な SAGA-GIS ライブラリやツールへのリンクを含む新しいオブジェクトを動的に作成すること60
  • raster_backend (ラスタデータを扱う際に用いる R パッケージ)、vector_backend (ベクトルデータを扱う際に用いる R パッケージ)、cores (処理に用いる CPU コアの最大数、デフォルトは all) など、一般的なパッケージオプションを設定すること
library(Rsagacmd)
saga = saga_gis(raster_backend = "terra", vector_backend = "sf")

この saga オブジェクトは、利用可能なすべての SAGA ツールへの接続を含んでいる。 これはライブラリ(ツールのグループ)のリストとして構成されており、ライブラリの内部にはツールのリストがある。 どのツールにも $ 記号でアクセスできる(TAB キーで自動補完することが可能)。

シード領域拡大アルゴリズムは,主に2つのステップで動作する (adams_seed_1994?; Böhner, Selige, and Ringeler 2006)。 まず、指定されたサイズのローカルウィンドウにおいて、最も分散の小さいセルを見つけることで、初期セル(seed)が生成される. 次に、領域成長アルゴリズムを用いて、seed の近傍画素をマージし、均質な領域を作成する。

sg = saga$imagery_segmentation$seed_generation

上記の例では、まず imagery_segmentation ライブラリを開き、次にその seed_generation ツールを使用した。 また、次のステップでツールのコード全体を再入力しないように、sg オブジェクトに割り当てる。61 sg と入力することで、ツールの簡単な概要と、パラメータ、説明、およびデフォルトのデータフレームが表示される。 また、tidy(sg) を使用すると、パラメータのテーブルだけを取り出すことができる。 seed_generation ツールは、少なくとも1つの入力 (features) つまりラスタデータを必要とする。 また、初期ポリゴンのサイズを指定する band_width などの追加パラメータを提供することができる。

ndvi_seeds = sg(ndvi, band_width = 2)
#plot(ndvi_seeds$seed_grid)

この出力は、3つのオブジェクトからなるリストである。variance – 局所分散のラスタマップ、 seed_grid – 生成されたシードを含むラスタマップ、 seed_points – 生成されたシードを含む空間ベクトルオブジェクト。

2つ目のSAGA GISツールは seeded_region_growing である。62 seed_region_growing ツールは、前のステップで計算した seed_gridndvi ラスタオブジェクトの 2 つの入力を必要とする。 さらに、入力の特徴量を標準化するための normalizeneighbour (4または8-neighborhood)、 method などのパラメータを指定することができます。 最後のパラメータには、0 または 1 を指定することができる(ラスタセルの値とその位置に基づいて領域を成長させるか、値のみを成長させるか)。 このメソッドの詳細な説明は、 Böhner, Selige, and Ringeler (2006) を参照。

ここでは、method1 に変更するだけです。つまり、出力される地域は、NDVI 値の類似性に基づいてのみ作成されることを意味する。

srg = saga$imagery_segmentation$seeded_region_growing
ndvi_srg = srg(ndvi_seeds$seed_grid, ndvi, method = 1)
plot(ndvi_srg$segments)

このツールは、3つのオブジェクトのリストを返す。このツールは、 segments, similarity, table という 3 つのオブジェクトのリストを返す。 similarity オブジェクトは、シードと他のセルとの類似性を示すラスタであり、table は入力シードに関する情報を格納したデータフレームである。 最後に、ndvi_srg$segmentsは、結果として得られた領域(Figure 10.4 の右側のパネル)を表すラスタである。 これをポリゴンに変換するには、as.polygons()st_as_sf() を使用する(Section 6.5)。

ndvi_segments = as.polygons(ndvi_srg$segments) |> 
  st_as_sf()
Normalized difference vegetation index (NDVI, left panel) and NDVi-based segments derived using t he seeded region growing algorithm for the Mongón study area.

FIGURE 10.4: Normalized difference vegetation index (NDVI, left panel) and NDVi-based segments derived using t he seeded region growing algorithm for the Mongón study area.

結果として得られるポリゴン(セグメント)は、類似した値を持つ領域を表す。 また、クラスタリング(k-meansなど)、地域化(SKATER など)、教師あり分類法など、さまざまな手法でさらに大きなポリゴンに集約することができる。 演習で試すことができる。

R には、似たような値を持つポリゴン(いわゆるセグメント)を作成するという目的を達成するための他のツールもある。 いくつかの画像分割アルゴリズムを実行できる SegOptim パッケージ (Gonçalves et al. 2019) や、地理空間データを扱うためにスーパーピクセルアルゴリズム SLIC を実装した supercells (Nowosad and Stepinski 2022) などが含まれている。

10.4 GRASS GIS

米国陸軍建設工学研究所(USA-CERL)は、1982年から1995年にかけて、地理資源解析支援システム(GRASS) の中核となるシステムを作成した [Table 10.1 ; Neteler and Mitasova (2008)]。 アカデミアは1997年からこの作業を継続した。 SAGA と同様、グラスも当初はラスタ処理に注力し、その後、GRASS 6.0以降、高度なベクタ機能を追加している (Bivand, Pebesma, and Gómez-Rubio 2013)

GRASS は、入力データを GRASS GIS データベースに格納する。 ベクタデータに関して、GRASS はデフォルトでトポロジカル GIS、すなわち隣接するフィーチャのジオメトリを一度だけ保存する。 ベクトル属性の管理にはデフォルトで SQLite を用い、属性はキーによってジオメトリ、すなわち GRASS GIS データベースにリンクされる(GRASS GIS vector management)。

GRASS を使う前に、GRASS GIS データベース (同じく R から) をセットアップする必要があるが、このプロセスに少し戸惑うかもしれない。 まず、GRASS のデータベースは専用のディレクトリを必要とし、そのディレクトリには location を置く必要がある(詳しくは grass.osgeo.orgGRASS GIS Database ヘルプページを参照)。 location には、1 つのプロジェクトまたは 1 つの領域のジオデータが格納される。 通常、1つの場所の中に、異なるユーザや異なるタスクを参照するいくつかのマップセットを存在させることができる。 各 location には、PERMANENT マップセット(自動的に作成される必須のマップセット)もある。 プロジェクトのすべてのユーザと地理データを共有するために、データベース所有者は PERMANENT マップセットに空間データを追加することができる。 さらに、PERMANENT マップセットには、ラスタデータの投影法、空間範囲、およびデフォルトの解像度が格納される。 まとめると、GRASS GIS データベースは多くの location を含み(1つのロケーションのデータはすべて同じ CRS を持つ)、それぞれの location は多くのマップセット(データセットのグループ)を格納することができる。 GRASS 空間データベースシステムの詳細は、Neteler and Mitasova (2008)GRASS GIS quick start 。 こkでは R から手軽に GRASS を使うため link2GI パッケージを使う。しかし、GRASS GIS データベースを順番に作ることもできる。 作り方は GRASS within R を参照。

ここでは、GIScience における最も興味深い問題の一つである巡回セールスマン問題 を用いた rgrass を紹介する。 ある巡回セールスマンが24件の顧客を訪問したいとする。 さらに、自宅を起点とし、最短距離で25カ所を回りたいということであった。 この問題に対する最適解は一つであるが、考えられる解をすべてチェックすることは、現代のコンピュータでは(ほとんど)不可能である (Longley 2015)。 この場合、可能な解の数は (25 - 1)! / 2、すなわち24の階乗を2で割った数に相当する(2で割るのは、順方向と逆方向を区別しないため)。 1回の繰り返しがナノ秒でも、9837145 年間に相当する。 幸いなことに、この想像を絶する時間のごく一部で実行できる、巧妙でほぼ最適なソリューションがある。 GRASS GIS は、これらの解決策の一つを提供する(詳細は、 v.net.salesmanを参照)。 今回の使用例では、ロンドンの街角にある最初の25の自転車ステーション(顧客の代わり)間の最短経路を見つけたい(最初の自転車ステーションは、巡回セールスマンの自宅に相当すると単純に仮定する)。

data("cycle_hire", package = "spData")
points = cycle_hire[1:25, ]

自転車ハイヤーポイントのデータの他に、この地域の道路網が必要である。 OpenStreetMap から osmdata パッケージの助けを借りてダウンロードすることができる(Section 8.2 も参照)。 そのために、道路網のクエリ(OSM 言語では “highway” とラベル付けされている)を points の bounding box に制限し、対応するデータを sf オブジェクトとして読み込む。 osmdata_sf() は、複数の空間オブジェクト(点、線、ポリゴンなど)を含むリストを返すが、ここでは線オブジェクトとその関連IDのみを保持する。63

library(osmdata)
b_box = st_bbox(points)
london_streets = opq(b_box) |>
  add_osm_feature(key = "highway") |>
  osmdata_sf() 
london_streets = london_streets[["osm_lines"]]
london_streets = dplyr::select(london_streets, osm_id)

これでデータが揃ったので、次に GRASS のセッションを開始する。 幸い、link2GI パッケージの linkGRASS() を使えば、たった一行のコードで GRASS 環境をセットアップできる。 空間オブジェクトは、空間データベースの投影と範囲を決定するものである。 まず、linkGRASS() は、あなたのコンピュータにインストールされている全ての GRASS を検索する。 ここでは ver_selectTRUE に設定しているので、見つかった GRASS-installation の中から対話的に一つを選択することができる。 もし、インストールが一つしかない場合は、linkGRASS() が自動的にそれを選択する。 次に、linkGRASS() は、GRASS GIS への接続を確立する。

library(rgrass)
link2GI::linkGRASS(london_streets, ver_select = TRUE)

GRASS の geoalgorithm を使用する前に、GRASS の空間データベースにデータを追加する必要があります。 幸いなことに、便利な関数 write_VECT() がこれを代行してくれる。 (ラスタデータには write_RAST() を使用する。) この例では、最初の属性カラムのみを使用して、道路と自転車レンタルポイントデータを追加し、GRASS で london_streetspoints という名前を付けている。

write_VECT(terra::vect(london_streets), vname = "london_streets")
write_VECT(terra::vect(points[, 1]), vname = "points")

rgrass パッケージは、入力と出力が terra オブジェクトであることを想定している。 したがって、 write_VECT() を使用するためには、 vect() 関数を使用して sf 空間ベクトルを terraSpatVector に変換する必要がある。64

現在、両方のデータセットが GRASS GIS のデータベースに存在しています。 ネットワークの解析を行うには、トポロジカルクリーンな道路ネットワークが必要です。 GRASS の "v.clean" は、重複、小角、ダングルの除去などを行う。 ここでは、後続のルーティングアルゴリズムが実際に交差点で右折または左折できるように、各交差点で改行し、その出力を streets_clean という名前の GRASS オブジェクトに保存している。

execGRASS(cmd = "v.clean", input = "london_streets", output = "streets_clean",
          tool = "break", flags = "overwrite")
To learn about the possible arguments and flags of the GRASS GIS modules you can you the help flag. For example, try execGRASS("g.region", flags = "help").

サイクリング・ステーションのいくつかのポイントは、正確に街路セグメント上に位置しない可能性がある。 しかし、それらの間の最短ルートを見つけるために、それらを最も近い道路セグメントに接続する必要がある。 "v.net"のconnect-operatorはまさにこれを行う。 その出力を streets_points_con に保存する。

execGRASS(cmd = "v.net", input = "streets_clean", output = "streets_points_con",
          points = "points", operation = "connect", threshold = 0.001,
          flags = c("overwrite", "c"))

得られたクリーンなデータセットは "v.net.salesman" アルゴリズムの入力となり、最終的にすべての自転車レンタルステーション間の最短経路を見つけることができる。 その引数の一つが center_cats で、これは入力として数値の範囲を必要とする。 この範囲は、最短ルートを計算するためのポイントを表している。 ここでは、すべての自転車ステーション間の経路を計算したいので、1-25に設定しておく。 巡回セールスマンアルゴリズムの GRASS ヘルプページを参照するには、 execGRASS("g.manual", entry = "v.net.salesman") を実行する。

execGRASS(cmd = "v.net.salesman", input = "streets_points_con",
          output = "shortest_route", center_cats = paste0("1-", nrow(points)),
          flags = "overwrite")

この結果を見るために、結果をRに読み込み、形状だけを残してsf-objectに変換し、mapviewパッケージを使って可視化しました(Figure 10.5 と Section 9.4)。

route = read_VECT("shortest_route") |>
  st_as_sf() |>
  st_geometry()
mapview::mapview(route) + points
Shortest route (blue line) between 24 cycle hire stations (blue dots) on the OSM street network of London.

FIGURE 10.5: Shortest route (blue line) between 24 cycle hire stations (blue dots) on the OSM street network of London.

その際、いくつか注意すべき点がある。

  • GRASSの空間データベースを使えば、より高速に処理できる。 つまり、最初に地理データを書き出しただけです。 そして、新しいオブジェクトを作成し、最終結果だけを R にインポートした。 現在利用可能なデータセットを調べるには、 execGRASS("g.list", type = "vector,raster", flags = "p") を実行してみてください。
  • また、R から既にある GRASS 空間データベースにアクセスすることも可能である。 R にデータをインポートする前に、いくつかの(空間)サブセット化を行いたい場合がある。 ベクトルデータには "v.select""v.extract" を使用する。 "db.select" を使用すると、対応するジオメトリを返さずにベクトルレイヤーの属性テーブルのサブセットを選択することができる。
  • また、実行中の GRASS のセッションから R を起動することもできます (詳細は Bivand, Pebesma, and Gómez-Rubio 2013 を参照)
  • GRASS で提供されて入るジオアルゴリズムの素晴らしいドキュメントは、 GRASS online help または execGRASS("g.manual", flags = "i") を参照。

10.5 いつ、何を使うべきか?

R-GIS のインターフェースは、個人の好みや作業内容、GIS の使い方に依存するため、一概にお勧めすることはできないし、研究分野にもよるだろう。 前述の通り、SAGA は大規模(高解像度)ラスタデータセットの高速処理に特に優れており、水文学者、気候学者、土壌学者に頻繁に利用されている (Conrad et al. 2015)。 一方、GRASS GIS は、トポロジーに基づく空間データベースをサポートする唯一の GIS であり、ネットワーク分析だけでなくシミュレーション研究にも特に有用である。 QGIS は、GRASS-GIS や SAGA-GIS と比較して、特に初めて GIS を使う方にとって使いやすく、おそらく最も人気のあるオープンソースのGISだと思われる。 したがって、qgisprocess は、ほとんどのユースケースに適切な選択である。 その主なメリットは

  • 複数の GIS に統一的にアクセスできるため、重複した機能を含む 1000 以上のジオアルゴリズム ( Table 10.1 ) を提供。例えば、QGIS、SAGA、GRASS などのジオアルゴリズムを使ってオーバーレイ操作を実行することが可能である。
  • データフォーマットの自動変換(SAGAは .sdat グリッドファイル、GRASS は独自のデータベースフォーマットを使用するが、対応する変換は QGIS が行う。)
  • 地理的な R オブジェクトを QGIS ジオアルゴリズムに自動的に渡し、R に戻すことができる。
  • 名前付き引数、デフォルト値の自動取得をサポートする便利な機能(rgrass からインスパイアされた)

もちろん、他の R-GIS ブリッジを使用した方が良いケースもある。 QGIS は、複数の GIS ソフトウェアパッケージへの統一インターフェースを提供する唯一の GIS であるが、対応するサードパーティのジオアルゴリズムのサブセットへのアクセスしか提供しない(詳細については、Muenchow, Schratz, and Brenning (2017) を参照)。 したがって、SAGA と GRASS の関数一式を使用するには、RSAGArgrass 以外は使わない方が良い。 最後に、地形データ、空間データベース管理機能(マルチユーザーアクセスなど)が必要な場合は、GRASSの利用を勧める。 また、ジオデータベースを用いてシミュレーションを行いたい場合 (Krug, Roura-Pascual, and Richardson 2010)qgisprocess が、呼び出しごとに常に新しい GRASS セッションを開始するので、rgrass を直接使用してみよう。

なお、スクリプティング・インターフェースを持つ GIS ソフトウェアパッケージは以下のように数多くあるが、これらを利用できる専用の R パッケージはない:gvSig、OpenJump、Orfeo Toolbox。65

10.6 その他のブリッジ

この章では、当面、デスクトップ GIS ソフトウェアへの R インタフェースに焦点を当てる。 これらの橋を強調したのは、専用の GIS ソフトウェアがよく知られており、地理データを理解するための一般的な「入り口」であることがある。 また、多くのジオアルゴリズムにアクセスすることができる。

その他の「ブリッジ」には、空間ライブラリへのインタフェース( Section 10.6.1 は R から GDAL CLI にアクセスする方法を示している)、空間データベース ( Section 10.6.2 参照)、ウェブマッピングサービス( Chapter 9 参照)などがある。 このセクションでは、可能なことのほんの一端を紹介する。 システムから他のプログラムを呼び出したり、他の言語と統合したり(特に Rcppreticulate を介して)、R の の柔軟性のおかげで、他の多くのブリッジが可能になっている。 目的は包括的なものではなく、章の冒頭で Sherman (2008) が引用した「柔軟性とパワー」にアクセスする他の方法を示すことである。

10.6.1 GDALへのブリッジ

Chapter 8 で述べたように、GDAL は多くの地理データ形式をサポートする低レベルのライブラリである。 GDAL は非常に効果的なので、ほとんどのGISプログラムは、車輪の再発明や特注の読み書きコードを使用するのではなく、地理データのインポートとエクスポートのためにバックグラウンドで GDAL を使用している。 しかし、GDAL が提供するのは、データ入出力だけではない。 ベクタデータとラスタデータの geoprocessing tools 、ラスタデータをオンラインで提供するためのタイルを作成する機能、ベクタデータの高速ラスタ化がある。 GDAL はコマンドラインツールであるため、R からは system() コマンドからアクセスすることができる。

以下のコードは、この機能を実現するものである。 linkGDAL() は、GDAL が動作しているコンピュータを検索し、実行ファイルの場所を PATH 変数に追加して、GDAL を呼び出せるようにする(Windows で通常必要になる)。

link2GI::linkGDAL()

これで、system() 関数を使用して、任意の GDAL ツールを呼び出すことができる。 例えば、ogrinfo は、ベクタデータセットのメタデータを提供する。 ここでは、このツールに2つのフラグを追加して呼び出する。 -al は全レイヤーの全フィーチャをリストアップし、-so は要約のみを取得する(完全なジオメトリのリストではない)。

our_filepath = system.file("shapes/world.gpkg", package = "spData")
cmd = paste("ogrinfo -al -so", our_filepath)
system(cmd)
#> INFO: Open of `.../spData/shapes/world.gpkg'
#>       using driver `GPKG' successful.
#> 
#> Layer name: world
#> Geometry: Multi Polygon
#> Feature Count: 177
#> Extent: (-180.000000, -89.900000) - (179.999990, 83.645130)
#> Layer SRS WKT:
#> ...

その他、よく使われるGDALのツールは以下の通り

  • gdalinfo: ラスタデータセットのメタデータを提供
  • gdal_translate: 異なるラスタファイルフォーマット間の変換
  • ogr2ogr: 異なるベクターファイルフォーマット間で変換
  • gdalwarp: ラスタデータセットの再投影、変換、切り抜き (clip)
  • gdaltransform: 座標変換

GDALツールの全リストとそのヘルプファイルは https://gdal.org/programs/

link2GI が提供する GDAL への「リンク」は、R やシステムの CLI からより高度な GDAL の作業を行うための基盤として利用することができるだろう。 TauDEM (http://hydrology.usu.edu/taudem) や Orfeo Toolbox (https://www.orfeo-toolbox.org/) は、コマンドラインインタフェースを提供する空間データ処理ライブラリ/プログラムである – 上記の例は、Rを介してシステムのコマンドラインからこれらのライブラリにアクセスする方法である。 これは、新しいRパッケージという形で、これらのライブラリへの適切なインタフェースを作成するための出発点となる可能性がある。

しかし、新しいブリッジを作成するプロジェクトに飛び込む前に、既存のRパッケージのパワーと、system() の呼び出しがプラットフォームに依存しない(一部のコンピュータで失敗する)可能性があることを認識しておくことが重要である。 さらに、sf は GDAL , GEOS , PROJ が提供するパワーのほとんどを Rcpp が提供する R/C++ インターフェイスを介して R にもたらし、system() の呼び出しを回避しているのである。

10.6.2 空間データベースへのブリッジ

空間データベース管理システム(空間DBMS)は、空間および非空間データを構造化して保存する。 大規模なデータの集合を、一意の識別子(主キーと外部キー)および暗黙のうちに空間を介して関連するテーブル(エンティティ)に整理することができる(たとえば、空間結合を考えてみてみよう)。 地理的なデータセットはすぐに大きくなったり、乱雑になったりする傾向があるため、この機能は便利である。 データベースは、空間および非空間フィールドに基づく大規模なデータセットの保存とクエリを効率的に行うことができ、マルチユーザーアクセスとトポロジーのサポートを提供する。

最も重要なオープンソースの空間データベースは PostGIS である (Obe and Hsu 2015)66 PostGIS のような空間 DBMS への R ブリッジは重要で、数ギガバイトの地理データを RAM にロードすることなく、、R セッションをクラッシュさせる可能性があるような巨大なデータストアにアクセスできる。 このセクションの残りの部分では、PostGIS in Action, Second Edition の “Hello real world” に基づいて、R から PostGIS を呼び出す方法を紹介する (Obe and Hsu 2015)67

QGIS Cloud (https://qgiscloud.com/) にある PostgreSQL/PostGIS データベースにアクセスしているため、この後のコードはインターネット接続している必要がある。68 最初のステップは、データベース名、ホスト名、およびユーザー情報を指定して、データベースへの接続を作成することである。

library(RPostgreSQL)
conn = dbConnect(drv = PostgreSQL(), 
                 dbname = "rtafdf_zljbqm", host = "db.qgiscloud.com",
                 port = "5432", user = "rtafdf_zljbqm", password = "d3290ead")

新しいオブジェクト conn は、Rセッションとデータベースの間のリンクを確立したに過ぎない。 データを保存することはない。

多くの場合、最初の質問は「データベースからどのテーブルが見つかるか」である。 これには、dbListTables()、次のように答えることができる。

dbListTables(conn)
#> [1] "spatial_ref_sys" "topology"        "layer"           "restaurants"    
#> [5] "highways" 

答えは、5つのテーブルである。 ここでは、restaurantshighways のテーブルのみを対象としている。 前者は米国内のファストフード店の位置を、後者は米国の主要な高速道路を表している。 テーブルで利用可能な属性について調べるには、dbListFields を実行する。

dbListFields(conn, "highways")
#> [1] "qc_id"        "wkb_geometry" "gid"          "feature"     
#> [5] "name"         "state"   

さて、利用可能なデータセットがわかったところで、いくつかのクエリを実行し、データベースに質問することができる。 クエリーは、データベースが理解できる言語(通常はSQL)で提供される必要がある。 最初のクエリは、highways テーブルから Maryland 州の US Route 1 ( MD ) を選択する。 なお、read_sf() は、データベースへのオープンな接続とクエリが提供されれば、データベースから地理データを読み込むことができる。 さらに、read_sf() は、どの列がジオメトリを表すかを知る必要がある(ここでは、wkb_geometry )。

query = paste(
  "SELECT *",
  "FROM highways",
  "WHERE name = 'US Route 1' AND state = 'MD';")
us_route = read_sf(conn, query = query, geom = "wkb_geometry")

この結果、MULTILINESTRING 型の us_route という名前の sf オブジェクトが生成される。

また、前述したように、非空間的な質問だけでなく、空間的な性質をもとにデータセットを問い合わせることも可能である。 これを示すために、次の例では選択した高速道路(Figure 10.6)の周囲に 35 km (35,000 m) のバッファを追加している。

query = paste(
  "SELECT ST_Union(ST_Buffer(wkb_geometry, 35000))::geometry",
  "FROM highways",
  "WHERE name = 'US Route 1' AND state = 'MD';")
buf = read_sf(conn, query = query)

なお、これはすでにお馴染みの関数(ST_Union()ST_Buffer() )を使った空間クエリであった。 また、sf パッケージにもあるが、こちらは小文字で書かれている( st_union() , st_buffer() )。 実際、sf パッケージの関数名は、PostGIS の命名規則にほぼ従っている。69

最後のクエリは、35 km のバッファーゾーン(Figure 10.6)内にあるすべてのハーディーズレストラン(HDE)を検索する。

query = paste(
  "SELECT *",
  "FROM restaurants r",
  "WHERE EXISTS (",
  "SELECT gid",
  "FROM highways",
  "WHERE",
  "ST_DWithin(r.wkb_geometry, wkb_geometry, 35000) AND",
  "name = 'US Route 1' AND",
  "state = 'MD' AND",
  "r.franchise = 'HDE');"
)
hardees = read_sf(conn, query = query)

空間 SQL クエリの詳細な説明は Obe and Hsu (2015) を参照。 最後に、次のようにデータベース接続を閉じるのがよい方法である。70

RPostgreSQL::postgresqlCloseConnection(conn)
Visualization of the output of previous PostGIS commands showing the highway (black line), a buffer (light yellow) and four restaurants (red points) within the buffer.

FIGURE 10.6: Visualization of the output of previous PostGIS commands showing the highway (black line), a buffer (light yellow) and four restaurants (red points) within the buffer.

PostGIS とは異なり、sfは空間ベクタデータのみをサポートしている。 PostGIS データベースに格納されたラスタデータを照会・操作するには、rpostgisパッケージ(Bucklin and Basille 2018)、または PostGIS インストールの一部に含まれる rastertopgsql などのコマンドラインツールを使用する必要がある。

このサブセクションでは、PostgreSQL/PostGIS の簡単な紹介にとどめる。 それでも、地理的および非地理的データを空間 DBMS で保存しながら、さらなる(地理)統計解析に必要なそれらのサブセットだけを R のグローバル環境にアタッチするという実践を奨励したい。 提示された SQL クエリのより詳細な説明と PostgreSQL/PostGIS 一般のより包括的な紹介は Obe and Hsu (2015) を参照。 PostgreSQL/PostGIS は、非常に難解なオープンソースの空間データベースである。 しかし、軽量なデータベースエンジンである SQLite/SpatiaLite や、バックグラウンドで SQLite を使用する GRASS も同様と言える(section (grass) 参照)。

If your datasets are too big for PostgreSQL/PostGIS and you require massive spatial data management and query performance, it may be worth exploring large-scale geographic querying on distributed computing systems. Such systems are outside the scope of this book but it worth mentioning that open source software providing this functionality exists. Prominent projects in this space include GeoMesa and Apache Sedona, formerly known as GeoSpark (Huang et al. 2017), which has and R interface provided by the apache.sedona package.

10.7 クラウドへのブリッジ

近年、インターネット上では、クラウド技術の利用が目立ってきている。 この中には、空間データの保存や処理に利用されることも含まれている。 Amazon Web Services, Microsoft Azure / Planetary Computer, Google Cloud Platform などの主要なクラウドコンピューティングプロバイダは、Sentinel-2 アーカイブのようなオープンな地球観測データの巨大なカタログを彼らのプラットフォーム上で提供している。 私たちは R を使い、これらのアーカイブから直接データに接続し、処理することができる。理想的には、同じクラウドや地域のマシンから接続することができる。

このような画像アーカイブをクラウド上でより簡単に、より効率的に利用するために、SpatioTemporal Asset Catalog (STAC), cloud-optimized GeoTIFF (COG) 画像形式、データキューブが有望視されていまる。 Section 10.7.1 では、これらの個々の開発について紹介し、Rからどのように利用できるかを簡単に説明する。

大規模なデータアーカイブをホストするだけでなく、地球観測データを処理するための多くのクラウドベースのサービスもここ数年で開始された。 その中には、R を含むプログラミング言語と様々なクラウドサービスとの間の統一的なインタフェースである OpenEO イニシアチブも含まれている。 OpenEO の詳細については、“openEO”の項を参照。

10.7.1 クラウドの STAC、COGs 、その他のデータキューブ

STAC(SpatioTemporal Asset Catalog)は、時空間データの汎用記述フォーマットで、画像、合成開口レーダー(SAR)データ、点群など、クラウド上の様々なデータセットの記述に使用されている。 STAC-API は、単純な静的カタログ記述の他に、カタログのアイテム(画像など)を空間、時間、その他のプロパティで照会するウェブサービスを提供している。 R では、rstacパッケージ(Simoes, Souza, et al. 2021)が STAC-API エンドポイントに接続し、アイテムを検索することができる。 以下の例では、Sentinel-2 Cloud-Optimized GeoTIFF (COG) dataset on Amazon Web Services から、事前に定義した関心領域と時間に交差するすべての画像を要求している。 結果は、見つかったすべての画像とそのメタデータ(雲量など)、および AWS 上の実際のファイルを指す URL を含んでいます。

library(rstac)
# Connect to the STAC-API endpoint for Sentinel-2 data
# and search for images intersecting our AOI
s = stac("https://earth-search.aws.element84.com/v0")
items = s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(7.1, 51.8, 7.2, 52.8), 
              datetime = "2020-01-01/2020-12-31") |>
  post_request() |> items_fetch()

クラウドストレージはローカルのハードディスクとは異なり、従来の画像ファイル形式はクラウドベースのジオプロセシングではうまく機能しない。 クラウドに最適化された GeoTIFF は、GeoTIFF の一種であり、クラウドストレージから画像の一部分のみを効率的に読み込むことができる。 そのため、画像の矩形部分や低解像度の画像の読み込みが非常に効率的になる。 GDAL(およびそれを使ったパッケージ)はすでに COG を扱うことができるので、R ユーザーであれば COG を扱うために何かをインストールする必要はない。ただし、データ提供者のカタログを閲覧する際には、COG が利用可能であることが大きなプラスになることを覚えておこう。

領域が大きい時、要求された画像を扱うのはまだ比較的困難である。それらは異なる地図投影を使用することがあり、空間的に重なることがあり、空間解像度はしばしばスペクトルバンドに依存する。 gdalcubes パッケージ (Appel and Pebesma 2019) は、個々の画像から抽象化し、画像コレクションを4次元データキューブとして作成し処理するために使用することができる。

以下のコードは、前回の STAC-API 検索で返された Sentinel-2 画像から、低解像度(250m)の最大 NDVI コンポジットを作成する最小限の例を示している。

library(gdalcubes)
# Filter images from STAC response by cloud cover 
# and create an image collection object
collection = stac_image_collection(items$features, 
                  property_filter = function(x) {x[["eo:cloud_cover"]] < 10})
# Define extent, resolution (250m, daily) and CRS of the target data cube
v = cube_view(srs = "EPSG:3857", extent = collection, dx = 250, dy = 250,
              dt = "P1D") # "P1D" is an ISO 8601 duration string
# Create and process the data cube
cube = raster_cube(collection, v) |>
  select_bands(c("B04", "B08")) |>
  apply_pixel("(B08-B04)/(B08+B04)", "NDVI") |>
  reduce_time("max(NDVI)")
# gdalcubes_options(parallel = 8)
# plot(cube, zlim = c(0, 1))

クラウドカバーによる画像のフィルタリングを行うために、画像コレクションを作成する際に各 STAC 結果アイテムに適用されるプロパティフィルタ関数を提供している。 この関数は、画像の利用可能なメタデータを入力リストとして受け取り、関数がTRUEを返す画像のみを考慮するような単一の論理値を返す。 この場合、10% 以上のクラウドカバーがある画像は無視する。 詳しくは、こちらのOpenGeoHubサマースクール2021で発表したチュートリアルを参照。

STAC、COGs、data cube を組み合わせて、衛星画像の(大規模)コレクションをクラウド上で解析するクラウドネイティブワークフローを形成する. これらのツールは、例えば、大規模な地球観測データの土地利用や土地被覆の分類を可能にする sits R パッケージのバックボーンを既に形成している。 このパッケージは、クラウドサービスで利用可能な画像コレクションからEOデータキューブを構築し、様々な機械学習アルゴリズムを用いてデータキューブの土地分類を実行するものである。 sits の詳細については、https://e-sensing.github.io/sitsbook/ をご覧いただくか、関連記事 (Simoes, Camara, et al. 2021) を参照。

10.7.2 openEO

地球観測データを処理するためのクラウドサービスは、大規模なデータアーカイブの他に、ここ数年で数多く提供されるようになった。 OpenEO (Schramm et al. 2021) は、データ処理のための共通言語を定義することによって、クラウドサービス間の相互運用性を支援するイニシアチブである。 最初のアイデアはr-spatial.org blog postで説明されており、ユーザーができるだけ少ないコード変更で簡単にクラウドサービス間を変更できるようにすることを目的としている。 標準化プロセスでは、データへのインタフェースとして多次元データキューブモデルを使用している。 8種類のバックエンドの実装が用意されており(https://hub.openeo.org)、ユーザーは R、Python、JavaScript、QGIS、Web エディタで接続し、コレクションに対してプロセスを定義(およびチェーン)することができる。 バックエンドによって機能や利用できるデータが異なるため、openeo R パッケージ (Lahn 2021) は接続されたバックエンドから利用できるプロセスとコレクションを動的にロードする。

その後、ユーザーは画像コレクションのロード、プロセスの適用と連鎖、ジョブの送信、結果の探索とプロットを行うことができる。

以下のコードは、openEO platform backend に接続し、利用可能なデータセット、プロセス、出力フォーマットを要求し、Sentinel-2データから最大 NDVI 画像を計算するプロセスグラフを定義し、最後にバックエンドにログインした後にグラフを実行する。 openEO プラットフォームのバックエンドには無料版があり、既存の機関やソーシャルプラットフォームのアカウントから登録することが可能である。

library(openeo)
con = connect(host = "https://openeo.cloud")
p = processes() # load available processes
collections = list_collections() # load available collections
formats = list_file_formats() # load available output formats
# Load Sentinel-2 collection
s2 = p$load_collection(id = "SENTINEL2_L2A",
                       spatial_extent = list(west = 7.5, east = 8.5,
                                             north = 51.1, south = 50.1),
                       temporal_extent = list("2021-01-01", "2021-01-31"),
                       bands = list("B04","B08")) 
# Compute NDVI vegetation index
compute_ndvi = p$reduce_dimension(data = s2, dimension = "bands",
                                  reducer = function(data, context) {
                                    (data[2] - data[1]) / (data[2] + data[1])
                                  })
# Compute maximum over time
reduce_max = p$reduce_dimension(data = compute_ndvi, dimension = "t",
                                reducer = function(x, y) {max(x)})
# Export as GeoTIFF
result = p$save_result(reduce_max, formats$output$GTiff)
# Login, see https://docs.openeo.cloud/getting-started/r/#authentication
login(login_type = "oidc", provider = "egi",
      config = list(client_id= "...", secret = "..."))
# Execute processes
compute_result(graph = result, output_file = tempfile(fileext = ".tif"))

10.8 演習

E1. Compute global solar irradiation for an area of system.file("raster/dem.tif", package = "spDataLarge") for March 21 at 11:00 AM using the r.sun GRASS GIS through qgisprocess.

E2. Compute catchment area and catchment slope of system.file("raster/dem.tif", package = "spDataLarge") using Rsagacmd.

E3. Continue working on the ndvi_segments object created in the SAGA GIS section. Extract average NDVI values from the ndvi raster and group them into six clusters using kmeans(). Visualize the results.

E4. Attach data(random_points, package = "spDataLarge") and read system.file("raster/dem.tif", package = "spDataLarge") into R. Select a point randomly from random_points and find all dem pixels that can be seen from this point (hint: viewshed can be calculated using GRASS GIS). Visualize your result. For example, plot a hillshade, the digital elevation model, your viewshed output, and the point. Additionally, give mapview a try.

E5. Use gdalinfo via a system call for a raster file stored on disk of your choice. What kind of information you can find there?

E6. Use gdalwarp to decrease the resolution of your raster file (for example, if the resolution is 0.5, change it into 1). Note: -tr and -r flags will be used in this exercise.

E7. Query all Californian highways from the PostgreSQL/PostGIS database living in the QGIS Cloud introduced in this chapter.

E8. The ndvi.tif raster (system.file("raster/ndvi.tif", package = "spDataLarge")) contains NDVI calculated for the Mongón study area based on Landsat data from September 22nd, 2000. Use rstac, gdalcubes, and terra to download Sentinel-2 images for the same area from 2020-08-01 to 2020-10-31, calculate its NDVI, and then compare it with the results from ndvi.tif.