Go to Qiita Advent Calendar 2024 Top
LoginSignup
4
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

R言語Advent Calendar 2024

Day 17
@ade_san

RからSentinel Hub APIをつかってSentinel-2衛星画像を取得する

Posted at

この記事はR言語 Advent Calendar 2024の17日目の記事です。まもなく冬至ということで、空に関する内容ということになりました（という後づけの理由）

Sentinel-2の衛星データを自動的にダウンロードしたい

  • 仕事で衛星画像の解析をすることがあり、せっかくなのでRでダウンロードから解析まで一気にやってしまおうとなりました。
  • Sentinel-2はESA(European Space Agency; 欧州宇宙機関)が運用している衛星
    • 可視光、近赤外光を含む複数の波長を宇宙から観測している
    • Copernicus Hubを通じて無償で提供されており、誰でも利用することが可能
      • 同サイトからはESAの他の衛星画像も利用可能（Sentinel-1)
  • Sentinel-2の衛星画像の取得手段は…
  • SentinelのAPIの使い方
    • Catalogue API, openEO API, Sentinel Hub APIと大きく分けられる
    • Copernicus Data Space Ecosystemにてアカウントを作成したら利用可能
      • ただし無償で利用可能なのは30000クレジット/月まで。有料分のAPIはSentinel Hubから購入することになりそう
    • 実は2023年10月まではCopernicus Open Access Hubから取得できていたのだが、このサイトが閉鎖され、現在の公式サイトとなっている。そのため、このサイトにアクセスするAPIは利用できなくなっている
    • そこで、新たなAPIに対応したRパッケージを探すことにした
  • RでAPIを利用した衛星画像の取得
    • {openeo}パッケージ: openEOサーバーを操作するためのパッケージ。openEOはESAが提供する地理空間データを解析するためのプラットフォーム
    • {CDSE}パッケージ: 今回紹介するパッケージ。
      • Sentinel Hub APIを用いて過去の衛星画像を検索、ダウンロードすることができるパッケージ
      • GeoTiff形式や、ただのJPEG形式でも取得できるほか、NDVIの形でも取得できる
      • また、タイルのうちさらに必要な箇所だけダウンロードできる。Copernicus Browserだとタイル単位でのダウンロードになるため、1枚1GB必要となることもある
    • ちなみに{openeo}パッケージはこの記事を書いている最中に見つけました。探してても見つからないのでなぜか突然みつけてしまう

API利用のための準備

  • 公式ドキュメントのBefore Startを見つつ軽く説明
    • Copernicus Data Space Ecosystemからアカウントを作成して、APIのトークンを取得。このトークンは一度しか表示されないのできっちり保存
    • トークンとユーザーIDは間違えて公開しないよう、Rスクリプトに平文で書くのではなく、.Renvironに保存して、Sys.getenv()で読み込むことが推奨されている
  • {sf}パッケージを利用するため、こちらも合わせてインストールをする

{CDSE}パッケージで何ができるの？

  • GetCollections()関数：{CDSE}パッケージにて利用可能なデータ一覧を確認
    • sentinel-2はlevel1Cと2Aが利用可能。
      • 両者の違いは値がどの程度補正がされているかによる。
      • 2Aは'Bottom of the atmosphere (BOA) reflectance, processed from L1C with Sen2Cor'ということで、大気下の反射率を表すように大気補正がされている。
    • その他、Sentinel-1, Sentinel-3も利用可能
  • SearchCatalog()関数: カタログから検索条件に合う衛星画像を検索。

    • 欲しい地域の範囲、時期から利用可能な衛星データのメタデータを取得する

    • 重要な情報として、tileCloudCoverという要素がある

      • 衛星画像のタイル（100km×100km）内に雲が何割含まれているか（雲量）を表している。この値が小さいほど雲が少なく、地表面が写っている割合が高いことを意味している
      • 事前に雲量が少ない画像のみを絞り込むことで、使える画像だけを取得することができる

    • Pasted image 20241129160240.png

    • SeasonalFilter()関数: 取得したカタログから特定の時系列のみを抽出する関数（今回は試していない…）

  • GetImage()関数: 文字通り画像を取得する関数
    • 引数の意味
      • format: "image/tiff", "image/png", "image/jpeg"の3つから選べる
      • file: この引数を設定すれば保存することができる。ファイル名指定はhere()を使う形でやるとやりやすい。
      • mosaicking_order: the order in which tiles are overlapped from which the output result is mosaickedらしい。？？
      • resolution(or pixels): 空間解像度の指定。resolutionならm単位、pixelsなら1-2500までの範囲。細かくするというより10mよりぼかしたいときにつかう？
      • buffer: width of the buffer to retrieve the image of enlarged areaらしい。？？
      • mask: the image should contain only pixels within Area of Interestらしい。？？
      • script: evalscriptで書かれたもので、どの波長のデータを取得し、どのような加工をするのかを指定したもの。{CDSE}パッケージではscriptsフォルダにいくつかのスクリプトが予め入っている。ちなみに自力で作成することも可能
    • 注意: 指定できるピクセル数は2500 pixelsまで
  • GetImageBy*関数群: 画像の取得をバッチ処理
    • GetImageByTimerange, GetImageByAOI, GetImageByBboxの3パターン。それぞれ時間、関心領域、バウンディングボックスでバッチ処理の範囲を決められる
    • サンプルコードでは、parallelパッケージを使った並列計算を実施している。これはapply系の並列計算になる様子(当たり前）
      • あまりapply系に親しみがないので今回は{purrr}パッケージのmap関数を使ってみる

実際にSentinel-2の衛星画像を取得する by map関数

## 画像の保存フォルダ名を指定。もし存在しない場合はフォルダを作成
folderName <-
  here("analysis", "processing", "img_sentinel2")
  
if (!dir.exists(folderName)) {dir.create(folderName)}

## 撮影日のうち、取得したい期間のみの撮影日のベクトルを取得
days <-
  imgs |>
    filter(between(acquisitionDate, ymd("2024-04-01"), ymd("2024-04-15"))) %>%
    .$acquisitionDate
 
## 保存先の名前一覧を作成する
savePathList <-
  tibble(
    days = days
  ) |>
    mutate(
      saveName = paste0("sentinel2_", days, ".tiff"),
      savePath = here("analysis", "img_sentinel2", saveName)
    )

## 衛星画像から取得したい内容を選択。CDSEパッケージに初期からあるevalscriptの一覧を取得
scriptList <-
  dir(
    system.file("scripts", package = "CDSE"), pattern = "\\.js$",
    full.names = T
    )
script_file <- scriptList[5]　## Sentinel2のLevel2Aのrowimageを取得するスクリプト

## 衛星画像を取得し、ラスタ形式(geotiff)で保存
map2(
  savePathList[[1]], savePathList[[3]],
  ~ GetImageByTimerange(
    time_range = .x, ## 取得したい時期
    file = .y, ## 保存先のファイル名
    aoi = aoi, ## 取得した領域のgeojson
    collection = "sentinel-2-l2a", ## 取得したいデータ指定:Sentinel-2の2Aレベルのもの
    script = script_file,
    format = "image/tiff", ## GeoTiffで取得
    resolution = 10, ## m単位での解像度指定:10m、Sentinel-2の最高解像度
    mosaicking_order = "leastCC",
    mask = TRUE,
    buffer = 100,
    client = OAuthClient
  )
)
  • 取得したい領域をGeoJson形式で取得→sf::read_sf()関数でその領域をsfオブジェクト化する
  • ちなみに同じ日付の衛星画像を取得しようとするとエラーが生じる（上書きができない）。たまにカタログに同一日の画像が2枚あることがある。その際は日付の一覧にunique()などで、重複がないようにしたらとりあえず上書きが発生しない
  • GetImageByTimerange()の引数について
    • mosaickingOrder:モザイク処理される重複タイルの順序の設定。ほとんどの場合でタイルは基本同じ軌道、取得のもの
      • mostResentがデフォルト。同じタイムスタンプのタイルが複数ある時は、SciHubからダウンロードされたもののうち後のもの（より現在に近い？もの）を使う
      • leastCCは雲量が最も少ないタイルから選択される
  • 北海道のある地域の衛星画像　Pasted image 20241129162928.png

Sentinel-2衛星画像についての補足

  • Sentinel Hubのサイトを元にまとめている

  • Sentinel-2のバンド構成

    • B1:aerosol(443nm, 60m)
    • B2:blue(490nm, 10m)
    • B3:green(560nm, 10m)
    • B4:red(665nm, 10m)
    • B5:red edge(705nm, 20m)
    • B6:(740nm, 20m)
    • B7:(783nm, 20m)
    • B8:NIR(842nm, 10m)
    • B8A:(865nm, 20m)
    • B9:(945nm, 60m)
    • B10:cirrus cloud detectionにも(1375nm, 60m)←大気圏より下の情報が含まれていないためL2Aからは除外
    • B11:SWIR1(1610nm, 20m)
    • B12:SWIR2(2190nm, 20m)

  • その他の取得可能な情報

    • "Sen2Cor"が大気補正用にESAが使っているもの。Sentinel-2は大気補正されたものが利用可能
    • AOT: "Aerosol Optical Thickness map, based on Sen2Cor processor"、エアロゾルの厚さを計算したもの
    • SNW(20m): "Snow probability, based on Sen2Cor processor"。雪である確率
    • CLD(20m): "Cloud probability, based on Sen2Cor processor" 雲である確率
    • SCL(20m): "Scene classification data, based on Sen2Cor processor, codelist"。以下の12種類の分類結果をマップにしたもの。このうち、7~10番が雲の領域を表している ![[Pasted image 20241210092447.png|300]]
    • 左の衛星画像は、SCLでは右のような地図となる。黄色、緑色の領域が雲の領域　Pasted image 20241217200653.png

  • その他の用語について

    • ダウンサンプリング/アップサンプリング
      • upsampling=ソースの解像度よりも高解像度化する場合（拡張と呼ばれる？）, downsampling=解像度を下げる場合。どちらも3つの処理方法がある。デフォルトがNEAREST, 他にはBILINEAR, BICUBICがある。
      • upsamplingのイメージがリンクされている
    • 単位：波長(Reflectance)とDN
      • DN→ソースデータのピクセルの実際の値。DNは難かしらの物理量を表すときのみ提供される。Sentinel-2では, DN=10000×REFLECTANCEがデフォルト
      • 非DN単位の値（？）は、ソース（DN）から32ﾋbit浮遊値から計算。非線形変換である可能性がありバンドの完全な値位置やステップサイズの予測は困難。
    • RawBandsのときのharmonizeの留意点
      • https://documentation.dataspace.copernicus.eu/APIs/SentinelHub/Data/S2L1C.html#harmonize-values
      • digital number(DN)の解釈が2022年1月のversion04.00以前から変更されているらしい
      • それぞれHarmovize valueを有効にすると…
        • REFLECTANCE: 負の値は0に打ち切られる。
        • DN: アップデート以前のデータと比較できるよう、DN = 10000 * REFLECTANCEと計算される。一方、無効になるとソースファイルに記載の通りの値になるので、処理時のベースラインが異なれば値も変化する。モザイク処理時に注意せよとのこと（画像をまたぐときのモザイク処理とか？）

evalscriptって何？

  • Sentinel Hubでデータの取得、前処理などを指定するためのスクリプト。javascript形式で書かれている
  • "RawBands.js"
  • "S2L2A_specifics.js
    • AOT, SCL, SNW, CLDと太陽角、衛星の角度といった基本情報。FLOAT32型として取得。
    • ちなみにSCLは本来UINT8型であるため、FLOAT32型だと小数点になる恐れがある（事実QGISでは小数点となっている）ものの、RではちゃんとUINT8型になっていた。
  • その他にあるのは…
    • "NDVI_CLOUDS_STAT.js"
    • "NDVI_dataMask_float32.js"
    • "NDVI_float32.js"
    • "NDVI_uint8.js"
    • "TrueColor.js"
    • "TrueColorS2L2A.js"
  • 自分でjavascriptが書ければ必要な要素だけ取得したり、変換をサーバー側でやってしまうことも可能ということ

まとめ

  • Sentinel-2の衛星データをRで取得するためのパッケージとして{CDSE}パッケージが便利
  • evalscriptを自分で書くことで、取得するデータや処理方法を指定することも可能
4
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?