🙌

『Python で学ぶ衛星データ解析基礎』を R でやってみる (1) 3章前半まで

2025/05/11に公開

本

この本がずっと積読になってたのでやってみようと思い。

3-1-4～3-1-5 APIでSentinel-2の画像を取得する

本で使われている Sentinel Hub はもう閉鎖されており、Sentilen-2 のデータをダウンロードするには別のサービスを使う必要があります。このうち、OpenEO を使う方法は以下の記事に書きました。

今回は、OpenEO 経由ではなく Copernicus Data Space Ecosystem（CDSE）を直接使ってみます。CDSE の API を使うには CDSE パッケージが便利です。使い方は以下で解説されていたので、詳しくはこっちを読んでください！、ということにしてここでは説明をサボろうと思います。

準備

公式ドキュメントにあるように、まずは CDSE にアカウントを作成します。

そして、OAuth アプリを作成し、client ID と client secret を適当なところに保存します。公式ドキュメントでは環境変数に保存するよう言っていますが、ここでは keyring パッケージを使ってみます。

keyring::key_set("CDSE-client-id")
keyring::key_set("CDSE-client-secret")

コレクションの確認

まずは、GetCollections() で API が提供している衛星データ一覧を確認しましょう。id が sentinel-2-l2a のものが使いたいコレクションです。

library(CDSE)

collections <- GetCollections(as_data_frame = TRUE)

collections
#>                      id                 title
#> 1        sentinel-2-l1c        Sentinel 2 L1C
#> 2    sentinel-3-olci-l2    Sentinel 3 OLCI L2
#> 3       sentinel-3-olci       Sentinel 3 OLCI
#> 4      sentinel-3-slstr      Sentinel 3 SLSTR
#> 5 sentinel-3-synergy-l2 Sentinel 3 Synergy L2
#> 6        sentinel-1-grd        Sentinel 1 GRD
#> 7        sentinel-2-l2a        Sentinel 2 L2A
#> 8        sentinel-5p-l2  Sentinel 5 Precursor
#> ...

衛星データの検索

衛星データの検索は SearchCatalog() を使います。この API は認証が必要なので、 GetOAuthClient() でつくった OAuth client を client オプションに渡します。この結果は sf オブジェクトになっているので、必要に応じて dplyr::arange() で並べ替えたり dplyr::filter() でさらに絞り込んだりするとよいでしょう。

本との違いを挙げると、

本では12 シーンがヒットしたようですが、ここでは2シーンだけでした。
被雲率で検索結果を絞り込んでいますが、CDSE パッケージにはそのようなオプションはないようです。filter 引数に CQL2 で絞り込み条件を書けるので、おそらく同等のことはできるんだと思います（まだ CQL2 をよく知らなくてわからない）。

ちなみに、bbox オプションはこの記事を書いている時点ではバグっていて動かないので、GitHub 版をインストールする必要があるかもしれません。bbox の代わりに aoi に対象範囲のポリゴンの sf オブジェクトを指定するのでも大丈夫です。

OAuthClient <- GetOAuthClient(
  id = keyring::key_get("CDSE-client-id"),
  secret = keyring::key_get("CDSE-client-secret")
)

catalog <- SearchCatalog(
  bbox = c(139.7081, 35.4818, 139.8620, 35.6594),
  from = "2021-11-10",
  to = "2021-11-20",
  collection = "sentinel-2-l2a",
  with_geometry = TRUE,
  client = OAuthClient
)

# 被雲率で並べ替え
catalog |> 
  dplyr::arrange(tileCloudCover)
#> Simple feature collection with 2 features and 11 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 138.7776 ymin: 35.13501 xmax: 140.0097 ymax: 36.14069
#> Geodetic CRS:  WGS 84
#>   acquisitionDate tileCloudCover areaCoverage   satellite
#> 1      2021-11-18          99.94          100 sentinel-2b
#> 2      2021-11-13           0.05          100 sentinel-2a
#>   acquisitionTimestampUTC acquisitionTimestampLocal
#> 1     2021-11-18 01:37:20       2021-11-18 10:37:20
#> 2     2021-11-13 01:37:24       2021-11-13 10:37:24
#>                                                            sourceId long.min
#> 1 S2B_MSIL2A_20211118T012939_N0500_R074_T54SUE_20230101T084609.SAFE 138.7776
#> 2 S2A_MSIL2A_20211113T012911_N0500_R074_T54SUE_20221231T081817.SAFE 138.7776
#>    lat.min long.max  lat.max                       geometry
#> 1 35.13501 140.0097 36.14069 POLYGON ((138.7776 36.12427...
#> 2 35.13501 140.0097 36.14069 POLYGON ((138.7776 36.12427...

衛星データの取得

ここでちょっとハードルが高いのは、CDSE には衛星データをそのままダウンロードする API はなく、OpenEO のように「衛星データに何らかの処理をしてその結果をダウンロードする」ということをする必要があります。その「何らかの処理」は evalscript という独自言語？で記述します。

よくある処理は、CDSE パッケージに同梱されているほか、以下のドキュメントにもまとまっています。

これを参考に、バンド2～4を65535諧調の値でダウンロードするスクリプトを書いたものが以下です。

script <- r"(
function setup() {
  return {
    input: [{
      bands: ["B02", "B03", "B04"]
    }],
    output: {
      id: "default",
      bands: 3,
      sampleType: SampleType.UINT16
    }
  }
}

function evaluatePixel(sample) {
    return [ sample.B04 * 65535, sample.B03 * 65535, sample.B02 * 65535 ]
}
)"

上のスクリプトを script 引数に指定してデータを取得するのが以下です。結果は terra パッケージの SpatRaster になっています。

image <- GetImage(
  bbox = c(139.7081, 35.4818, 139.8620, 35.6594),
  time_range = catalog$acquisitionDate,
  script = script,
  collection = "sentinel-2-l2a",
  format = "image/tiff",
  pixels = 2500,
  mosaicking_order = "leastCC", # 雲量が最も少ないシーンから値を取る
  client = OAuthClient
)

image
#> class       : SpatRaster 
#> dimensions  : 2500, 2500, 3  (nrow, ncol, nlyr)
#> resolution  : 6.156e-05, 7.104e-05  (x, y)
#> extent      : 139.7081, 139.862, 35.4818, 35.6594  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : file897c28a5e96 
#> colors RGB  : 1, 2, 3 
#> names       : file897c28a5e96_1, file897c28a5e96_2, file897c28a5e96_3

ちなみに、ダウンロードする画像の解像度は、pixels か resolution で指定しますが、resolution の挙動がよくわかりません。ドキュメントによれば単位はメートルらしいのですが、この領域は 10km くらいあるので、10m x 10m のタイルなら次元は 10,000 くらいになるはずです。しかし、実際には以下のような結果になります...。ということで、実質 pixels の上限値 2500 が解像度の上限になっているみたいです。

GetImage(
  ...
  resolution = 10,
  ...
)
#> class       : SpatRaster 
#> dimensions  : 1970, 1395, 3  (nrow, ncol, nlyr)
#> resolution  : 0.0001103226, 9.015228e-05  (x, y)
#> extent      : 139.7081, 139.862, 35.4818, 35.6594  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : file897c74ef54ec 
#> colors RGB  : 1, 2, 3 
#> names       : file897c74ef54ec_1, file897c74ef54ec_2, file897c74ef54ec_3

3-1-6 対象領域の画像を表示する

値の分布を確認

terra の使い方をまだよくわかっていないのですが、とりあえず各バンドの値の分布を見るだけなら hist() を使えばすぐ見れます。この謎のファイル名は、たぶんダウンロードして一時ファイルとして保存しているのが見えているだけで、それぞれ赤、青、緑、のはずです。

hist(image)

値の範囲を変えてプロット

本と同じく、0～3000 の範囲の値を 0～255 にスケールするのは clamp() と割り算を使って以下になるはずです。が、なんか見た目がだいぶ違うので、あってるか不安です...

plotRGB(clamp(image, 0, 3000) / 3000 * 255)

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

3-1-7 SpatioTemporal Asset Catalog (STAC) を利用する

R から STAC を使うには、rstac というパッケージが用意されています。STAC の公式サイトでも使い方が紹介されているので信用してよさそうです。

コレクションの取得

まず、rstac::collections() で衛星データ一覧を見て、Sentinel-2 L2A があるのを確認します。Element 84 の Earth Search は、本では v0 を使っていましたが、今は v1 があるようなのでそちらを使います。

stac_source <- rstac::stac("https://earth-search.aws.element84.com/v1")

available_collections <- stac_source |> 
    rstac::collections() |> 
    rstac::get_request()

available_collections
#> ###Collections
#> - collections (9 item(s)):
#>   - sentinel-2-pre-c1-l2a
#>   - cop-dem-glo-30
#>   - naip
#>   - cop-dem-glo-90
#>   - landsat-c2-l2
#>   - sentinel-2-l2a
#>   - sentinel-2-l1c
#>   - sentinel-2-c1-l2a
#>   - sentinel-1-grd
#> - field(s): collections, links, context

衛星データの検索

検索は rstac::stac_search() です。本とは datetime の指定の仕方が変わっていて、RFC 3339 形式で時刻まで指定しないといけないみたいです。

result <- stac_source |>
  rstac::stac_search(
    collections = "sentinel-2-l2a",
    datetime = "2021-04-10T00:00:00Z/2021-06-11T00:00:00Z",
    bbox = c(139.7081, 35.4818, 139.8620, 35.6594),
    limit = 100
  ) |>
  rstac::get_request()

result
#> ###Items
#> - matched feature(s): 21
#> - features (21 item(s) / 0 not fetched):
#>   - S2A_54SUE_20210606_2_L2A
#>   - S2A_54SUE_20210606_0_L2A
#>   - S2B_54SUE_20210601_1_L2A
#>   - S2B_54SUE_20210601_0_L2A
#>   - S2A_54SUE_20210527_1_L2A
#>   - S2A_54SUE_20210527_0_L2A
#>   - S2B_54SUE_20210522_1_L2A
#>   - S2B_54SUE_20210522_0_L2A
#>   - S2A_54SUE_20210517_1_L2A
#>   - S2A_54SUE_20210517_0_L2A
#>   - ... with 11 more feature(s).
#> - assets: 
#> aot, aot-jp2, blue, blue-jp2, coastal, coastal-jp2, granule_metadata, green, green-jp2, nir, nir-jp2, nir08, nir08-jp2, nir09, nir09-jp2, red, red-jp2, rededge1, rededge1-jp2, rededge2, rededge2-jp2, rededge3, rededge3-jp2, scl, scl-jp2, swir16, swir16-jp2, swir22, swir22-jp2, thumbnail, tileinfo_metadata, visual, visual-jp2, wvp, wvp-jp2
#> - item's fields: assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

サムネイルの取得

上で result を print した中に assets という項目がありますが、thumbnail もそこに入っているのがわかります。これらのデータを取得するのには 2 つ方法があります。

ひとつは、rstac::assets_download() を使ってデータをダウンロードすることです。出力先のディレクトリを指定するとそこに保存してくれます。

outdir <- "stac_data"
dir.create(outdir, showWarnings = FALSE)
rstac::assets_download(result, "thumbnail", outdir, overwrite = TRUE)

もうひとつは、rstac::assets_url() で URL を取得して、それを他の関数に渡すことです。ドキュメント的には、これは GDAL の virtual file system interface (/vsicurl とか)に渡すためのもののようですが、まあ普通にダウンロードしてもいいでしょう。

urls <- rstac::assets_url(result, "thumbnail")
curl::multi_download(urls, destfiles)

被雲率が最も低いものを探す

rstac::stac_search() の結果にはいろいろな情報があります。list になっていますが、rstac::items_as_tibble() とか rstac::items_as_sf() でいったんデータフレームにしてから絞り込むのが扱いやすそうです。

result |>
  rstac::items_as_tibble() |> 
  dplyr::filter(`eo:cloud_cover` < 10)

あるいは、rstac::items_filter() で直接絞り込むこともできます。

result |> 
  rstac::items_filter(filter_fn = \(x) x$properties$`eo:cloud_cover` < 10)

ちなみに、検索の際にも被雲率とかいろいろ条件を付けることができるみたいなので、もうちょっと書き方を勉強したいところです。

バンド2～4のデータをダウンロードする

基本的にはサムネイルと同じなので省略しますが、数分悩んだのは、B02 みたいな名前がどこにもないことでした。これは red・blue・green というわかりやすい名前になっているみたいでした。

# このダウンロードはちょっと重いので気を付ける
urls <- result |>
  rstac::items_filter(filter_fn = \(x) x$properties$`eo:cloud_cover` < 10) |> 
  rstac::assets_url(c("red", "blue", "green"))

『Python で学ぶ衛星データ解析基礎』を R でやってみる (1) 3章前半まで

本

3-1-4～3-1-5 APIでSentinel-2の画像を取得する

準備

コレクションの確認

衛星データの検索

衛星データの取得

3-1-6 対象領域の画像を表示する

値の分布を確認

値の範囲を変えてプロット

3-1-7 SpatioTemporal Asset Catalog (STAC) を利用する

コレクションの取得

衛星データの検索

サムネイルの取得

被雲率が最も低いものを探す

バンド2～4のデータをダウンロードする

続き

Discussion