🌱

RによりGoogle Earth Engineの植生指数の地図・時系列データを表示する

2024/02/25に公開

googleearthengine

衛星画像

landsat

tech

変更履歴・概要

2024/2/25 初版

できるようになること

RからGoogle Earth Engine（GEE）を操作して，衛星データを植生指数に変換して表示します。あわせて，標高データを利用して，海の領域をマスクしたり，陰影起伏図と重ねたりした画像も作成してみます。また，植生指数のデータを時系列データとして出力し，グラフとして表示してみます。

初期設定

RでGEEを操作する設定は，RでGoogle Earth Engineを操作できるようにするを，光学衛星データを表示する方法は，RでGoogle Earth Engineの光学衛星データを表示するを，それぞれ参考にしてください。

library(rgee)
library(rgeeExtra)
library(tidyverse)
library(leaflet)
ee_Initialize(project = 'プロジェクト名')
eemont <- reticulate::import("eemont")

標高データの利用

植生指数にふれる前に，標高データを利用して，海の領域のマスクや，陰影起伏図と重ねる処理の方法を紹介します。

海の領域をマスクした標高データの表示

定番はCGIAR/SRTM90_V4で以下のコードで海の領域でマスクした標高データを表示できます。

geom_point <- c(135.1830778, 34.69125701) # 兵庫県庁
ee_geom <- ee$Geometry$Point(geom_point)
SRTM90_V4 <- ee$Image('CGIAR/SRTM90_V4')
landMask <- SRTM90_V4$mask()
SRTM90_V4 <- SRTM90_V4$updateMask(landMask)
DSMVis <- list(
  min =0, max=1000, palette = terrain.colors(10)
)
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 11)
Map$addLayer(SRTM90_V4, DSMVis, 'Elevation')

ただ，データ取得年が2000年と少し古く，一部の領域では現状と一致していないようです。そこで，JAXAが整備しているALOS全球数値地表モデル (DSM) “ALOS World 3D - 30m (AW3D30)”を利用してみます。このデータセットはImageCollectionなので，データ表示にはmosaic()によりImageにしておきます。

geom_point <- c(135.1830778, 34.69125701) # 兵庫県庁
ee_geom <- ee$Geometry$Point(geom_point)
AW3D30 <- ee$ImageCollection('JAXA/ALOS/AW3D30/V3_2')
AW3D30.DSM <- AW3D30$select('DSM')
AW3D30.MSK <- AW3D30$select('MSK')$mosaic()
seaBitMask <- bitwShiftL(1, 1)
seamask <- AW3D30.MSK$bitwiseAnd(seaBitMask)$eq(0)
DSMVis <- list(
  min =0, max=1000, palette = terrain.colors(10)
)
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 11)
Map$addLayer(AW3D30.DSM$mosaic()$updateMask(seamask), DSMVis, 'Elevation') %>% 
  leaflet::addTiles(urlTemplate = "", attribution = '| Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

さらに海の領域のみを別の単色にしてみます。

sea <- seamask$Not()
sea <- sea$mask(sea)
DSM.notsea <- ee$ImageCollection(list(
  AW3D30.DSM$mosaic()$updateMask(seamask)$visualize(min =0, max=1000, palette = terrain.colors(10)),
  sea$visualize(palette = "#000044")
))$mosaic()
Map$addLayer(
  DSM.notsea,
  {},
  "Elevation"
) %>% leaflet::addTiles(urlTemplate = "", attribution = '| Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

陰影起伏図の作成

AW3D30データを用いて陰影起伏図を作成してみます。ee$Terrain$hillshade()を利用しますが，入力データはImageです。mosaic()でImageにしますが，CRSも変わってしまうため，以下のとおり元のデータのCRSを設定しておく必要があります。

proj <- AW3D30.DSM$first()$select(0)$projection()
AW3D30.shade <- ee$Terrain$hillshade(AW3D30.DSM$mosaic()$setDefaultProjection(proj))
Map$addLayer(AW3D30.shade, name = "hillshade") %>% 
  leaflet::addTiles(urlTemplate = "", attribution = '| Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

続いて，海の領域を単色化し，標高データと重ね合わせてみます。

DSMshade.notsea <- ee$ImageCollection(list(
  AW3D30.DSM$mosaic()$updateMask(seamask)$visualize(min =0, max=1000, palette = terrain.colors(10)),
  AW3D30.shade$updateMask(0.15)$visualize(palette=c("#000000", "#FFFFFF")),
  sea$visualize(palette = "#000044")
))$mosaic()
Map$addLayer(
  DSMshade.notsea,
  {},
  "ElevationShade"
) %>% leaflet::addTiles(urlTemplate = "", attribution = '| Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

植生指数の地図表示

rgeeExtraライブラリもしくはeemontで扱える植生指数はAwesome spectral Indicesにリストアップされています。

まずは，データの読み込みと前処理を行っておきます。今回はrgeeExtraライブラリで利用可能なee_ImageCollection_closestを利用して，指定日の直近のデータを抽出してみます。

L8_SR <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
  filterBounds(ee_geom) %>%
  ee_ImageCollection_closest("2019-09-01", 1, "Month") %>%
  ee$Image$scaleAndOffset()

続いて，植生指数への変換。rgeeExtraライブラリのspectralIndex()を利用します。もしくはeemontのee$Image$specralIndices()でも同じことが可能です。

L8_SR_SpectralIndex <- L8_SR %>%
  ee$Image$maskClouds() %>%
  ee$Image$spectralIndex(c('NDVI','SAVI'))%>%
  '[['(1)

地図として表示してみましょう。

# https://github.com/csaybar/rgee/blob/examples/Visualization/ndvi_symbology.R
visPalette <- c(
  "#d73027", "#f46d43", "#fdae61",
  "#fee08b", "#d9ef8b", "#a6d96a",
  "#66bd63", "#1a9850"
)
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 11)
map_index <- Map$addLayer(
  L8_SR_SpectralIndex[['NDVI']],
  list(min=0.0, max=1.0, palette=visPalette), 
  "NDVI"
) + 
  Map$addLegend(
    list(min=0.0, max=1.0, palette=visPalette),
    position="bottomright",color_mapping = "numeric") 
map_index %>% 
  leaflet::addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')

植生指数データとして，一定期間の中央値を計算することもできます。

geom_point <- c(135.1830778, 34.69125701) # 兵庫県庁
ee_geom <- ee$Geometry$Point(geom_point)
L8_SR <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
  filterDate("2019-01-01", "2020-01-01")$
  filterBounds(ee_geom)$
  filter(ee$Filter$lt('CLOUD_COVER', 50))$
  preprocess()
L8_SR_SpectralIndex <- L8_SR %>%
  ee$ImageCollection$spectralIndex(c('NDVI','SAVI') )%>%
  ee$ImageCollection$median()

標高データとの重ね合わせ

植生指数データのうち，海の領域をマスクしてみます。

L8_SR_SpectralIndex_seaMask <- L8_SR_SpectralIndex$updateMask(seamask)
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 11)
Map$addLayer(
  L8_SR_SpectralIndex_seaMask[['NDVI']],
  list(min=0.0, max=1.0, palette=visPalette), 
  "NDVI"
) %>% leaflet::addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey | Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

海の領域のみを別の単色にし，陰影起伏図を重ねてみます。

SI_notsea <- ee$ImageCollection(list(
  L8_SR_SpectralIndex$visualize(bands=c('NDVI'), min = 0.0, max=1.0,
                                palette = visPalette),
  AW3D30.shade$updateMask(0.15)$visualize(palette=c("#000000", "#FFFFFF")),
  sea$visualize(palette = "#000044")
))$mosaic()
map_index <- Map$addLayer(
  SI_notsea,
  {},
  "custom mosaic"
) +
  Map$addLegend(
    list(min=0.0, max=1.0, palette=visPalette),
    position="bottomright",color_mapping = "numeric") 
map_index %>%
  leaflet::addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey | Credit: ALOS DSM: Global 30m v3.2 (JAXA)')

植生指数のグラフ化

つづいて，植生指数のデータを時系列で抜き出して，グラフとして表示してみます。

抜き出す地点は兵庫県庁と六甲山の2地点とします。

pref_office <- c(135.1830778, 34.69125701) # 兵庫県庁
mt_rokko <- c(135.263611, 34.778056) # 六甲山

Google Earth Engineで読み込める形式に変換します。なお，各地点で100mのバッファーを発生させています。

ee_geom <- ee$FeatureCollection(list(
  ee$Feature(ee$Geometry$Point(pref_office)$buffer(100), list(name="pref_office")),
  ee$Feature(ee$Geometry$Point(mt_rokko)$buffer(100), list(name="mt_rokko"))
))

Landsat8の衛星データを読み込んで前処理。

L8_SR <- ee$ImageCollection('LANDSAT/LC08/C02/T1_L2')$
  filterBounds(ee_geom)$
  filterDate('2014-01-01','2024-01-01')$
  preprocess()

植生指数としてNDVIを算出します。

L8_SR_SpectralIndex <- L8_SR %>%
  ee$ImageCollection$spectralIndex(c('NDVI','SAVI') )

データの抽出は，eemontのgetTimeSeriesByRegions()を使います(https://eemont.readthedocs.io/en/latest/guide/timeSeries.html)。ここでは，上で発生させた各地点100ｍのバッファー内の平均値を求めています。地点が1つの場合や，複数地点をまとめて1つのデータとしたい場合はgetTimeSeriesByRegion()を使います。

L8_SR_SI <- L8_SR_SpectralIndex %>%
  ee$ImageCollection$getTimeSeriesByRegions(
    reducer = ee$Reducer$mean(),
    collection= ee_geom,
    bands = c('NDVI','SAVI'),
    scale = 30)

Rでデータ処理を行うため，ee_as_sf()を利用します（rgee-example）。合わせて，-9999はNAのためNAに変換しておきます。

L8_SR_SI_toR <- ee_as_sf(L8_SR_SI) %>%
  mutate_at(vars(c("NDVI", "SAVI")), na_if, -9999)
head(L8_SR_SI_toR)

## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 135.182 ymin: 34.69036 xmax: 135.2647 ymax: 34.77896
## Geodetic CRS:  WGS 84
##   NDVI SAVI                date        name reducer                       geometry
## 1   NA   NA 2014-01-12 01:35:59 pref_office    mean POLYGON ((135.1831 34.69216...
## 2   NA   NA 2014-01-12 01:35:59    mt_rokko    mean POLYGON ((135.2636 34.77896...
## 3   NA   NA 2014-01-28 01:35:52 pref_office    mean POLYGON ((135.1831 34.69216...
## 4   NA   NA 2014-01-28 01:35:52    mt_rokko    mean POLYGON ((135.2636 34.77896...
## 5   NA   NA 2014-02-13 01:35:37 pref_office    mean POLYGON ((135.1831 34.69216...
## 6   NA   NA 2014-02-13 01:35:37    mt_rokko    mean POLYGON ((135.2636 34.77896...

グラフとして表示。geom_rug()で欠損を含むデータ取得日を示しています。

ggplot(L8_SR_SI_toR, aes(x=date, y=NDVI, colour=name)) + 
  geom_smooth(linewidth=0.25, se = FALSE, span = 0.15) +
  geom_point(size=2) + 
  geom_rug(sides="b", colour="black") +
  ylim(c(0,1)) +
  scale_x_datetime(date_breaks="2 year", date_labels="%y/%m") +
  theme(legend.position="bottom")

年ごとに集計してみる

RでGoogle Earth Engineの光学衛星データをタイムラプス動画として表示するで紹介した，年ごとにデータを集計（減算）したデータの作成方法を利用して，植生指数データを年ごとに集計した上で，グラフ化してみます。

年ごとのデータ作成

years <- ee$List$sequence(2014, 2022)
createYearlyImage <- function(beginningYear){
  startDate <- ee$Date$fromYMD(beginningYear, 1, 1)
  endDate <- startDate$advance(1, 'year')
  yearFiltered <- L8_SR_SpectralIndex$filter(ee$Filter$date(startDate, endDate))
  total <- yearFiltered$median()$set('system:time_start', startDate$millis())
  return(total)
}
yearlyImages <- years$map(ee_utils_pyfunc(createYearlyImage))
yearlyCollection <- ee$ImageCollection$fromImages(yearlyImages)

2地点のデータの抽出

L8_SR_SI_yearly <- yearlyCollection %>%
  ee$ImageCollection$getTimeSeriesByRegions(
    reducer = ee$Reducer$mean(),
    collection = ee_geom,
    bands = c('NDVI'),
    scale = 30)
L8_SR_SI_yearly_toR <- ee_as_sf(L8_SR_SI_yearly) %>%
  mutate_at(vars(c("NDVI")), na_if,-9999)

グラフ化

ggplot(L8_SR_SI_yearly_toR, aes(x=date, y=NDVI, colour=name)) + 
  geom_point(size=2) + geom_line() +
  ylim(c(0,1)) +
  theme(legend.position="bottom")

Session info

sessionInfo()

## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Japanese_Japan.utf8  LC_CTYPE=Japanese_Japan.utf8    LC_MONETARY=Japanese_Japan.utf8
## [4] LC_NUMERIC=C                    LC_TIME=Japanese_Japan.utf8    
## 
## time zone: Etc/GMT-9
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] leaflet_2.1.2   lubridate_1.9.2 forcats_1.0.0   stringr_1.5.0   dplyr_1.1.2     purrr_1.0.1     readr_2.1.4    
##  [8] tidyr_1.3.0     tibble_3.2.1    ggplot2_3.4.2   tidyverse_2.0.0 rgeeExtra_0.0.1 rgee_1.1.5      here_1.0.1     
## [15] fs_1.6.2        rmarkdown_2.23 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0        farver_2.1.1            lazyeval_0.2.2          fastmap_1.1.1          
##  [5] webshot2_0.1.1          promises_1.2.0.1        digest_0.6.32           timechange_0.2.0       
##  [9] lifecycle_1.0.3         sf_1.0-12               geojsonio_0.11.2        ellipsis_0.3.2         
## [13] processx_3.8.1          terra_1.7-39            magrittr_2.0.3          compiler_4.3.1         
## [17] rlang_1.1.1             tools_4.3.1             utf8_1.2.3              yaml_2.3.7             
## [21] knitr_1.43              labeling_0.4.2          htmlwidgets_1.6.2       curl_5.0.0             
## [25] sp_2.0-0                classInt_0.4-9          reticulate_1.30         KernSmooth_2.23-21     
## [29] websocket_1.4.1         httpcode_0.3.0          withr_2.5.0             grid_4.3.1             
## [33] fansi_1.0.4             e1071_1.7-13            leafem_0.2.0            colorspace_2.1-0       
## [37] scales_1.2.1            crul_1.4.0              cli_3.6.1               crayon_1.5.2           
## [41] generics_0.1.3          rstudioapi_0.14         geojson_0.3.5           tzdb_0.4.0             
## [45] DBI_1.1.3               chromote_0.1.2          proxy_0.4-27            splines_4.3.1          
## [49] base64enc_0.1-3         vctrs_0.6.2             V8_4.3.3                webshot_0.5.4          
## [53] Matrix_1.5-4.1          jsonlite_1.8.4          geojsonsf_2.0.3         hms_1.1.3              
## [57] crosstalk_1.2.0         units_0.8-2             glue_1.6.2              codetools_0.2-19       
## [61] ps_1.7.5                leaflet.providers_1.9.0 stringi_1.7.12          gtable_0.3.3           
## [65] later_1.3.1             raster_3.6-20           munsell_0.5.0           pillar_1.9.0           
## [69] htmltools_0.5.5         jqr_1.2.3               R6_2.5.1                rprojroot_2.0.3        
## [73] evaluate_0.21           lattice_0.21-8          highr_0.10              png_0.1-8              
## [77] class_7.3-22            Rcpp_1.0.10             nlme_3.1-162            mgcv_1.8-42            
## [81] xfun_0.39               pkgconfig_2.0.3