8 地理データI/O

必須パッケージ

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

8.1 イントロダクション

この章では、地理データの読み書きの方法について説明する。 地理データ入力 (Import) はジオコンピューテーションに不可欠である。 実世界のアプリケーションはデータなしには不可能である。 データ出力 (Output) も重要で、あなたが研究の結果得られた価値ある新しいデータセットや改良されたデータセットを他の人が利用できるようにすることができる。 これらの入力/出力の処理をまとめて、データ I/O と呼ぶことができる。

地理データの入出力は、プロジェクトの最初と最後に数行のコードで行われることが多い。 簡単なワンステップであるため、見落とされがちである。 しかし、プロジェクトの初期に犯したミス(例えば、古いデータや何らかの欠陥のあるデータセットを使用すること)は、後々大きな問題につながる可能性があるため、どのデータセットが利用可能か、どこで見つけることができるか、どのように取得するのかを確認するためにかなりの時間をかける価値がある。 Section 8.2 では、このトピックについて合計で何テラバイトものデータを含む様々なジオポータルとその使用方法について説明する。 さらにデータへのアクセスを容易にするため、地理データをダウンロードするためのパッケージが多数開発されている。 これらのパッケージについては、Section 8.3 で説明する。

地理学のファイル形式は数多くあり、それぞれに長所と短所がある。 これらファイル形式については、Section 8.5 で説明する。 様々なファイル形式を実際に効率よく読み書きするための処理については、それぞれ Section 8.6、Section 8.7 で説明する。 最後の Section 8.8 では、ビジュアライゼーションに関する Chapter 9 に備えて、ビジュアル出力(地図)を保存するための方法を紹介する。

8.2 オープンデータの取得

インターネット上には膨大かつ増え続ける地理データがあり、その多くは無料でアクセス・利用することができる(ただし、提供者のクレジットを適切に表示することが必要)。 同じデータセットにアクセスする場所が複数あるという意味で、ある意味、データは多すぎるくらいにある。 一部のデータセットは品質が低い。 そこで、最初に最も重要な情報源をいくつか紹介する。 様々な「ジオポータル」(地理空間データセットを提供するウェブサービス、 Data.gov など)は、幅広いデータを提供しているが、特定の場所についてのみ提供している場合が多い(この話題については、最新の Wikipedia page で説明されている)。

グローバルなジオポータルの中には、この問題を克服しているものもある。 例えば、GEOSS portalCopernicus Open Access Hub には、全世界をカバーする多くのラスタデータセットが含まれている。 米国航空宇宙局(NASA)が運営するポータルサイト SEDAC や欧州連合の INSPIRE geoportal から、豊富なベクタデータセットにアクセスすることができ、世界や地域を網羅したデータを入手することができる。

ジオポータルは、ほとんどの場合空間的および時間的範囲などの特性に基づいてデータセットを照会できるグラフィカルなインターフェースを提供している。米国地質調査所の EarthExplorer はその代表例である。 ブラウザ上でインタラクティブにデータセットを探索することは、利用可能なレイヤーを理解する上で効果的な方法である。 しかし、データのダウンロードは、再現性と効率性の観点から、コードで行うのがベストである。 ダウンロードは、主に URL や API を経由して、様々な手法でコマンドラインから開始することができる(例: Sentinel API を参照)。 静的 URL にホストされているファイルは、download.file() でダウンロードすることができる。以下のコードは、 catalog.data.gov/dataset/national-parks から米国の国立公園のデータにアクセスする例である。

download.file(url = "https://irma.nps.gov/DataStore/DownloadFile/673366",
              destfile = "nps_boundary.zip",
              mode = "wb")
unzip(zipfile = "nps_boundary.zip")
usa_parks = read_sf(dsn = "nps_boundary.shp")

8.3 地理データパッケージ

地理データにアクセスするための R パッケージが多数開発されており、その一部を Table 8.1 で紹介している。 これらのパッケージは、1つまたは複数の空間ライブラリやジオポータルへのインターフェースを提供し、コマンドラインからのデータアクセスをさらに高速化することを目的としている。

TABLE 8.1: 地理データ取得の代表的 R パッケージ
パッケージ 説明
osmdata OpenStreetMap の小さなデータセットをダウンロードし、インポート。
osmextract OpenStreetMap の大きなデータセットをダウンロードし、インポート。
geodata 行政データ、標高データ、WorldClim データのダウンロードとインポート。
rnaturalearth Natural Earth ベクタ・ラスタデータ
rnoaa 米国海洋大気庁(National Oceanic and Atmospheric Administration, NOAA)の気候データをインポート。

Table 8.1 は、利用可能な地理データパッケージのごく一部に過ぎないことを強調しておく。 この他、tidycensustigris (USA)、cancensus (Canada)、eurostatgiscoR (European Union) あるいは idbr (international databases) など、様々な社会人口統計を取得する R パッケージが大量に存在している。Analyzing US Census Data (Walker 2022) には、こうしたデータを分析する方法がいくつか例示されている。 同様に、bcdata (Province of British Columbia)、geobr (Brazil)、RCzechia (Czechia)、rgugik (Poland) など、様々な地域や国の空間データにアクセスできる R パッケージが存在する。 その他の注目すべきパッケージは、R で Global Summary Daily Weather Data を提供する GSODR である(気象データソースの概要についてはパッケージの READMEを参照)。

各データパッケージは、データにアクセスするためのコードの書き方がそれぞれ異なる。 Table 8.1 の3つのパッケージについて、データを取得するコードチャンクを示すので、その違いを確認していただきたい。36 最初に、国の境界線は便利なことが多いので、rnaturalearth パッケージの ne_countries() 関数を用いて、以下のようにアクセスしてみよう。

library(rnaturalearth)
usa = ne_countries(country = "United States of America") # United States borders
class(usa)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
# alternative way of accessing the data, with geodata
# geodata::gadm("USA", level = 0, path = tempdir())

デフォルトでは、rnaturalearthSpatial* クラスのオブジェクトを返す。 この結果は、以下のように st_as_sf()sf オブジェクトに変換することができる。

usa_sf = st_as_sf(usa)

2つ目の例は、geodata パッケージを使用して、10分の空間分解能(赤道では約18.5 km)で全球の月別降水量の合計を含む一連のラスタをダウンロードしてみよう。 その結果、SpatRaster クラスのマルチレイヤオブジェクトが生成される。

library(geodata)
worldclim_prec = worldclim_global("prec", res = 10, path = tempdir())
class(worldclim_prec)

3つ目の例は、osmdata パッケージ (Padgham et al. 2018) を使って、OpenStreetMap (OSM) データベースから公園を検索してみる。 以下のコードチャンクに示すように、クエリーは関数 opq()(OpenStreetMap query の略)で始まり、その最初の引数は bounding box、つまり境界線(この場合はリーズの都市)を表すテキスト文字列である。 その結果は、どの OSM 要素(この場合は公園)に興味があるかを選択する関数に渡され、key-value ペアで表される。 次に、これらのデータは関数 osmdata_sf() に渡され、データのダウンロードと sf オブジェクトのリストへの変換が行われる(詳しくは vignette('osmdata') を参照)。

library(osmdata)
parks = opq(bbox = "leeds uk") |> 
  add_osm_feature(key = "leisure", value = "park") |> 
  osmdata_sf()

osmdata パッケージの制限は「容量制限」があり、大きな OSM データセット(例えば、大きな都市のすべての OSM データ)をダウンロードすることができない。 この制限を克服するために、osmextract パッケージが開発された。これは、あらかじめ定義された地域の OSM データベースの圧縮バージョンを含むバイナリ .pbf ファイルをダウンロードし、インポートすることができる。

OpenStreetMap は、クラウドソースによる膨大なグローバルデータベースであり、日々成長を続けている。また、OSM クエリの迅速な開発とテストを行うためのウェブサービス Overpass turbo から PostGIS データベースへのデータ取り込みを行うための osm2pgsql まで、データに容易にアクセスできるツールのエコシステムが充実している。 OSM から得られるデータセットの質は様々だが、データソースと OSM のエコシステムには多くの利点がある。データセットが世界中で利用でき、無料で、ボランティアの軍隊のおかげで常に改善されている。 OSM の利用は、「市民科学」とデジタルコモンズへの還元を促すものである(www.openstreetmap.org から、よく知る世界の一部を表すデータの編集を始めることができる)。 OSM データの活用例については、Chapter 10、Chapter 13、Chapter 14 を参照。

パッケージにデータセットが組み込まれていることがある。 この場合、アクセス方法は4つある。パッケージをアタッチする方法(パッケージが spData のように ‘lazy loading’ を使用している場合)は data(dataset, package = mypackage)、データセットを参照する方法は mypackage::dataset、生のデータファイルを参照する方法は system.file(filepath, package = mypackage) とする。 次のコードは、world データセット(親パッケージを library(spData) にアタッチしてロード済み)を使って、後者の2つのオプションを説明している。37

world2 = spData::world
world3 = read_sf(system.file("shapes/world.gpkg", package = "spData"))

最後の例、system.file("shapes/world.gpkg", package = "spData") は、spData パッケージの "shapes/" フォルダ内に格納されている world.gpkg ファイルへのパスを返す。

空間情報を得るもう一つの方法は、ジオコーディング(住所などの位置情報を座標に変換すること)である。 これは通常、オンラインサービスに問い合わせを行い、その結果として位置情報を取得するものである。 このようなサービスは数多く存在するが、使用するジオコーディングの方法、使用制限、コスト、API キーの要件などが異なっている。 R にはジオコーディングのためのパッケージがいくつかあるが、tidygeocoder は一貫したインタフェースで最も多くのジオコーディングサービスに接続することができるようである。 tidygeocoder のメイン関数は geocode で、アドレスを持つデータフレームを受け取り、"lat""long" として座標を追加する。 また、この関数は method の引数でジオコーディングサービスを選択することができ、多くの追加パラメータを持つ。

以下の例では、ロンドンのソーホー地区のビルにある John Snow(訳注:疫学的手法を導入しコレラの原因、感染経路を初めて特定した医師) の青い銘板の座標を検索している。

library(tidygeocoder)
geo_df = data.frame(address = "54 Frith St, London W1D 4SJ, UK")
geo_df = geocode(geo_df, address, method = "osm")
geo_df

得られたデータフレームは、st_as_sf() を用いて sf オブジェクトに変換することができる。

geo_sf = st_as_sf(geo_df, coords = c("lat", "long"), crs = "EPSG:4326")

また、このパッケージは、一組の座標に基づいて一連の情報(名前、住所など)を取得するために使用される逆ジオコーディングと呼ばれる逆の処理を実行することもできる。

8.4 地理ウェブサービス

空間データにアクセスするための Web API の標準化を目指して、Open Geospatial Consortium(OGC)は、Webサービス(OGC Web Services の略で OWS と総称)の仕様を多数策定している。 仕様には、Web Feature Service (WFS) 、Web Map Service (WMS) 、Web Map Tile Service (WMTS) 、Web Coverage Service (WCS) 、そして Web Processing Service (WPS) が含まれる。 PostGIS などの地図サーバーはこれらのプロトコルを採用し、クエリの標準化を進めている。 他のウェブ API と同様に、OWS API はデータを要求するために ? に続く ‘ベース URL’ と ‘エンドポイント’ および ‘URL クエリ引数’ を使う(httr パッケージ内の best-practices-api-packages vignette を参照)。

OWS のサービスには、さまざまな要望がある。 最も基本的なものの1つが getCapabilities であり、以下の httr 関数 GET()modify_url() で示されている。 このコードチャンクは、API のクエリを作成しデータ取得する方法を示している。この場合、国連食糧農業機関(Food and Agriculture Organization, FAO)が運営するサービスの能力を発見することもできる。

library(httr)
base_url = "http://www.fao.org"
endpoint = "/figis/geoserver/wfs"
q = list(request = "GetCapabilities")
res = GET(url = modify_url(base_url, path = endpoint), query = q)
res$url
#> [1] "https://www.fao.org/figis/geoserver/wfs?request=GetCapabilities"

上記のコードチャンクは、API リクエストを GET() 関数でプログラム的に構築する方法を示している。この関数は、ベース URL とクエリパラメータのリストを受け取り、簡単に拡張することができる。 リクエストの結果を、httr パッケージで定義されたクラス response のオブジェクト res に保存し、URL を含むリクエストの情報を含むリストとなる。 browseURL(res$url) を実行するとわかるように、結果はブラウザで直接読むこともできる。 リクエストの内容を抽出する一つの方法として、次のようなものがある。

txt = content(res, "text")
xml = xml2::read_xml(txt)
xml
#> {xml_document} ...
#> [1] <ows:ServiceIdentification>\n  <ows:Title>GeoServer WFS...
#> [2] <ows:ServiceProvider>\n  <ows:ProviderName>UN-FAO Fishe...
#> ...

WFS サービスからデータをダウンロードするには、GetFeature リクエストと特定の typeName (以下のコードチャンクに示す) が必要である。

利用できる名称は、アクセスする Web 機能サービスによって異なる。 GetCapabilities ウェブ技術を使ってプログラムで抽出することもでき、 (Nolan and Lang 2014)、ブラウザで出力された内容を手動でスクロールさせることもできる。

qf = list(request = "GetFeature", typeName = "area:FAO_AREAS")
file = tempfile(fileext = ".gml")
GET(url = base_url, path = endpoint, query = qf, write_disk(file))
fao_areas = read_sf(file)

write_disk() を使って、結果をメモリにロードされるのではなく、ディスクに書き込むようにすることで、sf でインポートできるようにすることに注意しておこう。 この例では、httr を使用して Web サービスに低レベルでアクセスする方法を示している。これは、Web サービスがどのように動作するかを理解するのに役立つ。 しかし、多くの日常的な作業には、より高度なインターフェースの方が適している場合があり、多くのRパッケージやチュートリアルは、まさにこの目的のために開発されたものである。 OWS サービスを利用する ows4R というパッケージが開発された。

8.5 ファイル形式

地理データセットは通常、ファイルまたは空間データベースとして保存される。 ファイル形式はベクタデータとラスタデータのどちらかを保存できるが、 PostGIS のような空間データベースは両方を保存できる ( Section 10.6.2 も参照)。 今日、ファイル形式の多様性は困惑するほどある。しかし、1960年代には、ハーバード大学で空間解析のための最初の広く配布されたプログラム(SYMAP)などの初期の GIS ソフトウェアが開発された。それ以降、多くの統合と標準化が行われてきた (Coppock and Rhind 1991)

GDAL(Geospatial Data Abstraction Library、「グードル」と発音する。「グー」の goo の oo 部分は、Object Oeirneted を表す)は、2000年のリリース以来、地理ファイルフォーマット間の非互換性に関連する多くの問題を解決した。 GDAL は、多くのラスタおよびベクタデータフォーマットの読み書きのための統一された高性能なインタフェースを提供する。38 GRASS、ArcGIS、QGIS など、多くのオープンおよびプロプライエタリな GIS プログラムは、GUI の背後に GDAL を使用して、地理データを取り込み、適切な形式で出力するという足回りの作業を行なっている。

GDAL は、200 以上のベクタおよびラスタデータフォーマットへのアクセスを提供する。 Table 8.2 では、よく使われる空間ファイルフォーマットについての基本情報を紹介している。

TABLE 8.2: 代表的な空間ファイル形式。
名称 拡張子 情報 タイプ モデル
ESRI Shapefile .shp (メインとなるファイル) よく使われているフォーマットで、少なくとも3つのファイルから構成される。ファイルサイズが 2GB を超えるもの、種類が混在するもの、名前が10文字以上のもの、列数が255以上のものはサポートされていない。 ベクタ Partially open
GeoJSON .geojson JSON 交換フォーマットを、シンプルフィーチャを含むように 拡張したもので、主に経度・緯度の座標を格納するために使用され、TopoJSON フォーマットによって拡張される。 ベクタ Open
KML .kml Google Earth で使用するために開発された、XML ベースの空間可視化フォーマット。ZIP 形式の KML ファイルは KMZ 形式。 ベクタ Open
GPX .gpx GPS データ交換のために作成された XML スキーマ。 Vector Open
FlatGeobuf .fgb ベクターデータを高速に読み書きできる単一ファイル形式。ストリーミング機能を持つ。 ベクタ Open
GeoTIFF .tif/.tiff 一般的なラスタフォーマット。空間メタデータを追加したTIFFファイル。 ラスタ Open
Arc ASCII .asc 最初の6行がラスタヘッダーで、その後にラスタセルの値が行と列に並んでいるテキスト形式。 ラスタ Open
SQLite/SpatiaLite .sqlite スタンドアローンのリレーショナルデータベースである SpatiaLite は、SQLite の空間拡張版。 ベクタとラスタ Open
ESRI FileGDB .gdb ArcGIS で作成された空間および非空間オブジェクト。可能なこと:複数のフィーチャクラス、トポロジー。GDAL によるサポートは限定される。 ベクタとラスタ Proprietary
GeoPackage .gpkg SQLite をベースとした軽量なデータベースコンテナで、プラットフォームに依存しないジオデータの交換を容易に行うことができます。 ベクタと(制限のある)ラスタ Open

ファイルフォーマットの標準化とオープンソースを保証する重要な発展は、1994年の Open Geospatial Consortium (OGC) の設立であった。 OGC は、シンプルフィーチャというデータモデル(Section 2.2.1 参照)を定義するだけでなく、例えば KML や GeoPackage などのファイルフォーマットで使用されているようなオープンスタンダードの開発も調整している。 OGC が推奨するオープンなファイルフォーマットは、プロプライエタリなフォーマットと比較して、いくつかの利点がある。標準が公開され、透明性が確保され、ユーザーがファイルフォーマットをさらに開発し、特定のニーズに合わせて調整する可能性が開かれることである。

ESRI Shapefile は最も一般的なベクタデータ交換フォーマットであるが、オープンフォーマットではない(仕様はオープン)。 1990年代初頭に開発されたもので、多くの制約がある。 まず、少なくとも3つのファイルで構成されるマルチファイル形式であること。 255列までしかサポートしておらず、列名は10文字まで、ファイルサイズは 2 GB までと制限されている。 さらに、ESRI Shapefile は、ポリゴンと複合ポリゴンの区別ができないなど、可能なすべてのジオメトリタイプをサポートしていない。39 このような制約があるにもかかわらず、長い間、有力な代替手段が見つかっていなかった。 最近では、GeoPackage が登場し、ESRI Shapefile に取って代わろうとしている。 Geopackage は、地理空間情報を交換するためのフォーマットで、OGC 規格の一つである。 GeoPackage 規格は、地理空間情報を SQLite の小さなコンテナに格納する方法についての規則を記述している。 したがって、GeoPackage は軽量な空間データベースコンテナであり、ベクタおよびラスタデータだけでなく、非空間データおよび拡張機能も格納することができる。 GeoPackage 以外にも、調べる価値のある地理空間データ交換フォーマットがある(Table 8.2)。

ラスタデータの形式としては、GeoTIFF 形式が主流である。 TIFF ファイル内に CRS などの空間情報を埋め込むことができる。 ESRI Shapefile と同様に1990年代に開発されたフォーマットで、オープンなフォーマットである。 さらに、GeoTIFF は現在も拡張・改良が続けられている。 GeoTIFF フォーマットに最近追加された最も重要なものの1つが、COG(Cloud Optimized GeoTIFF)と呼ばれるバージョンである。 COG として保存されたラスタオブジェクトは、HTTP サーバーでホストすることで、他の人がファイル全体をダウンロードすることなく、ファイルの一部だけを読むことができる(Section 8.6.2 と Section 8.7.2 を参照)。

その他にも、書籍の制限上、詳細な説明や Table 8.2 に掲載していない空間データフォーマットが多数存在する。 他のフォーマットを使用する必要がある場合は、ベクタおよびラスタドライバに関する GDAL のドキュメントを読むことを勧める。 さらに、空間データフォーマットの中には、ベクタやラスタ以外のデータモデル(タイプ)を格納できるものもある。 ライダー点群を格納するための LAS、LAZ 形式、多次元配列を格納するための NetCDF、HDF 形式が含まれる。

また、空間データは、CSV ファイルや Excel スプレッドシートなど、表形式(非空間)のテキスト形式で保存されることも多い。 例えば、GIS ツールを使わない人と空間サンプルを共有したり、空間データ形式を受け付けない他のソフトウェアとデータを交換したりする際に便利である。 しかし、この方法は、点よりも複雑な形状を保存するにはかなり困難であり、CRS に関する情報を直接保存できないなど、いくつかの問題が考えられる。

8.6 データ入力 (I)

sf::read_sf() (ベクタデータの読み込みに使うメイン関数) や terra::rast() (ラスタデータの読み込みに使うメイン関数) などのコマンドを実行すると、ファイルからデータを読み込むイベントの連鎖が無言で開始される。 さらに、さまざまな地理データを含む、あるいは異なるデータソースに簡単にアクセスできる R パッケージが数多く存在する。 これらはすべて、R にデータをロードするか、より正確には、RAM に保存されたワークスペースにオブジェクトを割り当てて、そこからアクセスできるようにするものである。R セッション中の .GlobalEnv からアクセスできる RAM に保存されている。

8.6.1 ベクタデータ

空間ベクタデータは、さまざまなファイル形式で提供されている。 .geojson.gpkg ファイルなど、よく使われる表現のほとんどは、裏で GDAL のベクタドライバ を使う sf 関数 read_sf()(または同等の st_read())で直接 R に取り込むことができる。 st_drivers() は、最初の2列に namelong_name を含むデータフレームを返す。続く列に、Table 8.3 の主要ファイルフォーマットについて図示しているように、データの書き込みやラスタデータの保存など GDAL(従って sf)で利用できる各ドライバの機能を返す。
以下のコマンドは、コンピュータにインストールされた GDAL の最初の3つのドライバを報告し(結果はインストールされた GDAL のバージョンによって異なる場合がある)、それらのフィーチャの要約を表示す。 なお、大半のドライバはデータの書き込みが可能であるが(87種類中51種類)、ベクタデータに加えてラスタデータを効率的に表現できるフォーマットは16種類しかない(詳しくは ?st_drivers())。

sf_drivers = st_drivers()
head(sf_drivers, n = 3)
summary(sf_drivers[-c(1:2)])
TABLE 8.3: ベクターデータを読み書きするための一般的な ドライバ/フォーマット。
name long_name write copy is_raster is_vector vsi
ESRI Shapefile ESRI Shapefile TRUE FALSE FALSE TRUE TRUE
GPX GPX TRUE FALSE FALSE TRUE TRUE
KML Keyhole Markup Language (KML) TRUE FALSE FALSE TRUE TRUE
GeoJSON GeoJSON TRUE FALSE FALSE TRUE TRUE
GPKG GeoPackage TRUE TRUE TRUE TRUE TRUE
FlatGeobuf FlatGeobuf TRUE FALSE FALSE TRUE TRUE

read_sf() の第一引数は dsn で、これはテキスト文字列または単一のテキスト文字列を含むオブジェクトであるべきである。 テキスト文字列の内容は、ドライバによって異なる可能性がある。 多くの場合、ESRI Shapefile(.shp)や GeoPackage 形式(.gpkg)と同様に、dsn はファイル名となる。 read_sf() は、ファイルの拡張子からドライバを推測する(下の例は、.gpkg の場合)。(訳注:日本語データを含むファイルは、文字エンコーディングを指定しないと文字化けを起こす。ほとんどの場合、CP932 を使用しているので、read_sf() に引数 options = "ENCODING=CP932" を設定するとよい。)

f = system.file("shapes/world.gpkg", package = "spData")
world = read_sf(f, quiet = TRUE)

ドライバによっては、dsn は、フォルダ名、データベースのアクセス認証情報、または GeoJSON 文字列表現として提供されることがある(詳細は、read_sf() のヘルプページの例を参照)。

ベクタドライバのフォーマットには、複数のデータレイヤを格納できるものがある。 デフォルトでは、read_sf()dsn で指定されたファイルの最初のレイヤーを自動的に読み込む。しかし、layer 引数を使用すると、他のレイヤーを指定することができる。

また、read_sf() 関数は、ファイルの一部だけを RAM に読み込むことができる方法が2つある。 1つ目は、query の引数に関連して、OGR SQL query text でデータのどの部分を読み取るかを指定できるようにしたものである。 以下の例では、タンザニアのみのデータを抽出している(Figure 8.1 :A)。 これは、"world" のレイヤから、name_long"Tanzania" と等しいすべての列 (SELECT *) を取得したい、と指定することで実現する。

tanzania = read_sf(f, query = 'SELECT * FROM world WHERE name_long = "Tanzania"')

利用可能な列の名前がわからない場合、'SELECT * FROM world WHERE FID = 1' でデータの1行だけを読み込むのが良い方法である。 FIDフィーチャ IDを表し、多くの場合行番号であるが、その値は使用するファイル形式に依存する。 例えば、FID は、ESRI Shapefile では 0 から始まり、他のファイルフォーマットでは 1 から始まるか任意とすることができる。

2つ目の仕組みは、wkt_filter の引数を使用する。 この引数は、データを抽出したい研究領域を表すよく知られたテキストを想定している。 タンザニアの国境線 50,000 m と交差するポリゴンをファイルから読み込む。 そのためには、(a) バッファを作成するか(Section 5.2.3)、(b) st_geometry()sf バッファオブジェクトを sfc ジオメトリオブジェクトに変換するか、(c) st_as_text() でジオメトリを WKT に変換する、のいずれかの方法で「フィルタ」を準備する必要がある。

tanzania_buf = st_buffer(tanzania, 50000)
tanzania_buf_geom = st_geometry(tanzania_buf)
tanzania_buf_wkt = st_as_text(tanzania_buf_geom)

さて、この「フィルタ」を wkt_filter の引数で適用してみよう。

tanzania_neigh = read_sf(f, wkt_filter = tanzania_buf_wkt)

Figure 8.1 :B に示すように、この結果はタンザニアとその 50 km バッファ内のすべての国を含めている。

クエリ(A)と wkt フィルターを用いて、ベクタデータのサブセットを読み込む(B)。

FIGURE 8.1: クエリ(A)と wkt フィルターを用いて、ベクタデータのサブセットを読み込む(B)。

当然ながら、一部のオプションは特定のドライバに固有のものである。40 例えば、表計算ソフトのフォーマット(.csv)に保存された座標を考えてみよう。 このようなファイルを空間オブジェクトとして読み込むには、当然、座標を表す列の名前(以下の例では、XY)を指定しなければならない。 これは、options パラメータの助けを借りて行うことができる。 可能なオプションについては、対応する GDAL ドライバの説明の「オープンオプション」セクションを参照。 カンマ区切り値(csv)形式は、http://www.gdal.org/drv_csv.html

cycle_hire_txt = system.file("misc/cycle_hire_xy.csv", package = "spData")
cycle_hire_xy = read_sf(cycle_hire_txt,
  options = c("X_POSSIBLE_NAMES=X", "Y_POSSIBLE_NAMES=Y"))

「XY」座標を記述する代わりに、1つの列でジオメトリ情報を記述することも可能である。 Well-known text (WKT)、well-known binary (WKB)、GeoJSON 形式がその例である。 例えば、world_wkt.csv のファイルには、世界の国々のポリゴンを表す WKT という列がある。 このことを示すために、今回も options パラメータを使用する。

world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt2 = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT", 
                    quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE)
サポートされているすべてのベクタファイル形式が、その座標参照系に関する情報を格納しているわけではない。 このような場合、st_set_crs() 関数を用いて不足する情報を追加することが可能である。 詳細は Section 2.4 も参照。

最後の例として、read_sf() が KML ファイルも読み込むことを紹介する。 KML ファイルは、地理情報を XML 形式で格納している。これは、アプリケーションに依存しない方法で Web ページを作成し、データを転送するためのデータ形式である (Nolan and Lang 2014)。 ここでは、ウェブから KML ファイルにアクセスする。 このファイルには、複数のレイヤが含まれている。 st_layers() は、利用可能なすべてのレイヤを表示する。 read_sf()layer パラメータの助けを借りて最初のレイヤ Placemarks を選択する。

u = "https://developers.google.com/kml/documentation/KML_Samples.kml"
download.file(u, "KML_Samples.kml")
st_layers("KML_Samples.kml")
#> Driver: KML 
#> Available layers:
#>              layer_name  geometry_type features fields crs_name
#> 1            Placemarks       3D Point        3      2   WGS 84
#> 2      Highlighted Icon       3D Point        1      2   WGS 84
#> 3                 Paths 3D Line String        6      2   WGS 84
#> 4         Google Campus     3D Polygon        4      2   WGS 84
#> 5      Extruded Polygon     3D Polygon        1      2   WGS 84
#> 6 Absolute and Relative     3D Polygon        4      2   WGS 84
kml = read_sf("KML_Samples.kml", layer = "Placemarks")

このセクションで紹介した例はすべて、地理データのインポートに sf パッケージを使用したものである。 高速で柔軟性があるが、特定のファイル形式については他のパッケージも見てみる価値があるだろう。 例として、geojsonsf パッケージがある。 ベンチマーク によると、.geojson を読むのに sf パッケージの10倍程度の速度がある。

8.6.2 ラスタデータ

ラスタデータは、ベクタデータと同様に多くのファイル形式があり、中にはマルチレイヤファイルをサポートするものもある。 terrarast() コマンドは、レイヤが1つだけのファイルが提供された場合、1つのレイヤで読み込む。

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

また、マルチレイヤファイルを読み込む場合にも有効である。

multilayer_filepath = system.file("raster/landsat.tif", package = "spDataLarge")
multilayer_rast = rast(multilayer_filepath)

これまでの例はすべて、ハードディスクに保存されているファイルから空間情報を読み取るものであった。 しかし、GDAL は HTTP/HTTPS/FTP の Web リソースなど、オンラインのリソースから直接データを読み込むことも可能である。 あとは、ファイルへのパスの前に /vsicurl/ というプレフィックスを付けるだけである。 試しに、2000年から2012年までの 500 m 解像度での全球の月別積雪確率に接続してみよう (T. Hengl 2021)。 12月の積雪確率は、COG(Cloud Optimized GeoTIFF)ファイル(Section 8.5)として、 に保存されている。 オンラインファイルを読むには、そのURLと /vsicurl/ プレフィックスを指定するだけである。

myurl = "/vsicurl/https://zenodo.org/record/5774954/files/clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif"
snow = rast(myurl)
snow
#> class       : SpatRaster 
#> dimensions  : 35849, 86400, 1  (nrow, ncol, nlyr)
#> resolution  : 0.00417, 0.00417  (x, y)
#> extent      : -180, 180, -62, 87.4  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif 
#> name        : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0

入力データが COG であるため、実際にはこのファイルを RAM に読み込むのではなく、値を取得せずに接続を作成している。 その値は、何らかの値に基づく操作(例えば、crop()extract())を適用した場合に読み取られる。 これにより、ファイル全体をダウンロードすることなく、データのごく一部だけを読み出すこともできるようになった。 例えば、レイキャビクの座標を指定し、extract() 関数を適用すると、12 月の積雪確率(70%)を求めることができる。

rey = data.frame(lon = -21.94, lat = 64.15)
snow_rey = extract(snow, rey)
snow_rey
#>   ID clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0
#> 1  1                                                         70

この方法では、大きな GeoTIFF ファイル全体をダウンロードするのではなく、一つの値だけをダウンロードすることになる。

上記の例は、単純な(しかし有用な)1つのケースを示しただけであるが、もっと探求すべきことがある。 また、/vsicurl/ のプレフィックスは、ラスタだけでなく、ベクタファイルフォーマットにも有効である。 ベクタファイルのURLの前に接頭辞を付けるだけで、read_sf()、オンラインストレージから直接ベクタを読み込むことができるようになる。

重要なのは、GDAL が提供する接頭辞は /vsicurl/ だけではないことである。ZIP アーカイブから空間ファイルを解凍せずに読み込むための /vsizip/ や、AWS S3 バケットにあるファイルをオンザフライで読み込むための /vsis3/ など、他にも多くの接頭辞が存在するのである。 詳しくは、https://gdal.org/user/virtual_file_systems.html

8.7 データ出力(O)

地理データの書き込みでは、あるフォーマットから別のフォーマットへの変換や、新しく作成したオブジェクトの保存が可能である。 データの種類(ベクタまたはラスタ)、オブジェクトのクラス(例: sf または SpatRaster )、保存される情報の種類と量(オブジェクトのサイズ、値の範囲など)に応じて、空間ファイルを最も効率的に保存する方法を知ることが重要である。 次の2つのセクションでは、その方法を説明する。

8.7.1 ベクタデータ

read_sf() の対語は write_sf().geojson.shp.gpkg などの最も一般的なものを含む、広範囲の地理ベクタファイルフォーマットに sf オブジェクトを書き込むことができる。 ファイル名から、write_sf() が自動的に使用するドライバを決定する。 また、書き込み速度はドライバに依存する。

write_sf(obj = world, dsn = "world.gpkg")

注意:同じデータソースに再度書き込もうとすると、この機能はファイルを上書きしてしまう。

write_sf(obj = world, dsn = "world.gpkg")

ファイルを上書きする代わりに、append = TRUE、ファイルに新しいレイヤーを追加することができる。これは、GeoPackageを含むいくつかの空間フォーマットでサポートされている。

write_sf(obj = world, dsn = "world_many_layers.gpkg", append = TRUE)

また、write_sf() と同等である、st_write() を使用することもできる。 ただし、デフォルトは異なる。ファイルを上書きしない(上書きしようとするとエラーを返す)、書き込まれたファイル形式とオブジェクトの短い要約を表示する、などである。

st_write(obj = world, dsn = "world2.gpkg")
#> Writing layer `world2' to data source `world2.gpkg' using driver `GPKG'
#> Writing 177 features with 10 fields and geometry type Multi Polygon.

layer_options 引数はまた、さまざまな目的で使用することができる。 そのひとつが、空間データをテキストファイルに書き出すことである。 これは、layer_options の中に GEOMETRY を指定することで可能である。 単純な点データセットの場合は AS_XY (座標のための新しい列を2つ作成する)、より複雑な空間データの場合は AS_WKT (空間オブジェクトのよく知られたテキスト表現を含む新しい列を1つ作成する) のいずれかになる。

write_sf(cycle_hire_xy, "cycle_hire_xy.csv", layer_options = "GEOMETRY=AS_XY")
write_sf(world_wkt, "world_wkt.csv", layer_options = "GEOMETRY=AS_WKT")

8.7.2 ラスタデータ

writeRaster() 機能は、SpatRaster のオブジェクトをディスク上のファイルに保存する。 この関数は、出力データ型とファイルフォーマットに関する入力を期待するが、選択されたファイルフォーマットに固有のGDALオプションも受け付ける(詳しくは ?writeRaster を参照)。

ラスタを保存する際、terra パッケージは以下の7つのデータ形式を提供する:INT1U、INT2S、INT2U、INT4S、INT4U、FLT4S、FLT8S。41 データ型は、ディスクに書き込まれるラスタオブジェクトのビット表現を決定する(Table 8.4)。 どのデータ型を使用するかは、ラスタオブジェクトの値の範囲による。 データ型が表現できる値が多いほど、ディスク上のファイルサイズは大きくなる。 符号なし整数(INT1U、INT2U、INT4U)はカテゴリデータに適しており、浮動小数点数(FLT4S、FLT8S)は通常連続データを表す。 writeRaster() は FLT4S をデフォルトとして使用する。 これはほとんどの場合において有効であるが、二値やカテゴリーデータを保存する場合、出力ファイルのサイズは不必要に大きくなる。 したがって、最小限の記憶容量を必要とし、なおかつすべての値を表現できるデータ型を使用することを勧める(summary() 関数で値の範囲を確認する)。

TABLE 8.4: terra パッケージが対応しているデータ型。
データタイプ 最小値 最大値
INT1U 0 255
INT2S -32,767 32,767
INT2U 0 65,534
INT4S -2,147,483,647 2,147,483,647
INT4U 0 4,294,967,296
FLT4S -3.4e+38 3.4e+38
FLT8S -1.7e+308 1.7e+308

デフォルトでは、出力ファイル形式はファイル名から導かれる。 ファイル名を *.tif とすると、以下のように GeoTIFF ファイルが作成される。

writeRaster(single_layer, filename = "my_raster.tif", datatype = "INT2U")

ラスタファイルフォーマットによっては、追加オプションがあり、writeRaster()options 引数に GDAL parameters を与えることで設定できる。 GeoTIFF ファイルは、デフォルトで terra で記述され、LZW 圧縮が施されている gdal = c("COMPRESS=LZW")。 圧縮を変更したり無効にしたりするには、この引数を変更する必要がある。

writeRaster(x = single_layer, filename = "my_raster.tif",
            gdal = c("COMPRESS=NONE"), overwrite = TRUE)

さらに、filetype = "COG" のオプションでラスタオブジェクトを COG (Cloud Optimized GeoTIFF, Section 8.5 ) として保存することができる。

writeRaster(x = single_layer, filename = "my_raster.tif",
            filetype = "COG", overwrite = TRUE)

8.8 ビジュアル出力

R は、多くの静的および対話的なグラフィックス形式をサポートしている。 静的プロットを保存する最も一般的な方法は、例えばグラフィックデバイスを開き、プロットを作成し、それを閉じることである。

png(filename = "lifeExp.png", width = 500, height = 350)
plot(world["lifeExp"])
dev.off()

この他の利用可能なグラフィックデバイスには、pdf()bmp()jpeg()tiff() がある。 出力プロットの幅、高さ、解像度など、いくつかのプロパティを指定することができる。

さらに、いくつかのグラフィックパッケージは、グラフィック出力を保存するための独自の関数を提供している。 例えば、tmap パッケージには、tmap_save() という関数がある。 オブジェクト名と新規ファイルへのファイルパスを指定することで、tmap オブジェクトをさまざまなグラフィックフォーマットまたは HTML ファイルに保存することができる。

library(tmap)
tmap_obj = tm_shape(world) + tm_polygons(col = "lifeExp")
tmap_save(tmap_obj, filename = "lifeExp_tmap.png")

一方、mapview パッケージで作成したインタラクティブマップは、mapshot() 関数を使用して HTML ファイルまたは画像として保存することができる。

library(mapview)
mapview_obj = mapview(world, zcol = "lifeExp", legend = TRUE)
mapshot(mapview_obj, file = "my_interactive_map.html")

8.9 演習

E1. List and describe three types of vector, raster, and geodatabase formats.

E2. Name at least two differences between the sf functions read_sf() and st_read().

E3. Read the cycle_hire_xy.csv file from the spData package as a spatial object (Hint: it is located in the misc folder). What is a geometry type of the loaded object?

E4. Download the borders of Germany using rnaturalearth, and create a new object called germany_borders. Write this new object to a file of the GeoPackage format.

E5. Download the global monthly minimum temperature with a spatial resolution of five minutes using the geodata package. Extract the June values, and save them to a file named tmin_june.tif file (hint: use terra::subset()).

E6. Create a static map of Germany’s borders, and save it to a PNG file.

E7. Create an interactive map using data from the cycle_hire_xy.csv file. Export this map to a file called cycle_hire.html.