💨

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

2025/05/14に公開

GIS

衛星データ

tech

本

前回

画像のダウンロード

本と同じ条件で検索するとこんな感じです。

# ここでは直接 bbox を返すだけでいいが、本と合わせてポリゴンを返すようにする
sel_square <- function(lon, lat, delta_lon, delta_lat) {
  poly <- sf::st_polygon(list(
    rbind(
      c(lon + delta_lon, lat + delta_lat),
      c(lon + delta_lon, lat - delta_lat),
      c(lon - delta_lon, lat - delta_lat),
      c(lon - delta_lon, lat + delta_lat),
      c(lon + delta_lon, lat + delta_lat)
    )
  ))

  sf::st_sfc(poly, crs = 4326)
}

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

aoi <- sel_square(130.53, 31.60, 0.06, 0.02)

items <- stac_source |>
  rstac::stac_search(
    collections = "sentinel-2-l2a",
    datetime = "2019-01-01T00:00:00Z/2020-12-31T23:59:59Z",
    bbox = sf::st_bbox(aoi),
    limit = 100
  ) |> 
  rstac::get_request()

ここで、本では eo:cloud_cover < 10 というのを検索条件に入れているのですが、やってみても効きませんでした。rstac がおかしいのかな？と思ったのですが、curl や rustac などほかのツールで試してもだめだったので、諦めました。/queryable というエンドポイントで CQL2 中で使える属性が確認できるのですが、これ見てもいける感じするんだけどなあ...

/queryables の結果

httr2::request("https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/queryables") |> 
  httr2::req_perform() |> 
  httr2::resp_body_json()
#> $`$schema`
#> [1] "https://json-schema.org/draft/2020-12/schema"
#> 
#> $`$id`
#> [1] "https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/queryables"
#> 
#> $type
#> [1] "object"
#> 
#> $title
#> [1] "Queryables for Collection sentinel-2-l2a"
#> 
#> $properties
#> $properties$`eo:cloud_cover`
#> $properties$`eo:cloud_cover`$`$ref`
#> [1] "https://stac-extensions.github.io/eo/v1.0.0/schema.json#/definitions/fields/properties/eo:cloud_cover"
#> 
#> 
#> 
#> $additionalProperties
#> [1] TRUE

上の条件だと広すぎて 400 件以上がヒットするのですが、 rstac::stac_search() の limit を上げるとエラーになってしまいます。仕方ないので、ダウンロードできた 100 件の中から条件がいいものを選ぶことにします。

# eo:cloud_cover の昇順に並び替えて、いちばん雲がないものを選ぶ
items |> 
  rstac::items_as_tibble() |> 
  arrange(`eo:cloud_cover`) |> 
  glimpse()

# URL を調べる
urls <- items |> 
  rstac::items_filter(filter_fn = \(x) x$properties$`s2:product_uri` == "S2B_MSIL2A_20201112T015919_N0500_R060_T52SFA_20230414T013628.SAFE") |> 
  rstac::assets_url(c("swir16", "nir", "red", "green", "blue"))

dir.create("data", showWarnings = FALSE)        # 保存ディレクトリ
destfiles <- file.path("data", basename(urls))  # ファイル名
curl::multi_download(urls, destfiles)           # ダウンロード

NDVI、MNDWI の計算

ダウンロードしたファイルを読み込みます。ちなみに、前回まで気付いてなかったんですが、rast() には複数のファイルを渡せます。このうち、バンド11は解像度が異なるので後回しにします。

library(terra)

tiffs <- rast(c(
  "data/B08.tif", 
  "data/B04.tif",
  "data/B03.tif",
  "data/B02.tif"
))

tiffs
#> class       : SpatRaster 
#> dimensions  : 10980, 10980, 4  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 6e+05, 709800, 3490200, 3600000  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 52N (EPSG:32652) 
#> sources     : B08.tif  
#>               B04.tif  
#>               B03.tif  
#>               B02.tif  
#> names       : B08, B04, B03, B02

トゥルーカラーでプロットするとこんな感じ。ちなみに、説明してませんでしたが、値を 0～255 にスケールし直すのは、stretch() を使わなくても geom_spatraster_rgb() に stretch という引数があります。ただ、minq・maxq にあたる引数はないので、細かく調整したい場合はやはり自分で stretch() で指定するのがいいでしょう。

library(tidyterra)
library(ggplot2)

ggplot() +
  geom_spatraster_rgb(data = tiffs[[2:4]], stretch = "lin")

出典：Element 84 提供の Sentinel-2 の衛星画像データを使って描画

これを鹿児島市周辺で切り抜きます。行政区画のデータを国土数値情報からダウンロードしてきます。

kagoshima <- sf::st_read("N03-20240101_46_GML/N03-20240101_46.geojson")

N03_004 に市町村名が入っているので、鹿児島市のみに絞り込みます。それを st_bbox() で bbox に変換、さらに bbox を st_as_sfc() で地物に変換します。また、切り抜くためには CRS を合わせないといけないので st_transform() で投影変換します。

aoi <- kagoshima |> 
  dplyr::filter(N03_004 == "鹿児島市") |> 
  sf::st_bbox() |> 
  sf::st_as_sfc() |> 
  sf::st_transform(crs(tiffs))

まとめて切り抜いておきます。

tiffs <- tiffs |> 
  crop(aoi)

NDVI の計算はこうです。1つにまとめたラスタから特定のバンドだけを取り出すには [[ ]] でアクセスできます。

ndvi <- (tiffs[[1]] - tiffs[[2]]) / (tiffs[[1]] + tiffs[[2]])

ggplot() +
  geom_spatraster(data = ndvi) +
  scale_fill_whitebox_c(palette = "purple")

出典：Element 84 提供の Sentinel-2 の衛星画像データを使って描画

MNDWI の計算をするために、バンド11を読み込みます。上に書いたように、バンド11は解像度が他と異なっています。このままでは計算できないので、リサンプリングが必要です。リサンプリングには resample() を使います。 method にリサンプリング方法がいろいろ指定できますが、ここでは本に合わせて average を使っておきます。

swir <- rast("data/B11.tif") |> 
  crop(aoi)

swir_resampled <- resample(swir, tiffs, method = "average")

以下が計算結果です。

mndwi <- (tiffs[[3]] - swir_resampled) / (tiffs[[3]] + swir_resampled)

ggplot() +
  geom_spatraster(data = mndwi) +
  scale_fill_whitebox_c(palette = "purple")

出典：Element 84 提供の Sentinel-2 の衛星画像データを使って描画

本

前回

画像のダウンロード

NDVI、MNDWI の計算

Discussion