🌱
RでHansenGlobalForestChangeデータと雲除去したSentinel画像を表示してみる
変更履歴・概要
2024/11/12 初版
2024/11/19 Sentinel-2のデータ表示でRGBの設定順が逆になっていたので修正
できるようになること
Rを用いてGoogle Earth Engine（GEE）のHansen Global Forest Changeデータを操作します。あわせて背景画像として，Sentinel-2について，2023年10月にGoogleがリリースしたCloud Score+データセットを利用した雲除去を行った画像を利用できるようにします。
初期設定
Cloud Score+データセットの説明にあるlinkcollection関数を利用するため， earthengine-apiをconda update -c conda-forge earthengine-apiにより，0.1.370から1.1.0にアップデートしました。
rgeeライブラリで検証されたバージョンではなく，ee_Initialize()でee_check_root_folderに関する エラーがでてGEEの認証ができません。
https://github.com/r-spatial/rgee/issues/355
https://github.com/r-spatial/rgee/issues/370
そこで，reticulateライブラリを利用してpythonコードを直接実行します。
reticulate::py_run_string("import ee")
reticulate::py_run_string("ee.Initialize(project='プロジェクト名')")
これだけでは，rgeeライブラリにある関数のうち，RへGoogleDrive等を経由してデータを読み込ませるee_as_sf()などが機能しません。原因は，これらの関数がee_Initialize()により生成されるrgee_sessioninfo.txtを読み込もうとするのですが，上記のreticulateライブラリから直接認証させた状態ではtxtファイルが生成されていないことによるものです。以下のファイルをRスクリプトファイルとして保存した上で，source()により読み込みます。
ee_sessioninfo.R
## ee_sessioninfo.R
## rgeeライブラリで`rgee_sessioninfo.txt`を生成する関数のみを抽出したファイル
# ee_utils.R
ee_utils_py_to_r <- function(x) {
p_r <- suppressWarnings(try(reticulate::py_to_r(x), silent = TRUE))
if (class(p_r) %in% 'try-error') {
return(x)
} else {
return(p_r)
}
}
# utils-auth.R
ee_source_python <- function(oauth_func_path) {
module_name <- gsub("\\.py$", "", basename(oauth_func_path))
module_path <- dirname(oauth_func_path)
reticulate::import_from_path(module_name, path = module_path, convert = FALSE)
}
ee_connect_to_py <- function(path, n = 5) {
ee_utils <- try(ee_source_python(oauth_func_path = path), silent = TRUE)
# counter added to prevent problems with reticulate
con_reticulate_counter <- 1
while (any(class(ee_utils) %in% "try-error")) {
ee_utils <- try(ee_source_python(path), silent = TRUE)
con_reticulate_counter <- con_reticulate_counter + 1
if (con_reticulate_counter == (n + 1)) {
python_path <- reticulate::py_discover_config()
message_con <- c(
sprintf("The current Python PATH: %s",
bold(python_path[["python"]])),
"does not have the Python package \"earthengine-api\" installed. Do you restarted/terminated",
"your R session after install miniconda or run ee_install()?",
"If this is not the case, try:",
"> ee_install_upgrade(): Install the latest Earth Engine Python version.",
"> reticulate::use_python(): Refresh your R session and manually set the Python environment with all rgee dependencies.",
"> ee_install(): To create and set a Python environment with all rgee dependencies.",
"> ee_install_set_pyenv(): To set a specific Python environment."
)
stop(paste(message_con,collapse = "\n"))
}
}
return(ee_utils)
}
ee_check_init <- function() {
# if EARTHENGINE_PYTHON is defined, then send it to RETICULATE_PYTHON
earthengine_python <- Sys.getenv("EARTHENGINE_PYTHON", unset = NA)
if (!is.na(earthengine_python))
Sys.setenv(RETICULATE_PYTHON = earthengine_python)
# Check ee_utils.py sanity
ee_current_version <- system.file("python/ee_utils.py", package = "rgee")
if (!file.exists(ee_current_version)) {
stop(
sprintf(
"The file %s does not exist in your system. Please re-install rgee: %s",
ee_current_version,
"remotes::install_github(\"r-spatial/rgee\")."
)
)
}
ee_utils <- ee_connect_to_py(path = ee_current_version, n = 5)
earthengine_version <- ee_utils_py_to_r(ee_utils$ee_getversion())
# is earthengine-api greater than 0.1.317?
#if (as.numeric(gsub("\\.","",earthengine_version)) < 01317) {
# warning(
# "Update your earthengine-api installations to v0.1.317 or greater. ",
# "Earlier versions are not compatible with recent ",
# "changes to the Earth Engine backend."
# )
# }
list(earthengine_version = earthengine_version, ee_utils = ee_utils)
}
ee_check_packages <- function(fn_name, packages) {
pkg_exists <- rep(NA, length(packages))
counter <- 0
for(package in packages) {
counter <- counter + 1
pkg_exists[counter] <- requireNamespace(package, quietly = TRUE)
}
if (!all(pkg_exists)) {
to_install <- packages[!pkg_exists]
to_install_len <- length(to_install)
error_msg <- sprintf(
"%s required the %s: %s. Please install %s first.",
bold(fn_name),
if (to_install_len == 1) "package" else "packages",
paste0(bold(to_install), collapse = ", "),
if (to_install_len == 1) "it" else "them"
)
stop(error_msg)
}
}
# utils-credentials.R
ee_create_credentials_drive <- function(user=NULL, ee_utils, quiet) {
# check googledrive R package installation
ee_check_packages("ee_Initialize", "googledrive")
# Check sanity of earth-engine and return ee_utils.py module
init <- ee_check_init()
ee_utils <- init$ee_utils
# setting drive folder
if (is.null(user)) {
ee_path <- ee_utils_py_to_r(ee_utils$ee_path())
ee_path_user <- ee_path
} else {
ee_path <- ee_utils_py_to_r(ee_utils$ee_path())
ee_path_user <- sprintf("%s/%s", ee_path, user)
}
# Load GD credentials (googledrive::drive_auth)
repeat {
full_credentials <- list.files(path = ee_path_user, full.names = TRUE)
drive_condition <- grepl(".*_.*@.*", basename(full_credentials))
# If the googledrive credential does not exist run googledrive::drive_auth
if (!any(drive_condition)) {
suppressMessages(
googledrive::drive_auth(
email = NULL,
cache = ee_path_user
)
)
} else {
drive_credentials <- full_credentials[drive_condition]
email <- sub("^[^_]*_", "", basename(drive_credentials)) # Obtain the email
suppressMessages(
googledrive::drive_auth(
email = email,
cache = ee_path_user
)
)
# This lines is to avoid that users have multiple token file.
# It delete the older one if the system detect two different token files.
new_full_credentials <- list.files(path = ee_path_user, full.names = TRUE)
new_drive_condition <- grepl(".*_.*@.*", basename(new_full_credentials))
if (sum(new_drive_condition) > 1) {
files_credentials_time <- file.info(new_full_credentials[new_drive_condition])$ctime
drive_credential_to_remove <- new_full_credentials[which.min(files_credentials_time)]
if (!quiet) {
message("Removing previous Google Drive Token ....")
}
file.remove(drive_credential_to_remove)
}
break
}
}
# Move credential to the main folder is user is set
if (!is.null(user)) {
# Clean previous and copy new GD credentials in ./earthengine folder
clean_drive <- list.files(ee_path, ".*_.*@.*", full.names = TRUE) %in% list.dirs(ee_path)
unlink(
list.files(ee_path, ".*_.*@.*", full.names = TRUE)[!clean_drive]
)
file.copy(
from = drive_credentials,
to = sprintf("%s/%s", ee_path, basename(drive_credentials)),
overwrite = TRUE
)
}
invisible(drive_credentials)
}
ee_create_credentials_gcs_ <- function(user, ee_utils) {
# check packages
ee_check_packages("ee_Initialize", "googleCloudStorageR")
# Check sanity of earth-engine and return ee_utils.py module
init <- ee_check_init()
ee_utils <- init$ee_utils
# setting gcs folder
if (is.null(user)) {
ee_path <- suppressWarnings(ee_utils_py_to_r(ee_utils$ee_path()))
ee_path_user <- ee_path
} else {
ee_path <- suppressWarnings(ee_utils_py_to_r(ee_utils$ee_path()))
ee_path_user <- sprintf("%s/%s", ee_path, user)
}
# gcs_credentials
full_credentials <- list.files(path = ee_path_user, full.names = TRUE)
gcs_condition <- grepl(".json", full_credentials)
if (!any(gcs_condition)) {
gcs_text <- paste(
sprintf("Unable to find a service account key (SAK) file in: %s", crayon::bold(ee_path_user)),
"To solve this problem:",
" 1) download it from your Google cloud console",
" 2) validate it using rgee::ee_utils_sak_validate [OPTIONAL].",
" 3) Use rgee::ee_utils_sak_copy to set the SaK in rgee.",
"A tutorial to obtain the SAK file is available at:",
"> https://r-spatial.github.io/rgee/articles/rgee05.html",
crayon::bold("As long as you haven't saved a SKA file, the following functions will not work:"),
"- rgee::ee_gcs_to_local()",
"- ee_extract(..., via = \"gcs\")",
"- ee_as_raster(..., via = \"gcs\")",
"- ee_as_stars(..., via = \"gcs\")",
"- ee_as_sf(..., via = \"gcs\")",
"- sf_as_ee(..., via = \"gcs_to_asset\")",
"- gcs_to_ee_image",
"- raster_as_ee",
"- local_to_gcs",
"- stars_as_ee",
sep = "\n"
)
gcs_info <- list(path = NA, message = gcs_text)
} else {
gcs_credentials <- full_credentials[gcs_condition]
googleCloudStorageR::gcs_auth(gcs_credentials)
if (!is.null(user)) {
unlink(list.files(ee_path, ".json", full.names = TRUE))
file.copy(
from = gcs_credentials,
to = sprintf("%s/%s", ee_path, basename(gcs_credentials)),
overwrite = TRUE
)
gcs_info <- list(path = gcs_credentials, message = NA)
} else {
gcs_info <- list(path = gcs_credentials, message = NA)
}
}
gcs_info
}
# utils-auth.R & ee_Initialize.R
ee_sessioninfo <- function(email = NULL,
user = NULL,
drive = NULL,
gcs = NULL) {
init <- ee_check_init()
ee_utils <- init$ee_utils
oauth_func_path <- system.file("python/ee_utils.py", package = "rgee")
utils_py <- ee_source_python(oauth_func_path)
sessioninfo <- sprintf(
"%s/rgee_sessioninfo.txt",
ee_utils_py_to_r(utils_py$ee_path())
)
if (is.null(email)) {
email <- NA
}
# Loading all the credentials: earthengine, drive and GCS.
drive_credentials <- NA
gcs_credentials <- list(path = NA, message = NA)
if (drive) {
ee_check_packages("ee_Initialize", "googledrive")
# drive init message
# if (!quiet) ee_message_02(init = TRUE)
# If the user is not NULL copy the drive credential in the subfolder
drive_credentials <- ee_create_credentials_drive(user, ee_utils, quiet = quiet)
# test_drive_privileges(user)
# if (!quiet) ee_message_02(init = FALSE)
}
if (gcs) {
ee_check_packages("ee_Initialize", "googleCloudStorageR")
# if (!quiet) ee_message_03(init=TRUE, gcs_credentials)
# Load GCS credentials
gcs_credentials <- ee_create_credentials_gcs(user, ee_utils)
# if (!quiet) ee_message_03(init=FALSE, gcs_credentials)
}
## rgee session file
# options(rgee.gcs.auth = gcs_credentials[["path"]])
# if (!quiet) ee_message_04(init = TRUE)
df <- data.frame(
email = email, user = user, drive_cre = drive_credentials, gcs_cre = gcs_credentials[["path"]],
stringsAsFactors = FALSE
)
write.table(df, sessioninfo, row.names = FALSE)
}
source("ee_sessioninfo.R")
ee_sessioninfo(email = "ndef", user = "GEEのユーザ名", drive=T, gcs=F)
利用するライブラリは変わりません。
library(rgee)
library(rgeeExtra)
library(tidyverse)
library(leaflet)
library(sf)
library(mapview)
# rmarkdownによりmapviewの結果を出力するための変更。Rで見る場合はTRUEのままでよい。
mapviewOptions(fgb=FALSE)
表示する地点
RでGoogle Earth Engineの植生指数時系列データからBFASTによる変化検知を試すで表示した，兵庫県姫路市の夢前スマートIC付近を表示させます。
geom_point <- c(134.692408, 34.963211)
# https://ja.wikipedia.org/wiki/%E5%A4%A2%E5%89%8D%E3%83%90%E3%82%B9%E3%82%B9%E3%83%88%E3%83%83%E3%83%97
Sentinel-2の雲除去
今回はSentinel-2を背景画像として表示させるにあたって，2023年10月にGoogleがリリースしたCloud Score+データセットを利用した雲除去を行います。詳しくは，以下の解説ページを参考にしてください。
Cloud Score+データセットはcsとcs_cdfという2つのバンドがあり，どちらかを選ぶことになります。上記ページによると，csバンドは鮮明なピクセルだけが必要な場合に最適で，cs_cdfは潜在的に，より多くの「使用可能な」ピクセルを使用できるとのことです。ここでは，csバンドを用いています。
なお，eemontを読み込んでからだと読み込み途中で止まるようなので，eemontが必要な場合には，Cloud Score+データセットを先に読み込んでおく方がいいようです。
ee_geom <- ee$Geometry$Point(geom_point)
csPlus <- ee$ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED')
QA_BAND <- 'cs'
CLEAR_THRESHOLD <- 0.60
S2_SR <- ee$ImageCollection("COPERNICUS/S2_SR_HARMONIZED")$
linkCollection(csPlus, list(QA_BAND))$
filterBounds(ee_geom)$
filterDate("2024-05-01", "2024-09-01")$
map(function(img) {img$updateMask(img$select(QA_BAND)$gte(CLEAR_THRESHOLD))})$
scaleAndOffset()$
median()
S2_SR_RGB <- S2_SR$
visualize(bands=c("B4", "B3", "B2"), min=0, max=0.25)
map_S2_SR_RGB <- Map$addLayer(S2_SR_RGB)
leaflet() %>%
addTiles(
urlTemplate = map_S2_SR_RGB$rgee$tokens,
attribution = 'Contains modified Copernicus Sentinel data | Cloud Score+ S2_HARMONIZED V1',
group="S2",
options = leaflet::tileOptions(opacity = 1)
)%>%
setView(geom_point[1], geom_point[2], zoom=14) %>%
addScaleBar(position="bottomleft")
RでGoogle Earth Engineの光学衛星データを表示するで紹介したQA60バンドによる雲の除去も同様に表示してみます。こちらでは，雲や雲の影と思われる箇所がありますが，上のCloud Score+データセットでは，明らかな雲や影は見当たりません。
maskS2clouds <- function(image) {
qa <- image$select('QA60')
#Bits 10 and 11 are clouds and cirrus, respectively.
cloudBitMask <- bitwShiftL(1, 10)
cirrusBitMask <- bitwShiftL(1, 11)
#Both flags should be set to zero, indicating clear conditions.
mask = qa$bitwiseAnd(cloudBitMask)$eq(0)$
And(qa$bitwiseAnd(cirrusBitMask)$eq(0))
return(image$updateMask(mask)$divide(10000))
}
S2_SR_QA <- ee$ImageCollection("COPERNICUS/S2_SR_HARMONIZED")$
filterDate("2024-05-01", "2024-09-01")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUDY_PIXEL_PERCENTAGE',20))$
map(maskS2clouds)$median()
S2_SR_QA_RGB <- S2_SR_QA$
visualize(bands=c("B4", "B3", "B2"), min=0, max=0.25)
map_S2_SR_QA_RGB <- Map$addLayer(S2_SR_QA_RGB)
leaflet() %>%
addTiles(
urlTemplate = map_S2_SR_QA_RGB$rgee$tokens,
attribution = 'Contains modified Copernicus Sentinel data',
group="S2",
options = leaflet::tileOptions(opacity = 1)
)%>%
setView(geom_point[1], geom_point[2], zoom=14) %>%
addScaleBar(position="bottomleft")
参考として，同一地点・期間のLandsat衛星の画像も表示しておきます。
cloudcover <- 50
# Landsat8
# https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2
L8_SR <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUD_COVER', cloudcover))$
preprocess()
# Landsat9
# https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC09_C02_T1_L2
L9_SR <- ee$ImageCollection("LANDSAT/LC09/C02/T1_L2")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUD_COVER', cloudcover))$
preprocess()
## 各衛星データのマージ
Landsat_SR <- L8_SR$merge(L9_SR)$
filterDate("2024-05-01", "2024-09-01")$
median()
Landsat_SR_RGB <- Landsat_SR$
visualize(bands=c("SR_B4", "SR_B3", "SR_B2"), min=0, max=0.25)
map_Landsat_SR_RGB <- Map$addLayer(Landsat_SR_RGB)
leaflet() %>%
addTiles(
urlTemplate = map_Landsat_SR_RGB$rgee$tokens,
attribution = 'Landsat image courtesy of the U.S. Geological Survey',
group="Landsat",
options = leaflet::tileOptions(opacity = 1)
) %>%
setView(geom_point[1], geom_point[2], zoom=14) %>%
addScaleBar(position="bottomleft")
Global forest changeデータセットを表示してみる
Global Forest Watchで公開されているHansen et al.(2003)の森林被覆・変化データがGoogleEarthEngineで利用できます。このデータセットはLandsat衛星画像を機械学習により分析したもので，現状ではオリジナル（バージョン 1.0）の2000～2012年のデータと，Landsat8を加えるなどの更新を行ったバージョン1.11の2011 ～2023年のデータをあわせたデータが公開されています。両者では利用している衛星データやアルゴリズムが異なり，統合してデータを利用する際には注意が必要です。
Introduction to Forest Change Analysis in Earth Engine
Hansen Global Forest Change v1.11 (2000-2023)
データセットは以下で読み込みます。
gfc <- ee$Image('UMD/hansen/global_forest_change_2023_v1_11')
緑が2000年時点から森林，赤が2022年までに森林が失われた（loss）したところ，青は2012年までに新たに森林になった（gain）ところ，マゼンダはlossとgainが両方おきたところを表示しています。ここでいうlossは森林（5m以上の植生があるところ）でおきた「stand-replacement disturbance」と定義されています。なお，gainはver1.11データセットでは更新がありませんので，lossとgainで対象期間が異なることに注意してください。
gfcVisParam <- list(
bands= c('loss', 'treecover2000', 'gain'),
max=c(1, 255, 1)
)
Map$setCenter(lon=geom_point[1], lat=geom_point[2], zoom=14)
Map$addLayer(gfc, gfcVisParam, 'forest cover, loss, gain')
つづいて，lossが観測された年を示してみます。赤いほど，最近のlossを示しています。
treeLossVisParam <- list(
bands= 'lossyear',
min= 0,
max= 23,
palette=c('yellow', 'red')
)
Map$addLayer(gfc, treeLossVisParam, 'tree loss year')
上記のSentinel-2を背景画像として，lossが観測された年を表示してみます。
lossyear <- gfc$
visualize(bands='lossyear', min=0, max=23, palette=c('yellow', 'red'))
map_lossyear <- Map$addLayer(lossyear) +
Map$addLegend(
list(min=0, max=23, palette=c('yellow', 'red')),
position="bottomright", color_mapping = "discrete")
leaflet() %>%
addTiles(
urlTemplate = map_S2_SR_RGB$rgee$tokens,
attribution = 'Contains modified Copernicus Sentinel data | Cloud Score+ S2_HARMONIZED V1',
group="S2",
options = leaflet::tileOptions(opacity = 1)
)%>%
addTiles(
urlTemplate = map_lossyear$rgee$tokens,
attribution = "| <a href='https://doi.org/10.1126/science.1244693' target='_blank'>Hansen Global Forest Change</a>",
group="lossyear",
layerId = "lossyear",
options = leaflet::tileOptions(opacity = 0.7)
) %>%
setView(geom_point[1], geom_point[2], zoom=14) %>%
addScaleBar(position="bottomleft") %>%
addLayersControl(
baseGroups = c("S2"),
overlayGroups = c("lossyear"),
position = "topleft")
Global forest changeデータセットを都道府県ごとに集計してみる
都道府県のポリゴンを読み込む
都道府県単位で森林損失面積などの集計を行うこともできます。
GooogleEarthEngineでは行政区域のデータとして，国連食糧農業機関（FAO）が開発した「Global Administrative Unit Layers（GAUL）」を利用できます。なお，ライセンス上，非商用目的に利用が制約されているので注意してください。データセットは，国，第1レベル（県など），第2レベル（市区町村など）に分かれており，今回は県のデータが必要なため，第1レベルのデータセットを読み込みます。
japan_prefectures <- ee$FeatureCollection("FAO/GAUL/2015/level1")$
filter(ee$Filter$eq('ADM0_NAME', 'Japan'))
#‘ 都道府県名が若干独特なので，以下により確認してください（結果は省略）。
adm1_names <- japan_prefectures$aggregate_array('ADM1_NAME')$getInfo()
print(adm1_names)
兵庫県を対象とするため，以下により抽出したデータを作っておきます。一応，確認のために表示してみます。
pref <- japan_prefectures$filter(ee$Filter$eq('ADM1_NAME', 'Hyoogo'))
Map$centerObject(eeObject = pref)
Map$addLayer(pref)%>%
addTiles(urlTemplate = "", attribution = '| Source of Administrative boundaries: The Global Administrative Unit Layers (GAUL) dataset,
implemented by FAO within the CountrySTAT and Agricultural Market Information System
(AMIS) projects')
なお，海岸線の位置など若干おかしなところがあるので，正確を期すのであれば，国土地理院などのデータを利用したほうが確実です。例えば，以前紹介した地球地図日本のデータをGEE上にアップロードして利用することもできます。
gm_jpn <- read_sf(gm_jpn.path) %>%
group_by(nam) %>% summarize()
pref <- gm_jpn %>% filter(nam == "Hyogo Ken") %>% sf_as_ee()
年ごとの森林の損失面積をグラフ表示してみる
lossyearデータを利用して，兵庫県内の年ごとの森林損失面積をグラフとして表示してみましょう。なお，Global Forest Watchには国・地域ごとのダッシュボードもあり，兵庫県のデータもみることができます。
lossImage <- gfc$select('loss')
lossAreaImage <- lossImage$multiply(ee$Image$pixelArea())
lossYear <- gfc$select('lossyear')
lossByYear <- lossAreaImage$addBands(lossYear)$reduceRegion(
reducer= ee$Reducer$sum()$group(
groupField= 1
),
geometry= pref,
scale= 30,
maxPixels= 1e9
)
# print(lossByYear$getInfo())
statsFormatted <- ee$List(lossByYear$get('groups'))$map(
ee_utils_pyfunc(
function(el) {
d <- ee$Dictionary(el)
losstable <- list(ee$Number(d$get('group'))$format("20%02d"),
d$get('sum'))
return(losstable)
}
)
)
statsDictionary <- ee$Dictionary(statsFormatted$flatten())
# print(statsDictionary$getInfo())
statsData <- statsDictionary$getInfo() %>%
enframe(name="year", value="value") %>%
mutate(value = as.numeric(unlist(value)))
ggplot(statsData,aes(x=year,y=value/10000)) + geom_bar(stat='identity') +
scale_y_continuous(name = "Area [ha]", expand = expansion(mult = c(0, .1))) +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme_classic()
都道府県ごとの森林の損失と増加を表示してみる
続いて，各都道府県の森林の損失と増加についてグラフ表示してみましょう。
# 森林の損失と増加を計算
treecoverImage <- gfc$select("treecover2000")
lossImage <- gfc$select("loss")
gainImage <- gfc$select("gain")
# 都道府県ごとの損失と増加の面積を計算
treecoverArea <- treecoverImage$reduceRegions(
collection = japan_prefectures,
reducer = ee$Reducer$sum(),
scale = 30
)
lossArea <- lossImage$reduceRegions(
collection = japan_prefectures,
reducer = ee$Reducer$sum(),
scale = 30
)
gainArea <- gainImage$reduceRegions(
collection = japan_prefectures,
reducer = ee$Reducer$sum(),
scale = 30
)
# 結果をデータフレームに変換
treecover_df <- ee_as_sf(treecoverArea) %>% sf::st_drop_geometry()
loss_df <- ee_as_sf(lossArea) %>% sf::st_drop_geometry()
gain_df <- ee_as_sf(gainArea) %>% sf::st_drop_geometry()
# データフレームをマージ
forest_change_df <- treecover_df %>%
left_join(loss_df, by = "ADM1_NAME") %>%
left_join(gain_df, by = "ADM1_NAME") %>%
select(c("ADM1_NAME", "sum.x", "sum.y", "sum"))
colnames(forest_change_df) <- c("Prefecture", "treecover", "Loss", "Gain")
森林の損失と増加の面積を棒グラフにプロット
ggplot(forest_change_df, aes(x = Prefecture)) +
geom_bar(aes(y = -Loss/10000), stat = "identity", fill = "red", alpha = 0.7) +
geom_bar(aes(y = Gain/10000), stat = "identity", fill = "green", alpha = 0.7) +
labs(title = "Forest Loss and Gain by Prefecture in Japan",
x = "Prefecture",
y = "Area [ha]") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
森林の損失と増加の面積割合を棒グラフにプロット
ggplot(forest_change_df, aes(x = Prefecture)) +
geom_bar(aes(y = -Loss/treecover * 100), stat = "identity", fill = "red", alpha = 0.7) +
geom_bar(aes(y = Gain/treecover * 100), stat = "identity", fill = "green", alpha = 0.7) +
labs(title = "Forest Loss and Gain by Prefecture in Japan",
x = "Prefecture",
y = "Area [proportion to treecover2000 (%)]") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
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
## [3] LC_MONETARY=Japanese_Japan.utf8 LC_NUMERIC=C
## [5] 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] mapview_2.11.2 sf_1.0-17 leaflet_2.2.2 lubridate_1.9.3
## [5] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
## [9] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [13] tidyverse_2.0.0 rgeeExtra_0.0.1 rgee_1.1.7 here_1.0.1
## [17] fs_1.6.4 rmarkdown_2.28
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.16.0 jsonlite_1.8.8
## [3] magrittr_2.0.3 farver_2.1.2
## [5] vctrs_0.6.5 memoise_2.0.1
## [7] askpass_1.2.1 base64enc_0.1-3
## [9] terra_1.7-78 webshot_0.5.5
## [11] htmltools_0.5.8.1 usethis_3.0.0
## [13] curl_5.2.1 raster_3.6-26
## [15] googleCloudStorageR_0.7.0 KernSmooth_2.23-22
## [17] htmlwidgets_1.6.4 plyr_1.8.9
## [19] cachem_1.0.8 mime_0.12
## [21] lifecycle_1.0.4 pkgconfig_2.0.3
## [23] webshot2_0.1.1 Matrix_1.6-5
## [25] R6_2.5.1 fastmap_1.1.1
## [27] shiny_1.9.1 digest_0.6.34
## [29] colorspace_2.1-0 ps_1.8.0
## [31] rprojroot_2.0.4 leafem_0.2.3
## [33] pkgload_1.4.0 crosstalk_1.2.1
## [35] geojson_0.3.5 googleAuthR_2.0.2
## [37] labeling_0.4.3 fansi_1.0.6
## [39] timechange_0.3.0 httr_1.4.7
## [41] compiler_4.3.3 gargle_1.5.2
## [43] proxy_0.4-27 remotes_2.5.0
## [45] bit64_4.0.5 withr_3.0.1
## [47] DBI_1.2.3 highr_0.11
## [49] pkgbuild_1.4.4 protolite_2.3.1
## [51] openssl_2.2.2 sessioninfo_1.2.2
## [53] classInt_0.4-10 tools_4.3.3
## [55] chromote_0.3.1 units_0.8-5
## [57] googledrive_2.1.1 odbc_1.4.2
## [59] zip_2.3.1 httpuv_1.6.15
## [61] glue_1.7.0 satellite_1.0.5
## [63] promises_1.3.0 grid_4.3.3
## [65] geojsonio_0.11.3 generics_0.1.3
## [67] gtable_0.3.5 leaflet.providers_2.0.0
## [69] tzdb_0.4.0 class_7.3-22
## [71] websocket_1.4.2 hms_1.1.3
## [73] sp_2.1-4 utf8_1.2.4
## [75] pillar_1.9.0 later_1.3.2
## [77] lattice_0.22-5 bit_4.0.5
## [79] tidyselect_1.2.1 miniUI_0.1.1.1
## [81] knitr_1.48 V8_5.0.1
## [83] crul_1.5.0 stats4_4.3.3
## [85] xfun_0.48 devtools_2.4.5
## [87] DT_0.33 stringi_1.8.4
## [89] geojsonsf_2.0.3 lazyeval_0.2.2
## [91] ggmap_4.0.0 yaml_2.3.10
## [93] httpcode_0.3.0 evaluate_1.0.0
## [95] codetools_0.2-19 cli_3.6.1
## [97] xtable_1.8-4 reticulate_1.39.0
## [99] munsell_0.5.1 processx_3.8.4
## [101] jquerylib_0.1.4 jqr_1.3.5
## [103] Rcpp_1.0.12 png_0.1-8
## [105] ellipsis_0.3.2 assertthat_0.2.1
## [107] blob_1.2.4 jpeg_0.1-10
## [109] profvis_0.4.0 urlchecker_1.0.1
## [111] bitops_1.0-9 scales_1.3.0
## [113] e1071_1.7-14 crayon_1.5.3
## [115] rlang_1.1.3
Discussion