💬

[R] openeo パッケージで Sentinel-2 のデータをダウンロードする

に公開

tl;dr

ダウンロードするだけなら OpenEO じゃなくて CDSE の方が楽そう

きっかけ

Pythonで学ぶ衛星データ解析基礎」を 勉強のために R でやってみよう、と思ったところ、さっそく壁にぶち当たります。本で紹介されている Sentinel Hub はもう閉鎖されていて、Sentilen-2 のデータをダウンロードするには別のサービスを使う必要があります。そのひとつが Copernicus Data Space Ecosystem(CDSE)で、R から CDSE の API を使う方法は以下の記事に解説されていました。

https://qiita.com/ade_san/items/be822f97d3ec270d9b98

ここで、この記事中でも紹介されていますが、CDSE 以外に、OpenEO というサービスもあり、こっちは公式の R クライアントが用意されているので、こっちの方がいいのでは?(CDSE API 用の R パッケージは個人が作っているもの)、という説があります。ということで軽く使い方を調べてみました。

https://open-eo.github.io/openeo-r-client/

OpenEO とは

OpenEO が何なのかというのは、全容を理解できていないのですが、About には「衛星データはローカルで処理するにはでかくなりすぎた。また、データへのアクセス方法や処理方法も多様化して複雑になりすぎた。そこで大統一クラウド GIS API ですよ!(超意訳)」みたいなことが書かれています。GIS サービスによっていろいろ仕様が違うのを、この OpenEO API を通せば統一的に使える、という世界観を目指しているみたいです。

けっこう前のブログですが、こっちの方が壮大さが伝わるかと思います。

https://r-spatial.org/2016/11/29/openeo.html

と聞いて、どう思ったでしょうか。あの...、その話、長くなります?? おれはデータをダウンロードしたいだけなんだが...、という気持ちになりましたか? そう、なんか薄々感じるように、めっちゃ高機能で便利なんですけど、「とりあえず衛星画像を数枚ダウンロードしたい」みたいな時にはちょっとオーバースペックな代物なんです。

まあとはいえ、かなり資金を集めて作られているだけあって、少なくとも R クライアントはそこそこちゃんとした品質のものになっています。覚えておいて損はないかもしれません。

準備

アカウント

Getting Started によれば EGI にアカウントを作ることになっていますが、私は作ってなくても使えました。よくわからないのですが、ここは実際に OpenEO のバックエンドとして使うサービスのアカウントがあればいいということなんだと思います。私は CDSE を使ったので、CDSE のアカウントがあれば大丈夫でした。

R パッケージのインストール

CRAN からインストールできます。

install.packages("openeo")

認証

途中まで認証なしで使えるんですが、どうせあとでログインすることになるので先にしておきましょう。ちなみに、connect() は重要です。明示的に関数に渡されませんが、active_connection() の結果が勝手に参照されることになります。なので、たぶん connect() を先に実行しておかないといろいろなものが動きません。

library(openeo)

# OpenEO のバックエンドとしてどのサービスを使うかを選ぶ。これは CSDE の URL。
con <- connect(host = "https://openeofed.dataspace.copernicus.eu")

login()

実行するとブラウザに飛ばされ、私の場合は以下のような認証画面が出ていました。許可すると認証が必要な API が使えるようになります。

使い方

さて、基本的な使い方がどんな感じかみていきましょう。

まず、定型文的にこれら4つのオブジェクトを作成します。

con <- connect(host = "https://...")

collections <- list_collections()
formats <- list_file_formats()

p <- processes()

connect() は上で説明したので飛ばします。

list_collections() は、利用可能なコレクション、つまりデータの種類をリストアップします。このオブジェクトは、どのコレクションを使うかを指定するために使われます。たとえば、describe_collection() は指定したコレクションのメタデータを表示してくれます。

describe_collection(colls$SENTINEL2_L2A)
#> SENTINEL2_L2A
#> Title:				Sentinel-2 L2A
#> Description:			SENTINEL-2 is a wide-swath, high-resolution, multi-spectral imaging mission, supporting Copernicus Land Monitoring studies...
#>  The level 2A data is atmospherically corrected using Sen2Cor.
#> Source:				ESA, CDSE
#> Platform:			sentinel-2a
#> Platform:			sentinel-2b
#> Constellation:			sentinel-2
#> Instrument:			---
#> Spatial extent (lon,lat):	(-180, -56), (180, 83)
#> Temporal extent:		2015-07-04T00:00:00Z / NA
#> Bands:
#>                                                       aliases center_wavelength common_name full_width_half_max gsd             name offset scale
#> 1                            IMG_DATA_Band_B01_60m_Tile1_Data            0.4427     coastal               0.021  60              B01      0 1e-04
#> ...

また、list_file_formats() は入出力でサポートされているデータフォーマットをリストアップします。これも同様に、入出力の形式を指定するために使います。

processes() は、データ処理のパイプラインを構築します。ここに様々な処理を指定しくわけですが、処理はその瞬間に実行されるわけではなく、最後に compute_result() のような特定の関数を呼び出した時にまとめて実行されるみたいです。dbplyr の lazy data frame みたいなものだと思うとよさそうです。

load_collection() は、指定した条件でデータを検索して読み込みます。いや、読み込みます、というのは正確ではなく、読み込むという処理を定義します。これを実行すると一瞬で返ってきて、「あれ、もっとデータ重いのでは?」と感じるかもしれませんが、この時点ではまだ実際のデータ読み込みは実行されていません。

data <- p$load_collection(
  id = colls$SENTINEL2_L2A,
  bands = c("B02", "B03", "B04"),
  spatial_extent = list(
    west = 139.7081,
    south = 35.4818,
    east = 139.8620,
    north = 35.6594
  ),
  temporal_extent = c("2021-11-10", "2021-11-20")
)

そして、本来は、これに続いて様々な計算を定義します。たとえば Geocomputation with R では、NDVI の計算が定義されています。このように、生データをダウンロードするのではなく、データの処理も含めてクラウド上で行って、最終的な成果物だけを落としてくるように設計されています。

# NDVI vegetation index を計算
compute_ndvi <- p$reduce_dimension(data, dimension = "bands",
                                   reducer = function(data, context) {
                                       (data[2] - data[1]) / (data[2] + data[1])
                                   })
# maximum over time を計算
reduce_max <- p$reduce_dimension(data = compute_ndvi, dimension = "t",
                                 reducer = function(x, y) {
                                     max(x)
                                 })

しかし、私がやりたいのはただデータをダウンロードすることなので、さっさと結果をまとめましょう。まず、save_result() はファイルに書き出すという処理を定義します。そして、compute_result() では、構築したデータ処理パイプラインを実行し、その結果を手元に持ってきます。

NetCDF ファイルとしてダウンロードしてくるのは以下のようなコードになります。この compute_result() の実行にはけっこう(数分程度?)時間がかかります。また、もし HTTP 401 Unauthorized というエラーが出た場合は login() でログインしなおしましょう。

result <- p$save_result(data = data, format = formats$output$netCDF)
compute_result(result, output_file = "tmp.nc")

stars 形式で直接 R に読み込みたい場合は以下のように指定します。

stars_df <- compute_result(result, as_stars = TRUE)

とりあえずプロットしてみるとそれっぽい図が描かれるので、あってそうです。

plot(stars_df)


(出典:Copernicus Data Space Ecosystem 提供の Sentinel-2 の衛星画像データを使って描画)

感想

ということで、OpenEO は高機能なのですが、ちょっと手元でデータと戯れたいだけみたいな時には格式高すぎるかもなあ、というのが正直な感想です。CDSE パッケージだと CDSE 依存になっちゃうわけですが、とりあえず今はそっちを使おうかなと思います。

Discussion