2 Rで地理データ

必須パッケージ

本書の最初の実用的な章であるため、ソフトウェア要件がある。 R の最新版(R 4.2.0 またはそれ以降のバージョン)がインストールされているコンピュータにアクセスする必要がある。 文章を読むだけでなく、各章のコードを実行して、ジオコンピューティングのスキルを身につけることを勧める。

R のスクリプト、出力、その他ジオコンピュテーションに関連するものを保存するために、コンピュータに新しいフォルダを作成することから始めるとよいだろう。 また、学習をサポートするため、ソースコードダウンロード またはクローンしておくと良い。 R コードを書く/実行する/テストする際には、 RStudio(ほとんどの人に推奨)または VS Code などの統合開発環境 (integrated development environment, IDE) をインストールすることを強く勧める6

R を初めて使う方は、R のコードを使ったジオコンピュテーションに入る前に、Garrett Grolemund の Hands on Programming with R や Claudia Engel の Introduction to R などの R 入門リソースに従うことを勧める。 作業内容を整理し(例:RStudio プロジェクト)、スクリプトに chapter-02-notes.R などのわかりやすい名前を付けて、学習しながらコードを記録するとよい。

セットアップが完了したら、いよいよコードを実行する。 以下のパッケージがまだインストールされていなければ、まず最初にこの章で使用する基礎的な R パッケージを、以下のコマンドでインストールする:7

install.packages("sf")
install.packages("terra")
install.packages("spData")
install.packages("spDataLarge", repos = "https://nowosad.r-universe.dev")
CRAN に書かれている R の解説を推奨する。 Mac や Linux を使っている場合、上の sf をインストールするコマンドは初回ではうまくいかないことがある。 Mac/Linux には、パッケージの README に記載されている「システム要件」がある。 また、ブログ rtask.thinkr.fr の記事 Installation of R 4.2 on Ubuntu 22.04.1 LTS など、他の OS 向けの解説もオンラインで見つけることができる。

本書の第1部を再現するために必要なパッケージは、以下のコマンドでインストールできる。 remotes::install_github("geocompr/geocompkg") 。 このコマンドは、remotes パッケージの関数 install_packages() を使用して、GitHub コードホスティング、バージョン、およびコラボレーション プラットフォームでホストされているソース コードをインストールする。 以下のコマンドを実行すると、本書全体を再現するために必要なすべての依存関係がインストールされる(警告:数分かかる場合がある)。 remotes::install_github("geocompr/geocompkg", dependencies = TRUE)

この章で紹介するコードを実行するために必要なパッケージは、以下のように library() 関数で「読み込み」(厳密には attach)をすることができる。

library(sf)         # ベクタデータ用クラスと関数
#> Linking to GEOS 3.10.3, GDAL 3.5.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(sf) の出力では、GEOS などの主要な地理ライブラリがどのバージョンで使用されているか表示される。これらのライブラリについては、Section 2.2.1 で説明する。

library(terra)     # ラスタデータ用クラスと関数

インストールされた他のパッケージには、この本で使用されるデータが含まれている。

library(spData)       # 地理データをロード
library(spDataLarge)  # 大きい地理データをロード

2.1 イントロダクション

この章では、基本的な地理データモデルであるベクタとラスタについて簡単に説明する。(訳注:本書では、GIS で用いられる vector は「ベクタ」と訳し、R のデータ形式の vector は「ベクトル」と訳す。) それぞれのデータモデルの背景にある理論や、データモデルが得意とする分野を紹介し、R での実際に演習する。

ベクタデータモデルは、点、線、ポリゴンを使って世界を表現する。 これらのデータセットには、境界が明確に定義されているため、ベクタデータセットは通常、高い精度を持つ(ただし、Section 2.5 で見るように、必ずしも正確ではない)。 ラスタデータモデルは、表面を一定の大きさのセルに分割している。 ラスタデータセットは、ウェブマッピングで使用される背景画像の基本であり、航空写真や衛星リモートセンシング装置の発明以来、地理データの重要なソースである。 ラスタは、空間的に特定のフィーチャ(feature、地物とも訳される)を特定の解像度で集約したもので、空間的に一貫性があり、拡張性がある(世界中の多くのラスタデータセットが利用可能)。

どちらを使うべきか? その答えは、アプリケーションの領域によって変わってくる。

  • 社会科学では、人間の居住地が不連続な境界線を持つ傾向があるため、ベクタデータが主流となる傾向がある
  • 環境科学において、多くの場合リモートセンシングデータへの依存から、ラスタが主流となっている

分野によっては重複する部分も多く、ラスタデータセットとベクタデータセットを併用することも可能である。 例えば、生態学者や人口統計学者などは、ベクタデータとラスタデータの両方を使うのが一般的である。 さらに、この2つの形式を相互に変換することも可能である(Chapter 6 参照)。 ベクタデータとラスタデータのどちらを使用するかは、次の章で説明するように、使用する前に基本的なデータモデルを理解することが重要である。 本書では、ベクタデータとラスタデータセットを扱うために、それぞれ sfterra パッケージを使用している。

2.2 ベクタデータ

という単語は、次の二つの意味があるため注意が必要。 本書では便宜上、地理学データ形式を「ベクタ」、R のクラスを「ベクトル」と表記を分ける。 前者はデータモデルを指し、後者は data.framematrix などと同じく R のクラスの一つである。 両者には、もちろん関連性もある。ベクタで重要な空間座標は、R ではベクトルオブジェクトで表される。

地理ベクタデータモデルは、座標参照系 (coordinate reference system, CRS) 内に位置する点に基づいている。 点は、バス停の位置のような独立したフィーチャ(地物)を表すことも、線やポリゴンのような複雑な幾何学的形状を形成するために連結されることもある。 点ジオメトリ(geometry、ISO 19125 では「形状」や「幾何形状」と訳されているが、本書ではジオメトリとする。)はほとんどの場合2次元のみで構成される(3次元の CRS には、海抜高度を表す \(z\) の値が追加される)。

このシステムにおいて、例えばロンドンは、座標 c(-0.1, 51.5) で表すことができる。 つまり、その位置は原点から東に -0.1 度、北に 51.5 度である。 この場合の原点は、地理的(‘lon/lat’)CRS の経度 (longitude) 0度(本初子午線)と緯度 (latitude) 0度(赤道)にある(Figure 2.1、左図)。 また、同じ点を投影した CRS では、 British National Grid の「東経/北緯」の値はおよそc(530000, 180000) で、ロンドンが CRS の原点から 530 km 、180 km に位置することを意味している。 これは視覚的にも確認することができ、幅100kmの灰色の格子線に囲まれた正方形の領域である「ボックス」が5個強、ロンドンを表す点と原点を分けている(Figure 2.1、右のパネル)。

National Grid の の起点が South West Peninsular の先の海上にあるため、英国内のほとんどの場所で正の東経と北緯の値を持つことになる。8 CRS については、Section 2.4 と Sections 7 で説明するが、この章では、座標が原点からの距離を表す2つの数値からなり、通常は \(x\)\(y\) の次元であることを知っていれば十分である。

原点(青丸)を基準にロンドン(赤X)の位置を表したベクトル(点)データの図解。左図は、緯度経度0°を原点とする地理的なCRS。右図は、South West Peninsula の西側の海を原点とする投影CRSを表している。原点(青丸)を基準にロンドン(赤X)の位置を表したベクトル(点)データの図解。左図は、緯度経度0°を原点とする地理的なCRS。右図は、South West Peninsula の西側の海を原点とする投影CRSを表している。

FIGURE 2.1: 原点(青丸)を基準にロンドン(赤X)の位置を表したベクトル(点)データの図解。左図は、緯度経度0°を原点とする地理的なCRS。右図は、South West Peninsula の西側の海を原点とする投影CRSを表している。

sf は、地理ベクタデータ用のクラスと、以下に説明する地理計算のための重要な低レベルライブラリへのコマンドラインインターフェイスを一貫した方法で提供する。

  • GDAL は、Chapter 8 でカバーされる広範な地理データ形式の読み取り、書き込み、操作のためのものである。
  • PROJ は、 Chapter 7 で扱う内容の根幹をなす座標系変換のための強力なライブラリ。
  • GEOS は、投影型 CRS を持つデータに対してバッファや重心の計算などの操作を行う平面ジオメトリエンジン、Chapter 5 でカバーする。
  • S2 は、Google が開発した C++ で書かれた球面幾何学エンジンである。 s2 パッケージ経由で、この章の Section 2.2.9 と Chapter 7 でカバーする。

これらのインターフェースに関する情報は、パッケージが最初にロードされたときに sf によって表示される。この章の最初にある library(sf) コマンドの下に表示されるメッセージ Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE は、リンクされている GEOS、GDAL、PROJ ライブラリのバージョン(コンピュータや時期によって変わる)と S2インターフェースがオンになっているかどうかという情報を教えてくれる。 今では当たり前のことであるが、そもそも再現性のあるジオコンピューティングを可能にするのは、様々な地理ライブラリとの緊密な連携に他ならない。

sf の特徴として、投影法未指定のデータで使用するデフォルトのジオメトリエンジンを変更することができる S2 のスイッチを切るには、sf::sf_use_s2(FALSE) というコマンドを使う。つまり、平面ジオメトリエンジン GEOS が、投影されていないデータに対するジオメトリ操作を含む、すべてのジオメトリ操作にデフォルトで使用されることになる。 Section 2.2.9 で見るように、平面幾何学は2次元の空間に基づいている。 GEOS などの平面ジオメトリエンジンは「平面」(投影)座標を、S2 などの球面ジオメトリエンジンは非投影(緯度経度)座標を想定している。

この章では、後続の章に備え、sf クラスを紹介する(GEOS インターフェースは Chapter 5 章で、GDAL インターフェースは Chapter 8 で説明する)。

2.2.1 シンプルフィーチャの紹介

シンプルフィーチャ(simple feature、ISO 19125 では、feature は「地物」、simple feature は「単純地物」と訳されているるが、本書ではシンプルフィーチャで統一する。)は、Open Geospatial Consortium (OGC) が開発・推奨するオープンスタンダードであり、その活動は後の章で再確認することになる(Section 8.5)。 シンプルフィーチャは、様々なジオメトリタイプを表現する階層的なデータモデルである。 仕様でサポートされている18種類のジオメトリのうち、地理研究の大部分で使用されているのは7種類のみである(Figure 2.2)を参照)。 これらのコアなジオメトリタイプは、R パッケージの sf で完全にサポートされている (Pebesma 2018)9

sf が完全にサポートするシンプルフィーチャ型

FIGURE 2.2: sf が完全にサポートするシンプルフィーチャ型

sf は、点、線、ポリゴン、およびそれらの「複合」バージョン(同じ種類のフィーチャをまとめて 1 つのフィーチャとするもの) という一般的なベクタフィーチャのタイプをすべて表現できる(ラスタデータクラスは sf ではサポートされていない)。 また、sf はジオメトリコレクションをサポートしており、1つのオブジェクトに複数のジオメトリタイプを格納することができる。 sf は、データクラスの sp (Pebesma and Bivand 2018)、GDAL と PROJ を通したデータ読み書きの rgdal (Bivand, Keitt, and Rowlingson 2018)、GEOS を通した空間演算の rgeos (Bivand and Rundel 2018) という 3 パッケージで提供されていたものと同じ機能(およびそれ以上)を提供する。

Chapter 1 の繰り返しになるが、R の地理パッケージは低レベルのライブラリとのインターフェイスに長い歴史があり、sf はこの伝統を受け継ぎ、最近のバージョンの幾何演算のための GEOS、地理データファイルの読み書きのための GDAL、投影座標参照系の表現と変換のための PROJ ライブラリと統一されたインターフェイスを提供する。 s2 すなわち「Google の球面幾何学ライブラリ s2 への R インタフェース」を通して、sf は高速で正確な「非平面形状の測定と操作」 (Bivand 2021) にアクセスすることができる。 2021年6月に公開された sf バージョン 1.0.0 からは、地理(経度・緯度)座標系を持つジオメトリに対してデフォルトs2 機能が使われるようになった。これは、sf 独自の特徴であり、Python パッケージ GeoPandas など、ジオメトリ演算に GEOS にしか対応していない空間ライブラリとは異なる。 s2 については、以降の章で説明する。

地理計算のための複数の強力なライブラリを単一のフレームワークに統合した sf の能力は、高性能ライブラリを用いた再現可能な地理データ解析の世界への「参入障壁」を低減する顕著な成果である。 sf の機能は、7つの vignette を含むウェブサイト(r-spatial.github.io/sf/)で詳しく説明されている。 これらは、以下のようにオフラインで見ることができる。(訳注:sf などを含めた日本語版の vignette を参照。)

vignette(package = "sf") # どのような vignette があるか表示
vignette("sf1")         # sf の紹介

最初の vignette で説明されているように、R のシンプルフィーチャオブジェクトはデータフレームに格納され、地理データは通常 ‘geom’ または ‘geometry’ という名前の特別な列を占める。 本章の冒頭で読み込んだ spData が提供する world データセットを使って、sf オブジェクトの内容とその動作について説明する。 world は、空間と属性の列を含む ‘sf データフレーム’ で、その名前は関数 names() によって返される(この例の最後の列は地理情報を含んでいる)。

class(world)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
names(world)
#>  [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
#>  [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"

この geom 列の内容は、sf オブジェクトに空間的な力を与える。world$geom は、国別ポリゴンのすべての座標を含む ‘list column’ である。 sf オブジェクトは、R の基本関数である plot() を使ってすぐにプロットすることができる。 次のコマンドは、Figure 2.3 を作成する。

plot(world)
sf パッケージを用いた世界の空間プロットで、各属性ごとのファセットを表示。

FIGURE 2.3: sf パッケージを用いた世界の空間プロットで、各属性ごとのファセットを表示。

多くの GIS プログラムの場合、地理的オブジェクトに対してデフォルトで単一のマップを作成する。一方、sf オブジェクトを plot() すると、データセットの各変数に対してマップが作成されることに注意しておこう。(訳注:このように、属性値を変えた同じ図を並べることを R では「ファセット」(facet) という。ファセットは、元々は切り子という意味。) この動作は、さまざまな変数の空間分布を調べるのに有効であり、Section 2.2.3 で詳しく説明する。

より広く言えば、地理的なオブジェクトを空間的な力を持つ通常のデータフレームとして扱うことは、特にデータフレームの扱いに慣れている場合、多くの利点がある。 例えば、よく使われる summary() 関数では、world オブジェクト内の変数の概要が表示され、便利である。

summary(world["lifeExp"])
#>     lifeExp                geom    
#>  Min.   :50.6   MULTIPOLYGON :177  
#>  1st Qu.:65.0   epsg:4326    :  0  
#>  Median :72.9   +proj=long...:  0  
#>  Mean   :70.9                      
#>  3rd Qu.:76.8                      
#>  Max.   :83.6                      
#>  NA's   :10

summary() コマンドでは、1つの変数しか選択していないが、ジオメトリのレポートも出力される。 これは、sf オブジェクトのジオメトリ列の「スティッキー」な動作を示している。つまり、Section 3.2 で見るように、ユーザーが意図的に削除しない限り、ジオメトリは保持されるということである。 この結果は、world に含まれる非空間的データと空間的データの両方を簡単に要約したものである。すべての国に対して、平均寿命は71歳(51歳未満から83歳以上の範囲、中央値は73歳)である。

上記のサマリー出力にある MULTIPOLYGON という単語は、world オブジェクトに含まれるフィーチャ(国)のジオメトリタイプを表している。 この表現は、インドネシアやギリシャのような島を持つ国には必要である。(訳注:日本も当然ながら当てはまる。) その他のジオメトリタイプについては、Section 2.2.4 で説明する。

このシンプルフィーチャオブジェクトの基本的な動作と内容を深く見てみると、「spatial data frame」と考えるのが妥当であることがわかる。

以下のコードは、world オブジェクトの最初の 2 行と最初の 3 列だけを含む sf オブジェクトを返す方法を示している。 この出力は、通常の data.frame と比較して、2つの大きな違いがある。それは、追加の地理的メタデータ(Geometry typeDimensionBounding box 、および Geodetic CRS という CRS 情報で始まる行の座標参照系情報)が含まれていることと、ここでは geom という名前の「幾何学列」が存在することである。

world_mini = world[1:2, 1:3]
world_mini
#> Simple feature collection with 2 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -18.3 xmax: 180 ymax: -0.95
#> Geodetic CRS:  WGS 84
#> # A tibble: 2 × 4
#>   iso_a2 name_long continent                                                geom
#>   <chr>  <chr>     <chr>                                      <MULTIPOLYGON [°]>
#> 1 FJ     Fiji      Oceania   (((-180 -16.6, -180 -16.5, -180 -16, -180 -16.1, -…
#> 2 TZ     Tanzania  Africa    (((33.9 -0.95, 31.9 -1.03, 30.8 -1.01, 30.4 -1.13,…

特に、「シンプル」であるべきクラスシステムにとっては、かなり複雑に見えるだろう。 しかし、このように物事を整理し、ベクタ地理データセットを扱うのに sf を使うには、それなりの理由がある。

sf パッケージがサポートする各ジオメトリタイプについて説明する前に、sf オブジェクトの構成要素を理解するために基本に立ち返ってみよう。 Section 2.2.5 は、シンプルフィーチャオブジェクトが、特殊なジオメトリ列を持つデータフレームであることを示している。 これらの空間列は、geom または geometry という列名になっている。world$geom は、上記の world オブジェクトの空間要素を指す。 これらのジオメトリ列は、クラス sfc の「リスト列」である(Section 2.2.7 を参照)。 sfc オブジェクトは、sfg(simple feature geometry、Section 2.2.6 で説明)というクラスの1つ以上のオブジェクトから構成されている。

シンプルフィーチャの空間成分の働きを理解するためには、シンプルフィーチャのフィーチャを理解することが必須である。 このため、現在サポートされているシンプルフィーチャのフィーチャタイプについて Section 2.2.4 で説明した後、sfg および sfc オブジェクトをベースとする sf オブジェクトを使用して R でこれらを表現する方法について説明する。

先ほどのコードでは、world_mini = world[1:2, 1:3] というコマンドで = を使って world_mini という新しいオブジェクトを作成している。 これは代入という。 同じ結果を得るために同等のコマンドは world_mini <- world[1:2, 1:3] である。 「矢印代入」の方が一般的だが、「イコール代入」の方が入力しやすく、Python や JavaScript などのよく使われる言語との互換性もあり教えやすいので、ここでは「イコール代入」を使っている。 どちらを使うかは、一貫性がある限り、好みの問題である(styler などのパッケージを使えば、スタイルを変更することが可能)。

2.2.2 なぜシンプルフィーチャなのか?

シンプルフィーチャは、QGIS や PostGIS など、多くの GIS アプリケーションのデータ構造の根幹をなすデータモデルとして広く支持されている。 データモデルを使用することで、例えば空間データベースからのインポートや空間データベースへのエクスポートなど、他のセットアップとの相互移行が可能になることが大きなメリットである。

R の観点からすると、「sp がすでに使われて検証もされているのに、なぜ sf パッケージを使うのか」という質問になる。 理由はいろいろある(シンプルフィーチャモデルの利点と連動している)。

  • 高速なデータの読み書きが可能
  • プロット性能の向上
  • sf オブジェクトは、ほとんどの操作でデータフレームとして扱うことができる
  • sf 関数名は比較的一貫性があり、直感的に理解できる (すべて st_ で始まる)。
  • sf 関数は |> 演算子と組み合わせることができ、R パッケージの tidyverse コレクションとうまく機能する。

sftidyverse パッケージのサポートは、地理ベクタデータセットを読み込むための read_sf() 関数で示される。これは、地理ベクタを読むための関数で、Section 8.6.1 で詳しく解説する。 関数 st_read() は、base R の data.frame に格納された属性を返す (長いメッセージを表示する。以下のコードチャンクには表示していない。)。一方、read_sf() は、データを tidyverse tibble として返し、表示は少ない。 以下、実際の例を紹介する(地理ベクタデータの読み方については Section 8.6.1 を参照)。

world_dfr = st_read(system.file("shapes/world.shp", package = "spData"))
#> Reading layer `world' from data source 
#>   `/Users/baba/Library/R/x86_64/4.2/library/spData/shapes/world.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 177 features and 10 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.6
#> Geodetic CRS:  WGS 84
world_tbl = read_sf(system.file("shapes/world.shp", package = "spData"))
class(world_dfr)
#> [1] "sf"         "data.frame"
class(world_tbl)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"

Chapter 3tidyverse 関数を使って sf オブジェクトを操作する方法を示す通り、sf は今や R での空間ベクタデータの分析に最適なパッケージである(spatstat パッケージは、多くの空間統計用の関数を提供するものの最適とは言えない。)。 spatstatterra は、空間統計のための多くの関数を提供するパッケージ・エコシステムであり、どちらもベクトル地理データクラスを持っているが、ベクタを扱う点においては sf と同じレベルの取り込みはしていない。 多くの人気パッケージが sf をベースに構築されており、前章の Section 1.4 にあるように、1日あたりのダウンロード数でその人気が上昇していることが示されている。 レガシーパッケージである rgeosrgdal から既存のパッケージやワークフローを移行するには時間がかかるが (Bivand 2021)、ロード時に「2023 年末までに引退させる」というメッセージによって表示されるように、できる限り早く移行するべきである。 まだ古いパッケージを使用している人は、「できるだけ早く GDAL と PROJ を使用した sf/stars/terra 関数に移行しよう」と言うことになる。

つまり、sf は将来性があるが、sp には将来性はないのである。 レガシークラスシステムに依存するワークフローでは、sf オブジェクトは、以下のように sp パッケージの Spatial クラスとの間で変換することができる。

library(sp)
world_sp = as(world, "Spatial") # sf オブジェクトを sp へ
# sp 関数 ...
world_sf = st_as_sf(world_sp)          # sp から sf

2.2.3 基本的な地図の作り方

基本的な地図は sfplot() で作成する。 デフォルトでは、Figure 2.4 の左側のパネルに示されているように、オブジェクトの各変数に対して1つのサブプロット、複数パネルのプロットが作成される。 プロットされるオブジェクトが単一の変数である場合、連続した色を持つ凡例または「キー」が生成される(右側のパネルを参照)。 色は col = でも設定できるが、この場合、連続したパレットや凡例は作成されていない。

plot(world[3:6])
plot(world["pop"])
sf を使ったプロット。複数変数(左)と単一変数(右)。sf を使ったプロット。複数変数(左)と単一変数(右)。

FIGURE 2.4: sf を使ったプロット。複数変数(左)と単一変数(右)。

プロットは、add = TRUE を設定することで、既存の画像にレイヤとして追加される。10 このことを示すために、また、属性と空間データ操作に関する Chapter 3 と Chapter 4 の内容を理解するために、次のコードチャンクはアジアの国々をフィルターして、1つのフィーチャに結合している。

world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)

世界地図の上にアジア大陸をプロットすることができるようになった。 add = TRUE が動作するためには、最初のプロットは1つのファセットだけでなければならないことに注意しておこう。 最初のプロットにキーがある場合は、reset = FALSE を使用する必要がある。

plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red")
世界各国の上にレイヤとしてアジアを加えたプロット。

FIGURE 2.5: 世界各国の上にレイヤとしてアジアを加えたプロット。

このようにレイヤを追加していくことで、レイヤ間の地理的な対応関係を検証することができる。 plot() 関数は、実行速度が速く、必要なコード行数も少ないのであるが、幅広いオプションを持つインタラクティブなマップを作成することはできない。 より高度な地図作成には、tmap などの専用可視化パッケージの利用を勧める(Chapter 9 参照)。

sfplot() メソッドでマップを修正する方法はいろいろある。 sf は R の基本的な描画メソッド plot()を拡張しているので、main = (地図のタイトルを指定する)などの引数は sf オブジェクトでも動作する( ?graphics::plot?par を参照)。11

Figure 2.6 は、この柔軟性を、世界地図の上に、直径(cex = で設定)が各国の人口を表す円を重ねることで表現している。 この図の非投影版は、以下のコマンドで作成できる(本章末の練習問題と、スクリプト 02-contplot.R を使って Figure 2.6 を再現することができる。)

plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)
plot(st_geometry(world_cents), add = TRUE, cex = cex)
国別大陸(塗りつぶし色で表現)と2015年の人口(円で表現、面積は人口に比例)

FIGURE 2.6: 国別大陸(塗りつぶし色で表現)と2015年の人口(円で表現、面積は人口に比例)

上記のコードでは、関数 st_centroid() を使って、あるジオメトリタイプ(ポリゴン)を別のタイプ(ポイント)に変換している( Chapter 5 参照)。引数 cex で見た目を変化させることができる。

sf の plot メソッドには、地理データに特有の引数もある。以下の例で、sf オブジェクトをコンテクストでプロットするためにexpandBB を使用してみよう。 expandBB は長さ4の数値ベクトルを取り、プロットの外接枠をゼロから相対的に下、左、上、右の順序で拡張する。 これを利用して、インドを巨大なアジアの隣国(東にある中国に重点を置く)の文脈でプロットすると、次のコードチャンクで Figure 2.7 を生成する。 (プロットにテキストを追加することについては、以下の演習を参照)。12

india = world[world$name_long == "India", ]
plot(st_geometry(india), expandBB = c(0, 0.2, 0.1, 1), col = "gray", lwd = 3)
plot(st_geometry(world_asia), add = TRUE)
インドの文脈で、expandBB論を実証。

FIGURE 2.7: インドの文脈で、expandBB論を実証。

プロットコードでインドを強調するために lwd を使用していることに注意。 さまざまな種類のジオメトリを表現するためのその他の視覚化技術については、Section 9.6 を参照。なお、ジオメトリについては次のセクションで開設する。

2.2.4 ジオメトリの種類

ジオメトリは、シンプルフィーチャを構成する基本的な要素である。 R のシンプルフィーチャは、sf パッケージがサポートする18種類のジオメトリタイプのいずれかを取ることができる。 この章では、最もよく使われる以下の 7 つのタイプに焦点を当てる: POINTLINESTRINGPOLYGONMULTIPOINTMULTILINESTRINGMULTIPOLYGONGEOMETRYCOLLECTIONPostGIS マニュアルで可能なフィーチャの全リストを検索してみよう。

一般に、シンプルフィーチャの符号化方式としては、WKB(Well-known binary)や WKT(Well-known text)が標準的である。 WKB の表現は通常、コンピュータで読みやすい16進数の文字列である。 このため、GIS や空間データベースでは、ジオメトリオブジェクトの転送や保存に WKB を使用している。 一方、WKT は、シンプルフィーチャを人間が読みやすいテキストマークアップで記述したものである。 どちらの形式も交換可能であり、ここでは当然 WKT で表す。

各ジオメトリタイプの基本は点である。 点(point)とは、2次元、3次元、4次元空間(詳しくは vignette("sf1") を参照、訳注:日本語版)の座標で、次のようなものである(Figure 2.8 の左図参照)。

  • POINT (5 2)

線(linestring、polyline)とは、例えば、点と点を結ぶ直線の列のことである(Figure 2.8 の中央を参照)。

  • LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

ポリゴン(polygon)とは、閉じた非交差環を形成する点の並びのことである。 「閉じた」とは、ポリゴンの最初と最後の点が同じ座標を持つことを意味する( Figure 2.8 の右図参照)13

  • 穴のないポリゴン POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
点、線、ポリゴンのジオメトリ説明図。

FIGURE 2.8: 点、線、ポリゴンのジオメトリ説明図。

ここまでは、1つのフィーチャに1つの幾何学的実体しかないジオメトリを作成した。 しかし、sf では、各ジオメトリタイプの「複合」バージョンを使用して、1つのフィーチャに複数のジオメトリを存在させることもできる。

  • 複合点 MULTIPOINT (5 2, 1 3, 3 4, 3 2)
  • 複合線 MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
  • 複合ポリゴン MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5), (0 2, 1 2, 1 3, 0 3, 0 2)))
MULTI*ジオメトリの説明図。

FIGURE 2.9: MULTI*ジオメトリの説明図。

最後に、ジオメトリコレクションには、(複合の)点や線を含むジオメトリの任意の組み合わせを含めることができる(Figure 2.10 を参照)。

  • ジオメトリコレクション GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2), LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))
ジオメトリコレクションの説明図

FIGURE 2.10: ジオメトリコレクションの説明図

2.2.5 sf クラス

シンプルフィーチャは、フィーチャと非フィーチャという2つの部分からなる。 Figure 2.11 は sf オブジェクトの作成方法を示している。ジオメトリは sfc オブジェクトから、属性は data.frame または tibble から取得される。 sf ジオメトリをゼロから構築する方法については、下記の Section 2.2.6 と Section 2.2.7 を参照。

sf オブジェクトの構成。

FIGURE 2.11: sf オブジェクトの構成。

フィーチャ以外の属性は、フィーチャの名称や、測定値、グループなどの属性を表す。 属性を説明するために、2017年6月21日のロンドンの気温が25℃であることを表してみたい。 この例では、ジオメトリ(座標)と、3つの異なるクラスを持つ属性(地名、気温、日付)が含まれている。14 クラス sf のオブジェクトは、属性 (data.frame) とシンプルフィーチャ列 (sfc) を組み合わせて、そのようなデータを表現するものである。 これらは、下図のように st_sf() で、ロンドンの例を作成する。

lnd_point = st_point(c(0.1, 51.5))                # sfg object
lnd_geom = st_sfc(lnd_point, crs = 4326)          # sfc object
lnd_attrib = data.frame(                          # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)   # sf object

何が起こったのだろうか?まず、座標を使ってシンプルフィーチャ(sfg)を作成した。 次に、ジオメトリをシンプルフィーチャ列(sfc)に変換し、CRS を設定した。 第三に、属性は data.frame に格納され、sfc のオブジェクトと st_sf() で結合された。 この結果、以下に示すように、sf オブジェクトが生成される(一部の出力は省略している)。

lnd_sf
#> Simple feature collection with 1 features and 3 fields
#> ...
#>     name temperature       date         geometry
#> 1 London          25 2017-06-21 POINT (0.1 51.5)
class(lnd_sf)
#> [1] "sf"         "data.frame"

その結果、sf のオブジェクトは、実際には sfdata.frame という2つのクラスを持っていることがわかった。 シンプルフィーチャとは、単純にデータフレーム(四角い表)であるが、空間属性がリスト列に格納されている。この列は、通常 Section 2.2.1 で説明されているように、geometry と呼ばれるものである。 この二面性が、「シンプルフィーチャ」のコンセプトの中心である。 ほとんどの場合、sfdata.frame のように扱われ、動作することができる。 シンプルフィーチャとは、要するにデータフレームに空間的な拡張を施したものである。

2.2.6 シンプルフィーチャ・ジオメトリ(sfg)

sfg クラスは、R のさまざまなシンプルフィーチャの種類を表する。点、線、ポリゴン(および複合点などの「複合」対応)、またはジオメトリコレクションである。

通常は、既存の空間ファイルを読み込むだけなので、自分でジオメトリを作成する面倒な作業は必要ない。 しかし、もし必要であれば、シンプルフィーチャオブジェクト ( sfg ) を一から作成するための関数群が用意されている。 これらの関数の名前はシンプルで一貫しており、すべて st_ という接頭辞で始まり、小文字でジオメトリタイプの名前で終わる。

sfg オブジェクトは、3 つの基本的な R データ型から作成することができる。

  1. 数値ベクトル:ひとつの点
  2. 行列:点の集合で、各行が点、多点、または線分を表す。
  3. リスト:行列、マルチ文字列、ジオメトリコレクションなどのオブジェクトの集合体

関数 st_point() は、数値ベクトルから単一の点を作成する。

st_point(c(5, 2))                # XY point
#> POINT (5 2)
st_point(c(5, 2, 3))             # XYZ point
#> POINT Z (5 2 3)
st_point(c(5, 2, 1), dim = "XYM") # XYM point
#> POINT M (5 2 1)
st_point(c(5, 2, 3, 1))          # XYZM point
#> POINT ZM (5 2 3 1)

その結果、長さ 2、3、4 のベクトルからそれぞれ XY(2次元座標)、XYZ(3次元座標)、XYZM(3次元に追加変数、通常は測定精度)の点タイプが生成されることがわかった。 XYM タイプは、dim 引数(dimension の略)で指定する必要がある。

これに対して、複合点(st_multipoint())と複合線(st_linestring())のオブジェクトの場合は、行列 (matrix) を使用する。

# rbind 関数により行列の作成が簡単になった。
## MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
st_multipoint(multipoint_matrix)
#> MULTIPOINT ((5 2), (1 3), (3 4), (3 2))
## LINESTRING
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)
#> LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

最後に、複合線、(複合)ポリゴン、ジオメトリコレクションの作成には、リストを使用する。

## POLYGON
polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
## 穴あきポリゴン
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5), (2 4, 3 4, 3 3, 2 3, 2 4))
## MULTILINESTRING
multilinestring_list = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
st_multilinestring((multilinestring_list))
#> MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
## MULTIPOLYGON
multipolygon_list = list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
                         list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
st_multipolygon(multipolygon_list)
#> MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5)), ((0 2, 1 2, 1 3, 0 3, 0 2)))
## GEOMETRYCOLLECTION
geometrycollection_list = list(st_multipoint(multipoint_matrix),
                              st_linestring(linestring_matrix))
st_geometrycollection(geometrycollection_list)
#> GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2),
#>   LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))

2.2.7 シンプルフィーチャ列 (sfc)

1つの sfg オブジェクトには、1つのシンプルフィーチャだけが含まれている。 シンプルフィーチャ列(simple feature column, sfc)は、sfg オブジェクトのリストであり、さらに、使用中の座標参照系に関する情報を含めることができる。 例えば、単純な2つのフィーチャを2つのフィーチャを持つ1つの物体に結合するには、st_sfc() 関数を使用することができる。 sfcsf データフレームのジオメトリ列を表すので、これは重要である。

# sfc POINT
point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc
#> Geometry set for 2 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
#> CRS:           NA
#> POINT (5 2)
#> POINT (1 3)

ほとんどの場合、sfc オブジェクトは、同じジオメトリタイプのオブジェクトを含んでいる。 したがって、ポリゴン型の sfg オブジェクトをシンプルフィーチャ列に変換すると、ポリゴン型の sfc オブジェクトもできてしまい、st_geometry_type() で確認することができる。 同様に、multilinestring のジオメトリ列は、multilinestring 型の sfc オブジェクトになる。

# sfc POLYGON
polygon_list1 = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
polygon1 = st_polygon(polygon_list1)
polygon_list2 = list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
polygon2 = st_polygon(polygon_list2)
polygon_sfc = st_sfc(polygon1, polygon2)
st_geometry_type(polygon_sfc)
#> [1] POLYGON POLYGON
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
# sfc MULTILINESTRING
multilinestring_list1 = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
multilinestring1 = st_multilinestring((multilinestring_list1))
multilinestring_list2 = list(rbind(c(2, 9), c(7, 9), c(5, 6), c(4, 7), c(2, 7)), 
                            rbind(c(1, 7), c(3, 8)))
multilinestring2 = st_multilinestring((multilinestring_list2))
multilinestring_sfc = st_sfc(multilinestring1, multilinestring2)
st_geometry_type(multilinestring_sfc)
#> [1] MULTILINESTRING MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

また、異なるジオメトリタイプの sfg オブジェクトから、sfc オブジェクトを作成することも可能である。

# sfc GEOMETRY
point_multilinestring_sfc = st_sfc(point1, multilinestring1)
st_geometry_type(point_multilinestring_sfc)
#> [1] POINT           MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

前述したように、sfc のオブジェクトは、さらに座標参照系(CRS)の情報を格納することができる。 デフォルト値は、NA (Not Available) である。st_crs() で確認できる。

st_crs(points_sfc)
#> Coordinate Reference System: NA

sfc オブジェクトのすべてのジオメトリは、同じ CRS を持つ必要がある。 CRS は st_sfc() (または st_sf() ) の crs 引数で指定できる。これは crs = "EPSG:4326" のようなテキスト文字列として提供される CRS 識別子を取る (他の CRS 表現とこの意味の詳細については Section 7.2 を参照)。

# 'EPSG' CRS code を参照する識別子で CRS を設定
points_sfc_wgs = st_sfc(point1, point2, crs = "EPSG:4326")
st_crs(points_sfc_wgs) # print CRS (only first 4 lines of output shown)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#> ...

2.2.8 sfheaders パッケージ

sfheaders は、sf オブジェクトの構築、変換、操作を高速化する R パッケージである (Cooley 2020)。 これは、ベクトル、行列、データフレームから sf オブジェクトを、 sf ライブラリに依存せずに高速に構築することに重点を置いており、 ヘッダファイルを通じてその基礎となる C++ コードを公開する (sfheaders という名前はこのためである。) 。 このアプローチにより、他の人がコンパイルされた高速に動作するコードを使って拡張することが可能になる。 sfheaders のコア関数すべてに対して、対応する C++ の実装があり、その説明は Cpp vignette で説明されている。 多くの人にとって、R の関数は、パッケージの計算速度の恩恵を受けるのに十分すぎるほどだろう。 sfheaderssf とは別に開発されたが、完全に互換性があり、前のセクションで説明したタイプの有効な sf オブジェクトを作成することを目的としている。

sfheaders の最も単純な使用例は、sfgsfcsf オブジェクトの構築例とともに、以下のコードチャンクで示されている。

  • sfg_POINT に変換されたベクトル
  • sfg_LINESTRING に変換された行列
  • sfg_POLYGON に変換されたデータフレーム

まず、最も単純な sfg オブジェクトを作成する。これは、v というベクタに割り当てられた、1 つの座標ペアである。

v = c(1, 1)
v_sfg_sfh = sfheaders::sfg_point(obj = v)
v_sfg_sfh # sf をロードせずに表示
#>      [,1] [,2]
#> [1,]    1    1
#> attr(,"class")
#> [1] "XY"    "POINT" "sfg" 

上の例では、sf が読み込まれていないときに sfg オブジェクト v_sfg_sfh が出力され、その基本的な構造を示している。 sf が読み込まれた場合(今回のように)、上記のコマンドの結果は sf のオブジェクトと区別がつかない。

v_sfg_sf = st_point(v)
print(v_sfg_sf) == print(v_sfg_sfh)
#> POINT (1 1)
#> POINT (1 1)
#> [1] TRUE

次の例は、sfheaders が行列とデータフレームから sfg オブジェクトを作成する方法を示している。

# matrices
m = matrix(1:8, ncol = 2)
sfheaders::sfg_linestring(obj = m)
#> LINESTRING (1 5, 2 6, 3 7, 4 8)
# data.frames
df = data.frame(x = 1:4, y = 4:1)
sfheaders::sfg_polygon(obj = df)
#> POLYGON ((1 4, 2 3, 3 2, 4 1, 1 4))

オブジェクト vmdf を再利用して、以下のように簡単な特徴列 (sfc) を作ることもできる(出力は示さない)。

sfheaders::sfc_point(obj = v)
sfheaders::sfc_linestring(obj = m)
sfheaders::sfc_polygon(obj = df)

同様に、sf オブジェクトも以下のように作成することができる。

sfheaders::sf_point(obj = v)
sfheaders::sf_linestring(obj = m)
sfheaders::sf_polygon(obj = df)

いずれの例も CRS(座標参照系)は定義されていない。 もし、sf 関数を使った計算や幾何演算をする予定があるのなら、CRS を設定することを勧める(詳しくは、Chapter 7 参照)。

df_sf = sfheaders::sf_polygon(obj = df)
st_crs(df_sf) = "EPSG:4326"

また、sfheaderssf オブジェクトの「分解」と「再構築」を得意としている。つまり、ジオメトリ列を、各頂点の座標とジオメトリフィーチャ(および複数フィーチャ)の ID に関するデータを含むデータフレームに変換するのである。 Chapter 5 で取り上げた、異なるタイプのジオメトリ列を「キャスト」する際に、高速かつ信頼性の高い処理を行うことができる。 パッケージの documentation やこの本のために開発されたテストコードでのベンチマークでは、このような操作に対して sf パッケージよりもはるかに高速であることが示されている。

2.2.9 S2による球面幾何学操作

球面幾何学エンジンは、世界が丸いという事実に基づいているが、2点間の直線やポリゴンに囲まれた面積の計算など、地盤計算のための簡単な数学的手続きは、平面(投影)幾何学を前提としている。 sf バージョン 1.0.0 以降、R は Google の S2 球面幾何エンジンとのインターフェース s2 インターフェースパッケージにより、球面幾何の操作を「すぐに」サポートする。 S2 は、離散型グローバルグリッドシステム (Discrete Global Grid System, DGGS) の例として、おそらく最もよく知られている。 もう一つの例は、H3 グローバルな六角形の階層的な空間インデックスである (Bondaruk, Roberts, and Robertson 2020)

S2 には、e66ef376f790adf8a5af7fca9e6e422c03c9143f のような文字列を使って地球上の任意の場所を記述するのに役立つ場合もあるかもしれない。むしろ、S2 に対する sf のインターフェースの主な利点は、距離、バッファ、面積計算などの計算のためのドロップイン関数を提供している点である。これは sf のドキュメントで説明されており、コマンド vignette("sf7")で開くことができる(訳注:日本語版)。

sf は、S2 に対して、オンとオフの2つのモードで動作させることができる。 S2 ジオメトリエンジンはデフォルトでオンになっている。以下のコマンドで確認することができる。

sf_use_s2()
#> [1] TRUE

ジオメトリエンジンをオフにした結果の例を以下に示す。この章で作成した india オブジェクトの周りにバッファを作成する(S2 をオフにしたときに発せられる警告に注意)。

india_buffer_with_s2 = st_buffer(india, 1)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
india_buffer_without_s2 = st_buffer(india, 1)
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
S2 ジオメトリエンジンをオフにした場合の結果の例。インド周辺のバッファを同じコマンドで作成したが、紫色のポリゴンオブジェクトは S2 をオンにして作成したため、バッファは 1 m になった。大きい薄緑色のポリゴンは S2 をオフにして作成したため、バッファは不正確な経度/緯度の単位となった。

FIGURE 2.12: S2 ジオメトリエンジンをオフにした場合の結果の例。インド周辺のバッファを同じコマンドで作成したが、紫色のポリゴンオブジェクトは S2 をオンにして作成したため、バッファは 1 m になった。大きい薄緑色のポリゴンは S2 をオフにして作成したため、バッファは不正確な経度/緯度の単位となった。

本書では、特に明記していない限り、S2 がオンになっていることを前提に説明する。 以下のコマンドで再度オンにする。

sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
sf の S2 使用は多くの場合妥当だが、場合によっては、R セッションの間、あるいはプロジェクト全体において S2 をオフにする正当な理由があることもある。 sf の GitHub リポジトリの issue 1771 に書かれているように、デフォルトの動作は、S2をオフにした場合(そして古いバージョンの sf で)動作するコードを失敗させる可能性がある。 こうした極端な事例には、S2 の厳密な定義に従っていないポリゴンに対する操作が含まれる。 もし、「#> Error in s2_geography_from_wkb ...」のようなエラーメッセージが表示されたら、S2 をオフにした後に、エラーメッセージを発生させたコマンドを再度試してみる価値があるかもしれない。 プロジェクト全体で S2 を無効にするには、プロジェクトのルートディレクトリ(メインフォルダ)に .Rprofile というファイルを作成し、その中に sf::sf_use_s2(FALSE) というコマンドを記述すればよい。

2.3 ラスタデータ

空間ラスタデータモデルは、連続した格子状のセル(ピクセルとも呼ばれる; Figure 2.13 :A)で世界を表現する。 このデータモデルでは、各セルが同じ一定の大きさを持つ、いわゆる正方形のグリッドを指すことが多く、本書では正方形のグリッドにのみ焦点を当てることにする。 しかし、他にも回転格子、せん断格子、直線格子、曲線格子など、いくつかの種類の格子が存在する(Pebesma and Bivand (2022) の第1章、または Tennekes and Nowosad (2022) の第2章を参照のこと)。

ラスタデータモデルは通常、ラスタヘッダと、等間隔に並んだセル(ピクセルとも呼ばれる; Figure 2.13 :A)を表す行列(行と列を持つ)から構成される。15 ラスタヘッダは、座標参照系、範囲、原点を定義する。 原点は、多くの場合、行列の左下隅の座標である(ただし、terra パッケージでは、デフォルトで左上隅が使用される(Figure 2.13:B))。 ヘッダは、列数、行数、セルサイズの解像度によって範囲を定義する。 したがって、原点から出発して、セルの ID を使うか(Figure 2.13:B)、明示的に行と列を指定することで、1つのセルに対して簡単にアクセスし、変更することができる。 この行列表現により、直方体のベクタポリゴンのように、各セルの角の4点の座標を明示的に保存する必要がなくなる(実際には、原点という1つの座標しか保存しない)。 これと地図代数(Section 4.3.2)により、ラスタ処理はベクタデータ処理よりもはるかに効率的で高速に処理することができる。 しかし、ベクタデータとは異なり、ラスタレイヤの1つのセルには1つの値しか入れることができない。 値は、数値またはカテゴリ(Figure 2.13:C)である。

ラスタデータの種類:(A)セルID、(B)セル値、(C)色付きラスタマップ。

FIGURE 2.13: ラスタデータの種類:(A)セルID、(B)セル値、(C)色付きラスタマップ。

ラスタ地図は通常、標高、気温、人口密度、スペクトルデータなどの連続的な現象を表現する。 土壌や土地被覆クラスなどの離散的なフィーチャも、ラスタデータモデルで表現することができる。 ラスタデータセットの両方の使い方を説明するために、Figure 2.14、ラスタデータセットでは、離散的なフィーチャの境界が不鮮明になることがあることを示している。 アプリケーションの性質によっては、離散的なフィーチャのベクタ表現がより適切な場合もある。

連続ラスタとカテゴリラスタの例。

FIGURE 2.14: 連続ラスタとカテゴリラスタの例。

2.3.1 ラスタデータを扱うための R パッケージ

過去 20 年の間に、ラスタデータセットを読み込んで処理するためのパッケージがいくつか開発された。 Section 1.5 にあるように、その中でも特に raster は 2010 年にリリースされ、R のラスタ機能を一変させ、terrastars が開発されるまでこの分野の最高峰のパッケージとなった。 最近開発されたこの2つのパッケージは、ラスタデータを扱うための強力で高性能な機能を備えており、両者の使用例にはかなりの重複がある。 この本では、古くて(多くの場合)遅い raster に代わる terra に焦点を当てる。 terra のクラスシステムについて学ぶ前に、このセクションでは terrastars の類似点と相違点について説明する。この知識は、様々な状況下でどちらが最も適切かを判断するのに役立つ。

まず、terra は最も一般的なラスタデータモデル(通常のグリッド)に焦点を当て、stars はあまり一般的ではないモデル(通常のグリッド、回転グリッド、シアーグリッド、レクチリニアグリッド、カーヴィリニアグリッドなど)も保存できるようになっている。 通常、terra は1層または多層のラスタを扱うが16stars パッケージはラスタデータキューブを保存する方法を提供する。一方、stars パッケージは、ラスタデータキューブ(多くのレイヤ(バンドなど)、多くの時間(月など)、多くの属性(センサータイプ A とセンサータイプ B など)を持つラスタオブジェクト)を保存する方法を提供する。 重要なのは、どちらのパッケージでも、データキューブのすべてのレイヤまたは要素が同じ空間寸法および範囲を持っている必要があることである。 第二に、どちらのパッケージもラスタデータをすべてメモリに読み込むか、メタデータだけを読み込むかを選択できる。これは通常、入力ファイルのサイズに応じて自動的に行われる。 しかし、両者はラスタ値の保存方法が大きく異なる。 terra はC++のコードをベースにしており、ほとんど C++ のポインタを使用している。 stars は、小さいラスタでは配列のリストとして、大きいラスタでは単なるファイルパスとして値を格納する。 第三に、stars の関数は sf のベクタオブジェクトや関数と密接に関係しており、一方 terra はベクタデータに対して独自のオブジェクトクラス、すなわち SpatVector を使用している。 第四に、両パッケージは、そのオブジェクトに対して様々な機能がどのように働くかについて、異なるアプローチを持っている。 terra パッケージは、ほとんどが多数の組み込み関数に依存しており、各関数は特定の目的(例えば、リサンプリングやトリミングなど)を持っている。 一方、stars はいくつかの組み込み関数(通常 st_ で始まる名前)を使用し、既存のR関数(例えば split()aggregate() )に対する独自のメソッドを持ち、また既存の dplyr 関数(例えば filter()slice() )に対するメソッドも持つ。

オブジェクトを terra から stars に変換する(st_as_stars() を使用)、またはその逆(rast() を使用)が簡単である点が重要である。 また、stars パッケージの最も包括的な紹介として、Pebesma and Bivand (2022) を読むことを勧める。

2.3.2 terra 入門

terra パッケージ は、R のラスタオブジェクトをサポートする。 前身の raster(同じ開発者 Robert Hijmans による)と同様、ラスタデータセットの作成、読み込み、書き出し、操作、処理を行うための豊富な関数群を提供する。 terra の機能は、より成熟した raster パッケージとほぼ同じであるが、いくつかの相違点がある。terra の関数は通常、raster の対応する関数よりも計算効率が高い。 一方、raster クラスシステムは人気があり、他の多くのパッケージで使用されている。 古いスクリプトやパッケージとの後方互換性を確保するために、2種類のオブジェクトをシームレスに変換することができる。例えば、raster パッケージの raster()stack()brick() などの関数がある(地理データを扱うための R パッケージの進化については、前の章を参照)。

ラスタデータを操作するための関数に加え、terra はラスタデータセットを扱うための新しいツールを開発するための基盤となる低レベルの関数を多数提供している。 また、terra では、メインメモリに収まりきらないような大きなラスタデータセットを扱うことができる。 この場合、terra はラスタを小さなチャンクに分割し、ラスタファイル全体を RAM にロードする代わりに、それらを繰り返し処理することが可能である。

terra の概念を説明するために、spDataLarge のデータセットを使ってみよう。 これは、ザイオン国立公園(米国ユタ州)のエリアをカバーする数個のラスタオブジェクトと1個のベクタオブジェクトから構成されている。 例えば、srtm.tif はこの地域のデジタル標高モデルである(詳しくは、そのドキュメント ?srtm を参照)。 まず、my_rast という名前の SpatRaster オブジェクトを作成しよう。

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
class(my_rast)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"

コンソールにラスタの名前を入力すると、ラスタヘッダ(寸法、解像度、範囲、CRS)といくつかの追加情報(クラス、データソース、ラスタ値の要約)が出力される。

my_rast
#> class       : SpatRaster 
#> dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
#> resolution  : 0.000833, 0.000833  (x, y)
#> extent      : -113, -113, 37.1, 37.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : srtm.tif 
#> name        : srtm 
#> min value   : 1024 
#> max value   : 2892

dim() は行、列、レイヤの数、ncell() はセル(ピクセル)の数、res() は空間解像度、ext() は空間範囲、crs() は座標参照系(ラスタ再投影は Section 7.8 で扱っている)をそれぞれ報告する専用の関数である。 inMemory() は、ラスタデータがメモリに格納されているか、ディスクに格納されているかを報告する。

help("terra-package") は、terra 関数の完全なリストを返す。

2.3.3 基本的な地図の作り方

sf パッケージと同様に、terra も独自のクラスに対して plot() メソッドを提供する。

plot(my_rast)
異本的なラスタのプロット。

FIGURE 2.15: 異本的なラスタのプロット。

R でラスタデータをプロットすることはこのセクションの範囲外であるが、他にもいくつかのアプローチがある。

  • plotRGB() terra パッケージの関数で、SpatRaster オブジェクトの 3 つのレイヤを基に 赤-緑-青 のプロットを作成する。
  • ラスタおよびベクタオブジェクトの静的および対話型マップを作成するための tmap などのパッケージ(Chapter 9 を参照)
  • 関数、例えば rasterVis パッケージの levelplot() は、経時変化を視覚化する一般的な手法であるファセットを作成するためのものである。

2.3.4 ラスタクラス

SpatRaster クラスは、terra のラスタオブジェクトを表する。 R でラスタオブジェクトを作成する最も簡単な方法は、ディスクまたはサーバーからラスタファイルを読み込むことである( Section 8.6.2 .

single_raster_file = system.file("raster/srtm.tif", package = "spDataLarge")
single_rast = rast(raster_filepath)

terra パッケージは、GDAL ライブラリの助けを借りて、多数のドライバをサポートしている。 ファイルからのラスタは、通常、ヘッダとファイル自体へのポインタを除いて、完全に RAM に読み込まれるわけではない。

ラスタは、同じ rast() 関数を使用して、ゼロから作成することもできる。 次のコードチャンクで、新しい SpatRaster オブジェクトが生成してみよう。 出来上がったラスタは、本初子午線と赤道を中心とした36個のセル(nrowsncols で指定された6列と6行)から構成されている(xminxmaxyminymax パラメータを参照)。 ラスタオブジェクトのデフォルトの CRS は WGS84 であるが、crs 引数で変更することができる。 これは、解像度の単位が度であることを意味し、ここでは 0.5(resolution)に設定している。 各セルには値(vals)が、1 はセル 1、2 はセル 2、といったように割り当てられている。 rast() は、(matrix() とは異なり)左上から順番にセルを埋めていく。つまり、一番上の行には 1 から 6、二番目の行には 7 から 12 の値が含まれる。

new_raster = rast(nrows = 6, ncols = 6, resolution = 0.5, 
                  xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
                  vals = 1:36)

ラスタオブジェクトを作成する他の方法については、?rast を参照。

SpatRaster クラスは、複数のレイヤを処理することもできる。これは通常、1 つのマルチスペクトル衛星ファイルまたは時系列のラスタに対応する。

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)
multi_rast
#> class       : SpatRaster 
#> dimensions  : 1428, 1128, 4  (nrow, ncol, nlyr)
#> resolution  : 30, 30  (x, y)
#> extent      : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 12N (EPSG:32612) 
#> source      : landsat.tif 
#> names       : landsat_1, landsat_2, landsat_3, landsat_4 
#> min values  :      7550,      6404,      5678,      5252 
#> max values  :     19071,     22051,     25780,     31961

nlyr() は、SpatRaster オブジェクトに格納されているレイヤの数を取得する。

nlyr(multi_rast)
#> [1] 4

複数レイヤのラスタオブジェクトの場合、terra::subset() でレイヤを選択することができる。17 第2引数にレイヤ番号またはレイヤ名を指定する。

multi_rast3 = subset(multi_rast, 3)
multi_rast4 = subset(multi_rast, "landsat_4")

逆の操作として、複数の SpatRaster オブジェクトを 1 つにまとめることも c 関数で可能である。

multi_rast34 = c(multi_rast3, multi_rast4)
ほとんどの SpatRaster オブジェクトはラスタの値を保存せず、ファイル自体へのポインタを保存する。 これは重大な副作用がある。つまり、".rds"".rda" ファイルに直接保存したり、クラスタコンピューティングで使用したりすることができない。 (1) wrap() 関数を用いて、Rオブジェクトとして保存したりクラスタコンピューティングで利用できる特殊な一時的オブジェクトを作成する、または (2) writeRaster() を用いて通常のラスタとして保存する。

2.4 地理的および投影された座標参照系

ベクタとラスタの空間データ型は、空間データに内在する概念を共有している。 このうち最も基本的なものは座標参照系(CRS)で、データの空間要素が地球(または他の天体)の表面とどのように関連しているかを定義するものである。 CRSには、本章の冒頭で紹介したように、地理的なものと投影的なものがある( Figure 2.1 参照)。 このセクションでは、各タイプについて説明し、CRSの設定、変換、クエリについて深く掘り下げた Chapter 7 の基礎を築く。

2.4.1 地理座標系

地理座標系は、地球上の任意の位置を経度と緯度という2つの値で特定する(Figure 2.17 の左図参照)。 経度(longitude)は、本初子午線面からの角距離で東西方向の位置である。 緯度(latitude)は赤道面の北または南の角度距離である。 そのため、地理的な CRS の距離はメートル単位ではない。 このことは、Section 7 で示されているように、重要な結果をもたらす。

地理座標系における地球の表面は、球面または楕円面で表現される。 球体モデルは、地球がある半径の完全な球体であると仮定したもので、単純であるという利点があるが、同時に不正確である。 楕円体型モデルは、赤道半径と極半径の2つのパラメータで定義される。 地球は圧縮されているので、これが適している。すなわち、赤道半径は極半径より約 11.5 km 長い (Maling 1992)18

楕円体は、CRS のより広い構成要素である測地系(datum)の一部である。 どの楕円体を使うか、直交座標と地表の位置の正確な関係などの情報が含まれている。 測地系には、地心座標系(WGS84 など)とローカル測地系(NAD83 など)がある。 この2種類の測地系の例は、Figure 2.16 で見ることができる。 黒線は地心測地系を表し、その中心は地球の重心に位置し、特定の場所に最適化されていない。 紫色の破線で示されるローカル測地系は、楕円面を特定の位置の地表に合わせるようにずらしたものである。 これにより、例えば大きな山脈による地表の局所的な変動を、ローカル CRS で説明することができる。 これは Figure 2.16 に見られるように、フィリピンの地域にはローカル測地系が適合しているが、地球表面の他の大部分とはずれているのである。 Figure 2.16 の二つの測地系は、ジオイド(地球平均海面のモデル)の上に載せられている。19

ジオイドの上に表示された地心座標系およびローカル測地系データ(フォールスカラーと、スケールファクター1万による垂直方向の誇張)。ジオイドの画像は Ince et al. (2019) の作品から流用したものである。

FIGURE 2.16: ジオイドの上に表示された地心座標系およびローカル測地系データ(フォールスカラーと、スケールファクター1万による垂直方向の誇張)。ジオイドの画像は Ince et al. (2019) の作品から流用したものである。

2.4.2 投影座標参照系

すべての投影 CRS は、前節で説明した地理的 CRS に基づいており、地図投影に依存して、地球の 3 次元表面を投影 CRS の東方向距離(easting)と北方向距離(Northing)(x と y)の値に変換している。 投影 CRS は、暗黙のうちに平坦な表面上のデカルト座標に基づいている(Figure 2.17 の右パネル)。 投影 CRS には、原点、X軸とY軸、メートルなど線形単位がある。

この推移は、何らかのデフォルメを加えないとできない。 そのため、地表の面積、方向、距離、形状など、地表の一部の性質が歪んでしまう。 投影座標系は、これらの特性のうち1つか2つしか保持できない。 等面積投影は面積を、方位投影は方向を、等距離投影は距離を、コンフォーマ投影は局所形状を保持する。

投影の種類は、円錐型、円筒型、平面型(方位角型)の3つに大別される。 円錐図法では、地球表面が1本の接線または2本の接線に沿って円錐に投影される。 この投影では接線に沿って歪みが最小になり、接線からの距離に応じて歪みが大きくなる。 そのため、中緯度地域の地図に最も適している。 円筒投影は、表面を円筒に写すものである。 この投影は、1本の接線または2本の接線に沿って地表に接することによっても作成することができる。 円筒形の投影は、全世界の地図を作成するときに最もよく使われる。 平面投影は、地球上のある点または接線に沿った平面上にデータを投影するものである。 極域の地図作成によく使われる。 sf_proj_info(type = "proj") は、PROJ ライブラリでサポートされている利用可能な投影のリストを提供する。

さまざまな投影法、その種類、特性、適性についての簡単な要約は、“Map Projections” (1993)https://www.geo-projections.com/ で見ることができる。 Chapter 7 では、CRS について説明し、ある CRS から別の CRS に投影する方法について拡張する。 今は、以下を知っておけば十分である。

  • 座標系は地理的なオブジェクトの重要な構成要素であること
  • データがどの CRS であるか、また、地理的(lon/lat)か投影(通常メートル)かを知ることは重要であり、R による空間およびジオメトリ操作の処理方法に影響を及ぼす。
  • sf オブジェクトの CRS は関数 st_crs() で問い合わせることができ、terra オブジェクトの CRS は関数 crs() を用いて問い合わせることができる。
ベクタデータ型の地理座標系(WGS 84、左)と投影座標系(NAD83 / UTMゾーン12N、右)の例。

FIGURE 2.17: ベクタデータ型の地理座標系(WGS 84、左)と投影座標系(NAD83 / UTMゾーン12N、右)の例。

2.5 単位

CRS の重要な特徴は、空間単位の情報を含んでいることである。 家の寸法がフィートなのかメートルなのかが大事なように、地図でも同じことが言える。 ページや画面上の距離と地上での距離の関係を示すために、地図上にスケールバーやその他の距離表示を追加することは、地図製作の良い習慣である。 同様に、ジオメトリデータやセルを測定する単位を正式に指定して文脈を示し、その後の計算が文脈に沿って行われるようにすることが重要である。

新しい特徴として、sf オブジェクトのジオメトリデータは、単位をネイティブにサポートしている。 これは、sf における距離、面積、その他の幾何学的計算が、units パッケージ (Pebesma, Mailund, and Hiebert 2016) で定義された units 属性に付属する値を返すことを意味している。 これは、単位の違いによる混乱を防ぎ(ほとんどの CRS はメートル、一部の CRS はフィート)、次元に関する情報を提供するのに有利である。 ルクセンブルクの面積を計算する以下のコードで、その様子を見てみよう。

luxembourg = world[world$name_long == "Luxembourg", ]
st_area(luxembourg) # requires the s2 package in recent versions of sf
#> 2.41e+09 [m^2]

出力は平方メートル(m2)の単位で、結果が2次元空間を表していることがわかる。 この情報は属性として保存され(興味のある読者は attributes(st_area(luxembourg)) で確認できる)、人口密度(単位面積当たりの人口、通常は 1 km2)など単位を使用する後の計算に反映させることができる。 単位を報告することで混乱を防ぐことができる。 luxembourg の例で言えば、単位が明記されていないままだと、単位がヘクタールであると誤って判断してしまう可能性がある。 この膨大な数字を消化しやすいサイズに変換するために、結果を100万(1平方キロメートルの中の平方メートル数)で割れば良いと思うかもしれない。

st_area(luxembourg) / 1000000
#> 2409 [m^2]

しかし、その結果を再び平方メートルとするのは誤りである。 解決策としては、units パッケージで正しい単位を設定することである。

units::set_units(st_area(luxembourg), km^2)
#> 2409 [km^2]

ラスタデータの場合、単位も同様に重要である。 しかし、今のところ単位をサポートする空間パッケージは sf だけであり、ラスタデータを扱う場合は分析単位の変更(例えば、ピクセル幅を帝国単位から10進数に変換するなど)に注意して取り組む必要がある。 my_rast オブジェクト(上記参照)は、単位に十進角を持つ WGS84 投影を使用している。 その結果、その分解能も10進数で示されるが、res() 関数は単に数値ベクタを返すだけなので、このことを知っている必要がある。

res(my_rast)
#> [1] 0.000833 0.000833

UTM 投影を使うと、単位が変わってしまう。

repr = project(my_rast, "EPSG:26912")
res(repr)
#> [1] 83.5 83.5

ここでも、res() コマンドは単位を持たない数値ベクタを返すので、UTM 投影の単位がメートルであることを知っておく必要がある。

2.6 演習

E1. Use summary() on the geometry column of the world data object that is included in the spData package. What does the output tell us about:

  • Its geometry type?
  • The number of countries?
  • Its coordinate reference system (CRS)?

E2. Run the code that ‘generated’ the map of the world in Section 2.2.3 (Basic map making). Find two similarities and two differences between the image on your computer and that in the book.

  • What does the cex argument do (see ?plot)?
  • Why was cex set to the sqrt(world$pop) / 10000?
  • Bonus: experiment with different ways to visualize the global population.

E3. Use plot() to create maps of Nigeria in context (see Section 2.2.3).

  • Adjust the lwd, col and expandBB arguments of plot().
  • Challenge: read the documentation of text() and annotate the map.

E4. Create an empty SpatRaster object called my_raster with 10 columns and 10 rows. Assign random values between 0 and 10 to the new raster and plot it.

E5. Read-in the raster/nlcd.tif file from the spDataLarge package. What kind of information can you get about the properties of this file?

E6. Check the CRS of the raster/nlcd.tif file from the spDataLarge package. What kind of information you can learn from it?