RでGoogle Earth Engineの光学衛星データを表示する
変更履歴・概要
2024/3/2 改訂
実行環境をR4.3.3に更新しています。画像中に使用した衛星データの引用ができるように見直しました。leafsyncはrmarkdownによる出力ができていないので，ローカルのhtmlへの保存のみ追記しました。（出力方法が分かりましたら，コメント等頂ければありがたいです）。
2024/2/4 初版
画像の出力ができていない箇所，使用した衛星データの引用が画像中に出力できていない箇所がいくつかあります（出力方法が分かりましたら，コメント等頂ければありがたいです）。
できるようになること
RからGoogle Earth Engine（GEE）を操作して，光学衛星データを表示させます。GEEでは10～数十mレベルの中解像度光学センサを持つ，LandsatとSentinel-2のデータを利用できます。
初期設定
RでGEEを操作する設定は，RでGoogle Earth Engineを操作できるようにするを参考にしてください。
library(rgee)
library(rgeeExtra)
ee_Initialize(project = 'プロジェクト名')
衛星データの選択
衛星ごとに見てみましょう。なお，衛星データは前処理のレベルごとにデータが公開されている場合がありますので，必要な処理がされているデータを用いることになります。前処理については，宙畑のページに分かりやすく整理してあるので参考にしてください。ここでは，大気補正済データ（Level2）のasset IDを示します。
取得する衛星データ(1) Landsat
Landsatはアメリカが打ち上げた地球観測衛星で，最初のLandsat-1は1972年に打ち上げられました。Landsat-4以降は30mの解像度で，16日に1回地表を撮影しています。2024年現在で運用中なのはLandsat-7（ただしScan Line Correctorの故障によりSLC-offのデータ（What is Landsat 7 ETM+ SLC-off data?）），Landsat-8，Landsat-9の3機です（Landsat Missions Timeline）。
GEEで取得可能なデータや特性等は以下のページで確認してください。
https://developers.google.com/earth-engine/datasets/catalog/landsat
2024年時点ではLandsatデータはCollection2の利用が推奨されています。データの種類については@NagasakiNoHitoさんの説明が分かりやすいです。
Landsat-8の大気補正済データ（Tier1）を以下により取得できます。
L8_SR <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")
取得する衛星データ(2) Sentinel-2
Sentinel-2は欧州（ESA）が打ち上げた地球観測衛星で，10mの解像度で10日に1回地表を撮影しています。2024年現在では2機体制で運用されています。GEEではLevel 1Cデータが2015年6月から，Level 2Aデータが2017年3月から利用可能です。GEEで取得可能なデータや特性等は以下のページで確認してください。
https://developers.google.com/earth-engine/datasets/catalog/sentinel-2
S2_SR <- ee$ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
衛星データの取得～表示：Landsat-8を例に
衛星データをみるためには，衛星データの指定→範囲等によるフィルタリング→（前処理）→（減算）→（表示パラメータの設定をしてから）表示という手順を踏みます。順に書くと以下のとおりになります。
とりあえず，2022年の兵庫県庁周辺のLandsat-8データを取得して，表示させてみましょう。
衛星データの指定：Landsat-8の大気補正済データ（Tier1）
L8_SR <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")
データ取得日の範囲の指定。例では2022/1/1～2023/1/1の1年間を指定しています。
L8_SR <- L8_SR$filterDate("2022-01-01", "2023-01-01")
データ取得範囲の指定（この場合はgeom_pointが重なるデータ）
geom_point <- c(135.1830778, 34.69125701) # 兵庫県庁
ee_geom <- ee$Geometry$Point(geom_point)
L8_SR <- L8_SR$filterBounds(ee_geom)
雲が多いデータの除去
L8_SR <- L8_SR$filter(ee$Filter$lt('CLOUD_COVER', 50))
ここまでで取得する衛星データの指定は終わりです。取得したデータは以下により確認できます。rgeeライブラリのee_get_date_ic()を用いても確認できます。
L8_SR$aggregate_array('system:id')$getInfo()
## [1] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220102" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220408"
## [3] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220729" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220915"
## [5] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221001" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221102"
## [7] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221118" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221220"
前処理。ここでは，rgeeExtraライブラリのpreprocess()で処理しています。
L8_SR <- L8_SR$preprocess()
減算する（1つの画像にまとめる）ため，取得された各ピクセルの中央値を計算。なお，L8_SR$reduce(ee$Reducer$median())でも中央値が出力されますが，出力バンド名に_medianが追加されます。減算については，GEE公式ページの記載を参考にしてください。
L8_SR_median <- L8_SR$median()
表示させるバンド等の設定
L8_SR_visParams <- list(
bands = c("SR_B4", "SR_B3", "SR_B2"),
min = 0.0, max=0.3, gamma=1.2)
表示
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 12)
Map$addLayer(
L8_SR_median,
L8_SR_visParams,
"true-color: Red-Green-Blue"
) %>%
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_median <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
filterDate("2022-01-01", "2023-01-01")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUD_COVER', 50))$
preprocess()$median()
Landsat-8データの前処理について
上記の例ではrgeeExtraライブラリのpreprocess()を用いました。処理の内容はscale factorの補正と品質評価（QA）バンドによる雲と影のピクセルの除去です。Digital Earth AFRICA のページも分かりやすいと思います。
ここでは，コードを書いて処理してみたいと思います。なお，インターネット上のコードの多くはCollection 1データに対する前処理となっており，Collection 2データのQAバンド名やBitとは異なっているので，注意してください。
preprocessL8C2 <- function(image) {
cloudsBitMask <- bitwShiftL(1, 3)
cloudShadowBitMask <- bitwShiftL(1, 4)
cirrusBitMask <- bitwShiftL(1, 2)
qa <- image$select('QA_PIXEL')
mask <- qa$bitwiseAnd(cloudShadowBitMask)$eq(0)$
And(qa$bitwiseAnd(cloudsBitMask)$eq(0))$
And(qa$bitwiseAnd(cirrusBitMask)$eq(0))
opticalBands <- image$select('SR_B.')$multiply(0.0000275)$add(-0.2)
thermalBands <- image$select('ST_B.*')$multiply(0.00341802)$add(149.0)
return(image$addBands(opticalBands, NULL, TRUE)$
addBands(thermalBands, NULL, TRUE)$
updateMask(mask))
}
L8_SR_median <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
filterDate("2022-01-01", "2023-01-01")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUD_COVER', 50))$
map(preprocessL8C2)$median()
衛星データの指定例
rgeeExtraのee_ImageCollection_closest()を使うことで，以下の例のように指定日の直近のデータを指定することもできます。
L8_SR_1 <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
filter(ee$Filter$eq("WRS_PATH", 110))$
filter(ee$Filter$eq("WRS_ROW", 36)) %>%
ee_ImageCollection_closest("2019-03-15", 1, "Month") %>%
ee$ImageCollection$first() %>%
ee$Image$preprocess()
特定の期間で雲が最も少ないデータを選ぶこともできます。
L8_SR_2 <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
filter(ee$Filter$eq("WRS_PATH", 110))$
filter(ee$Filter$eq("WRS_ROW", 36))$
filterDate("2020-01-01", "2020-04-01")$
sort("CLOUD_COVER")$first() %>%
ee$Image$preprocess()
衛星データの取得～表示：Sentinel-2の場合
Sentinel-2も同じように表示できます。
Sentinel-2は雲の除去がQA60バンドを利用する方法と，COPERNICUS/S2_CLOUD_PROBABILITYを利用する方法があります。ここでは，QA60バンドを利用して除去してみます。
https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR_HARMONIZED
なお，rgeeExtraライブラリのpreprocess()はSentinel-2データでも処理できるようになっていますが，HARMONIZED collectionに対しては雲の除去がうまく処理できませんでした。
データの読み込み・前処理
geom_point <- c(135.1830778, 34.69125701) # 兵庫県庁
ee_geom <- ee$Geometry$Point(geom_point)
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 <- ee$ImageCollection("COPERNICUS/S2_SR_HARMONIZED")$
filterDate("2022-01-01", "2023-01-01")$
filterBounds(ee_geom)$
filter(ee$Filter$lt('CLOUDY_PIXEL_PERCENTAGE',20))
S2_SR_median <- S2_SR$map(maskS2clouds)$median()
表示の設定から表示まで
S2_SR_visParams <- list(
bands = c("B4", "B3", "B2"),
min = 0.0, max=0.3)
Map$setCenter(lon = geom_point[1], lat = geom_point[2], zoom = 12)
Map$addLayer(
S2_SR_median,
S2_SR_visParams,
"true-color composite: Red-Green-Blue"
) %>%
addTiles(urlTemplate = "", attribution = '| Contains modified Copernicus Sentinel data (2022)')
メタデータ等の表示
取得データのメタデータや値の表示方法をいくつか紹介します。
取得しているデータのメタデータを表示
ee_print(L8_SR)
## Registered S3 method overwritten by 'geojsonsf':
## method from
## print.geojson geojson
## ──────────────────────────────────────────────────────────────────────────── Earth Engine ImageCollection ──
## ImageCollection Metadata:
## - Class : ee$ImageCollection
## - Number of Images : 8
## - Number of Properties : 27
## - Number of Pixels* : 8887057872
## - Approximate size* : 132.81 GB
## Image Metadata (img_index = 0):
## - ID : LANDSAT/LC08/C02/T1_L2/LC08_110036_20220102
## - system:time_start : 2022-01-02 01:35:00.952
## - system:time_end : 2022-01-02 01:35:00.952
## - Number of Bands : 19
## - Bands names : SR_B1 SR_B2 SR_B3 SR_B4 SR_B5 SR_B6 SR_B7 SR_QA_AEROSOL ST_B10 ST_ATRAN ST_CDIST ST_DRAD ST_EMIS ST_EMSD ST_QA ST_TRAD ST_URAD QA_PIXEL QA_RADSAT
## - Number of Properties : 94
## - Number of Pixels* : 1110882234
## - Approximate size* : 16.60 GB
## Band Metadata (img_band = 'SR_B1'):
## - EPSG (SRID) : WGS 84 / UTM zone 53N (EPSG:32653)
## - proj4string : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs
## - Geotransform : 30 0 393285 0 -30 3946815
## - Nominal scale (meters) : 30
## - Dimensions : 7518 7777
## - Number of Pixels : 58467486
## - Data type : DOUBLE
## - Approximate size : 894.73 MB
## ────────────────────────────────────────────────────────────────────────────────────────────────────────────
## NOTE: (*) Properties calculated considering a constant geotransform and data type.
データ内のイメージ数
L8_SR$size()$getInfo()
## [1] 8
取得したイメージの一覧
L8_SR$aggregate_array('system:id')$getInfo()
## [1] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220102" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220408"
## [3] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220729" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20220915"
## [5] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221001" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221102"
## [7] "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221118" "LANDSAT/LC08/C02/T1_L2/LC08_110036_20221220"
データ取得日の表示
eedate_to_rdate(L8_SR$first()$get("system:time_start"))
## [1] "2022-01-02 01:35:00 GMT"
データ取得日の範囲の表示
range <- L8_SR$reduceColumns(
ee$Reducer$minMax(),
list("system:time_start")
)
col_min <- eedate_to_rdate(range$get("min"))
col_max <- eedate_to_rdate(range$get("max"))
cat("Date range: ", as.character(col_min), as.character(col_max))
## Date range: 2022-01-02 01:35:00.952 2022-12-20 01:35:11.125
イメージのプロパティ名の表示。L8_SR$first()はrgeeExtraライブラリを利用していればL8_SR[[1]]としてもよい。
L8_SR$first()$propertyNames()$getInfo()
## [1] "DATUM" "REFLECTANCE_ADD_BAND_3"
## [3] "REFLECTANCE_ADD_BAND_4" "REFLECTANCE_ADD_BAND_5"
## [5] "REFLECTANCE_ADD_BAND_6" "REFLECTANCE_ADD_BAND_7"
## [7] "system:footprint" "REFLECTIVE_SAMPLES"
## [9] "system:version" "GROUND_CONTROL_POINTS_VERSION"
## [11] "SUN_AZIMUTH" "DATA_SOURCE_TIRS_STRAY_LIGHT_CORRECTION"
## [13] "UTM_ZONE" "DATE_ACQUIRED"
## [15] "ELLIPSOID" "system:time_end"
## [17] "DATA_SOURCE_PRESSURE" "LANDSAT_PRODUCT_ID"
## [19] "STATION_ID" "TEMPERATURE_ADD_BAND_ST_B10"
## [21] "DATA_SOURCE_REANALYSIS" "REFLECTANCE_MULT_BAND_7"
## [23] "system:time_start" "REFLECTANCE_MULT_BAND_6"
## [25] "L1_PROCESSING_LEVEL" "PROCESSING_SOFTWARE_VERSION"
## [27] "L1_DATE_PRODUCT_GENERATED" "ORIENTATION"
## [29] "REFLECTANCE_MULT_BAND_1" "WRS_ROW"
## [31] "REFLECTANCE_MULT_BAND_3" "REFLECTANCE_MULT_BAND_2"
## [33] "TARGET_WRS_ROW" "REFLECTANCE_MULT_BAND_5"
## [35] "REFLECTANCE_MULT_BAND_4" "THERMAL_LINES"
## [37] "TIRS_SSM_POSITION_STATUS" "GRID_CELL_SIZE_THERMAL"
## [39] "IMAGE_QUALITY_TIRS" "TRUNCATION_OLI"
## [41] "NADIR_OFFNADIR" "CLOUD_COVER"
## [43] "REQUEST_ID" "EARTH_SUN_DISTANCE"
## [45] "GEOMETRIC_RMSE_VERIFY" "TIRS_SSM_MODEL"
## [47] "COLLECTION_CATEGORY" "SCENE_CENTER_TIME"
## [49] "GRID_CELL_SIZE_REFLECTIVE" "SUN_ELEVATION"
## [51] "ALGORITHM_SOURCE_SURFACE_TEMPERATURE" "TEMPERATURE_MAXIMUM_BAND_ST_B10"
## [53] "CLOUD_COVER_LAND" "GEOMETRIC_RMSE_MODEL"
## [55] "ROLL_ANGLE" "COLLECTION_NUMBER"
## [57] "DATE_PRODUCT_GENERATED" "L1_REQUEST_ID"
## [59] "DATA_SOURCE_OZONE" "SATURATION_BAND_1"
## [61] "DATA_SOURCE_WATER_VAPOR" "SATURATION_BAND_2"
## [63] "SATURATION_BAND_3" "IMAGE_QUALITY_OLI"
## [65] "SATURATION_BAND_4" "LANDSAT_SCENE_ID"
## [67] "SATURATION_BAND_5" "MAP_PROJECTION"
## [69] "SATURATION_BAND_6" "SENSOR_ID"
## [71] "SATURATION_BAND_7" "SATURATION_BAND_8"
## [73] "WRS_PATH" "SATURATION_BAND_9"
## [75] "TARGET_WRS_PATH" "L1_PROCESSING_SOFTWARE_VERSION"
## [77] "TEMPERATURE_MULT_BAND_ST_B10" "L1_LANDSAT_PRODUCT_ID"
## [79] "PROCESSING_LEVEL" "ALGORITHM_SOURCE_SURFACE_REFLECTANCE"
## [81] "GROUND_CONTROL_POINTS_MODEL" "SPACECRAFT_ID"
## [83] "TEMPERATURE_MINIMUM_BAND_ST_B10" "GEOMETRIC_RMSE_MODEL_Y"
## [85] "REFLECTIVE_LINES" "GEOMETRIC_RMSE_MODEL_X"
## [87] "THERMAL_SAMPLES" "system:asset_size"
## [89] "DATA_SOURCE_AIR_TEMPERATURE" "GROUND_CONTROL_POINTS_VERIFY"
## [91] "system:index" "REFLECTANCE_ADD_BAND_1"
## [93] "REFLECTANCE_ADD_BAND_2" "DATA_SOURCE_ELEVATION"
## [95] "WRS_TYPE" "system:id"
## [97] "system:bands" "system:band_names"
crsの表示
L8_SR$first()$projection()$getInfo()
## $type
## [1] "Projection"
##
## $crs
## [1] "EPSG:32653"
##
## $transform
## [1] 30 0 393285 0 -30 3946815
解像度（nominal scale）の表示
L8_SR$first()$projection()$nominalScale()$getInfo()
## [1] 30
バンド名の表示
L8_SR$first()$bandNames()$getInfo()
## [1] "SR_B1" "SR_B2" "SR_B3" "SR_B4" "SR_B5" "SR_B6"
## [7] "SR_B7" "SR_QA_AEROSOL" "ST_B10" "ST_ATRAN" "ST_CDIST" "ST_DRAD"
## [13] "ST_EMIS" "ST_EMSD" "ST_QA" "ST_TRAD" "ST_URAD" "QA_PIXEL"
## [19] "QA_RADSAT"
バンドの値の表示例
L8_SR$first()$reduceRegion(
reducer=ee$Reducer$max(),
geometry=ee_geom$buffer(distance=100),
scale=30
)$get('SR_B1')$getInfo()
## [1] 0.1006575
Landsatデータのpath/row
L8_SR$first()$get('WRS_PATH')$getInfo() #110
## [1] 110
L8_SR$first()$get('WRS_ROW')$getInfo() #36
## [1] 36
衛星画像のカラー合成
光学センサーの衛星データは，波長域の異なる反射光の強さを測定したものになります。カラー合成は，衛星データの各波長域（バンド）のうち３つを，人が目で見ることのできる赤，緑，青でそれぞれ着色して合成することで，カラー画像を作成する作業です。着色するバンドの選び方でいくつかのパターンがあります。詳しくは，以下のページを参考にしてください。
https://ja.wikipedia.org/wiki/衛星画像
http://rs.aoyaman.com/img_pro/b1.html
https://earthobservatory.nasa.gov/features/FalseColor/page6.php
https://gisgeography.com/sentinel-2-bands-combinations/
トゥルーカラー画像：赤-緑-青
人の目の配色と同じように表現した画像
L8_SR_trueParams <- list(
bands = c("SR_B4", "SR_B3", "SR_B2"),
min = 0.0, max=0.3)
L8_SR_Mtrue <- Map$addLayer(
L8_SR_median,
L8_SR_trueParams,
"true-color: Red-Green-Blue"
)
L8_SR_Mtrue %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
ナチュラルカラー画像：赤-近赤外-緑
近赤外線帯域のデータを使い、森林（植物）を緑色で強調した画像
L8_SR_natParams <- list(
bands = c("SR_B4", "SR_B5", "SR_B3"),
min = 0.0, max=0.3)
L8_SR_Mnat <- Map$addLayer(
L8_SR_median,
L8_SR_natParams,
"natural-color: Red-NIR-Green"
)
L8_SR_Mnat %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
フォルスカラー画像：近赤外-赤-緑
近赤外線帯域のデータを使い、森林（植物）を赤色で表したカラー画像です。
L8_SR_falseParams <- list(
bands = c("SR_B5", "SR_B4", "SR_B3"),
min = 0.0, max=0.3)
L8_SR_Mfalse <- Map$addLayer(
L8_SR_median,
L8_SR_falseParams,
"false-color: NIR-Red-Green"
)
L8_SR_Mfalse %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
フォルスカラー画像：短波長赤外1-近赤外-緑
短波赤外線，近赤外線帯域のデータを使い、水（黒色）や森林（植生：緑色）の分布状況が分かるカラー画像です。
L8_SR_false2Params <- list(
bands = c("SR_B6", "SR_B5", "SR_B3"),
min = 0.0, max=0.3)
L8_SR_Mfalse2 <- Map$addLayer(
L8_SR_median,
L8_SR_false2Params,
"false-color: SWIR-NIR-Green"
)
L8_SR_Mfalse2 %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
Rによる表示あれこれ
rgeeライブラリによる便利な表示をいくつか紹介します。
データの比較
各データの重ね合わせ
+でデータを重ね合わせることができます。ここでは，トゥルーカラー画像（赤-緑-青）にフォルスカラー画像（短波長赤外1-近赤外-緑）を重ねています。図の左上のボックスのようなボタンで各データの表示・非表示を選択できます。
L8_SR_false3Params <- list(
bands = c("SR_B6", "SR_B5", "SR_B3"),
min = 0.0, max=0.3, opacity=0.5)
L8_SR_Mfalse3 <- Map$addLayer(
L8_SR_median,
L8_SR_false3Params,
"false-color: SWIR-NIR-Green"
)
m <- L8_SR_Mtrue + L8_SR_Mfalse3
m %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
2つのデータの比較（Side by side view）
静止画での出力では動きませんが，中央付近のスライダーを動かすことができます。なお，leaflet.extras2ライブラリが必要なので，インストールしてください。
m <- L8_SR_Mfalse | L8_SR_Mfalse2
m %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
2つ以上のデータの比較
leafsyncライブラリのsync()により，各パネルのデータを同期させて表示できます。
library(leafsync)
m <- leafsync::sync(L8_SR_Mtrue,L8_SR_Mnat,L8_SR_Mfalse,L8_SR_Mfalse2)
m
rmarkdownによる出力では静止画像を出力できなかったので図をお示しできませんが，試してみてください。ローカルのhtmlファイルには以下で保存できます。
htmltools::save_html(m, "m.html")
ポリゴンや点などの重ね合わせ
ポリゴンや点などを重ねることもできます。
GEEにデータを送って表示
まずはGEEにデータを送って表示させます。
geom_rect <- list(c(135.17, 34.72), c(135.17, 34.68),
c(135.22, 34.68), c(135.22, 34.72), c(135.17, 34.72))
ee_geom_rect <- ee$Geometry$Polygon(geom_rect)
m <- L8_SR_Mtrue +
Map$addLayer(ee_geom_rect, list(color="FF0000"), 'bbox')
m %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
GEEにデータを送らずにR上で表示させる方法
mapviewライブラリを用いることで，ポリゴンデータや点を重ねます（Earth Engine with Rの該当ページ[1]）。
library(sf)
library(leaflet)
library(mapview)
# rmarkdownによりmapviewの結果を出力するための変更。Rで見る場合はTRUEのままでよい。
## https://github.com/r-spatial/mapview/issues/351
## https://www.r4wrds.com/intro/m_intro_rmarkdown
mapviewOptions(fgb=FALSE)
データを作成。
KobeAirport <- st_sfc(st_point(c(135.223889, 34.632778)), crs =4326)
bbox_temp <- st_sfc(st_polygon(list(rbind(c(135.17, 34.72), c(135.17, 34.68),
c(135.22, 34.68), c(135.22, 34.72), c(135.17, 34.72)))),
crs=4326)
mapview(x = bbox_temp, map = L8_SR_Mtrue %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey'), color="white", lwd=4, alpha.regions=0)
点データの重ね合わせ。ただし，mapviewではzoomの設定が引き継がれないようです。
mapview(x = KobeAirport, map = L8_SR_Mtrue %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey'))
インタラクティブマップであれば，表示後に変えられるのですが，静止画での出力では不便なので，leafletOptions()でZoomの設定を変更しておきます。なお，SetViewの設定は引き継がれないようです。
# https://stackoverflow.com/questions/71321030/locking-the-zoom-in-on-mapview
map1 <- leaflet(options=leafletOptions(minZoom = 12, maxZoom = 12)) %>%
addTiles() %>%
setView(geom_point[1], geom_point[2], 12) %>%
addTiles(
urlTemplate = L8_SR_Mtrue$rgee$tokens,
layerId = "true-color: Red-Green-Blue",
options = leaflet::tileOptions(opacity = 1)
) %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
mapview(x = KobeAirport, map = map1)
ちなみに，leafletライブラリで作成することもできます。この場合SetViewの設定は引き継がれるようです。
leaflet(KobeAirport) %>%
addTiles() %>%
addCircleMarkers(col="red") %>%
setView(geom_point[1], geom_point[2], 12) %>%
addTiles(
urlTemplate = L8_SR_Mtrue$rgee$tokens,
layerId = "true-color: Red-Green-Blue",
options = leaflet::tileOptions(opacity = 1)
) %>%
addTiles(urlTemplate = "", attribution = '| Landsat-8 image courtesy of the U.S. Geological Survey')
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] mapview_2.11.2 sf_1.0-15 leafsync_0.1.0 leaflet_2.2.1 rgeeExtra_0.0.1 rgee_1.1.7
## [7] here_1.0.1 fs_1.6.3 rmarkdown_2.25
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 dplyr_1.1.4 fastmap_1.1.1 lazyeval_0.2.2
## [5] leaflet.extras2_1.2.2 webshot2_0.1.1 promises_1.2.1 digest_0.6.34
## [9] lifecycle_1.0.4 geojsonio_0.11.3 ellipsis_0.3.2 processx_3.8.3
## [13] terra_1.7-71 magrittr_2.0.3 compiler_4.3.3 rlang_1.1.3
## [17] tools_4.3.3 utf8_1.2.4 yaml_2.3.8 knitr_1.45
## [21] htmlwidgets_1.6.4 sp_2.1-3 classInt_0.4-10 curl_5.2.0
## [25] reticulate_1.35.0 KernSmooth_2.23-22 websocket_1.4.1 httpcode_0.3.0
## [29] withr_3.0.0 grid_4.3.3 stats4_4.3.3 fansi_1.0.6
## [33] e1071_1.7-14 leafem_0.2.3 colorspace_2.1-0 scales_1.3.0
## [37] crul_1.4.0 cli_3.6.1 crayon_1.5.2 generics_0.1.3
## [41] rstudioapi_0.15.0 geojson_0.3.5 DBI_1.2.2 chromote_0.2.0
## [45] proxy_0.4-27 stringr_1.5.1 base64enc_0.1-3 vctrs_0.6.5
## [49] V8_4.4.2 webshot_0.5.5 Matrix_1.6-5 jsonlite_1.8.8
## [53] geojsonsf_2.0.3 crosstalk_1.2.1 jquerylib_0.1.4 units_0.8-5
## [57] glue_1.7.0 codetools_0.2-19 ps_1.7.6 leaflet.providers_2.0.0
## [61] stringi_1.8.3 later_1.3.2 raster_3.6-26 munsell_0.5.0
## [65] tibble_3.2.1 pillar_1.9.0 htmltools_0.5.7 satellite_1.0.5
## [69] R6_2.5.1 jqr_1.3.3 rprojroot_2.0.4 evaluate_0.23
## [73] lattice_0.22-5 highr_0.10 png_0.1-8 class_7.3-22
## [77] Rcpp_1.0.12 xfun_0.42 pkgconfig_2.0.3
-
アクセスできなくなっています ↩︎
Discussion