🌱
RによりGoogle Earth Engineの植生指数の地図・時系列データを表示する
変更履歴・概要
2024/9/14 更新 実行環境をR4.3.3に更新しています。addLegendの引数を見直して，凡例のタイトルを表示できるようにしました。
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),
name='NDVI', position="bottomright", color_mapping = "numeric", opacity=1)
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),
name='NDVI', position="bottomright", color_mapping = "numeric", opacity=1)
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.3 (2024-02-29 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.2.2 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [8] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 rgeeExtra_0.0.1 rgee_1.1.7 here_1.0.1
## [15] fs_1.6.4 rmarkdown_2.26
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.16.0 jsonlite_1.8.8 magrittr_2.0.3 farver_2.1.1
## [5] vctrs_0.6.5 memoise_2.0.1 base64enc_0.1-3 terra_1.7-71
## [9] webshot_0.5.5 htmltools_0.5.8.1 usethis_2.2.3 curl_5.2.1
## [13] raster_3.6-26 KernSmooth_2.23-22 htmlwidgets_1.6.4 plyr_1.8.9
## [17] cachem_1.0.8 mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
## [21] webshot2_0.1.1 Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
## [25] shiny_1.8.1.1 digest_0.6.34 colorspace_2.1-0 ps_1.7.6
## [29] rprojroot_2.0.4 leafem_0.2.3 pkgload_1.3.4 crosstalk_1.2.1
## [33] geojson_0.3.5 labeling_0.4.3 fansi_1.0.6 timechange_0.3.0
## [37] httr_1.4.7 mgcv_1.9-1 compiler_4.3.3 proxy_0.4-27
## [41] remotes_2.5.0 bit64_4.0.5 withr_3.0.0 DBI_1.2.2
## [45] pkgbuild_1.4.4 highr_0.10 protolite_2.3.0 sessioninfo_1.2.2
## [49] classInt_0.4-10 tools_4.3.3 chromote_0.3.1 units_0.8-5
## [53] odbc_1.4.2 httpuv_1.6.15 glue_1.7.0 nlme_3.1-164
## [57] promises_1.3.0 grid_4.3.3 sf_1.0-16 geojsonio_0.11.3
## [61] generics_0.1.3 gtable_0.3.5 leaflet.providers_2.0.0 tzdb_0.4.0
## [65] class_7.3-22 websocket_1.4.1 hms_1.1.3 sp_2.1-4
## [69] utf8_1.2.4 pillar_1.9.0 later_1.3.2 splines_4.3.3
## [73] lattice_0.22-5 bit_4.0.5 tidyselect_1.2.1 miniUI_0.1.1.1
## [77] knitr_1.46 V8_4.4.2 crul_1.4.2 xfun_0.43
## [81] devtools_2.4.5 DT_0.33 stringi_1.8.3 ggmap_4.0.0
## [85] lazyeval_0.2.2 yaml_2.3.8 geojsonsf_2.0.3 evaluate_0.23
## [89] codetools_0.2-19 httpcode_0.3.0 cli_3.6.1 xtable_1.8-4
## [93] reticulate_1.35.0 munsell_0.5.1 processx_3.8.4 jquerylib_0.1.4
## [97] jqr_1.3.3 Rcpp_1.0.12 png_0.1-8 ellipsis_0.3.2
## [101] blob_1.2.4 jpeg_0.1-10 profvis_0.3.8 urlchecker_1.0.1
## [105] bitops_1.0-7 scales_1.3.0 e1071_1.7-14 crayon_1.5.2
## [109] rlang_1.1.3
Discussion