Google Earth Engine（GEE）——R 语言 Google 地球引擎20个基本案例分析_gee地图制作案例分析

1. 避免将客户端函数和对象与服务器函数和对象混合

Earth Engine 服务器对象是具有以ee（例如ee$Image、ee$Reducer）开头的构造函数的对象，并且此类对象上的任何方法都是服务器函数。任何不是以这种方式构造的对象都是客户端对象。客户端对象可能来自 R Earth Engine 客户端（例如 Map）或 R 语言（例如 date、data.frame、c()、list()）。

为避免意外行为，请勿在脚本中混合使用客户端和服务器功能，如此处、此处和此处讨论的那样。有关地球引擎中客户端与服务器的深入解释，请参阅此页面和/或本教程。以下示例说明了混合客户端和服务器功能的危险：

   错误— 此代码不起作用！

# 这种状态不会执行
for (i in seq_len(table$size())) {
  print('No!') 
}

你能发现错误吗？请注意，这table$size()是服务器对象上的服务器方法，不能与客户端功能（如seq_len函数）一起使用。

您可能希望使用 for 循环的一种情况是 with 显示结果Map，因为 Map 对象和方法是客户端。

   好- 使用客户端函数显示地球引擎空间对象。

l8_ts <- sprintf(
  "LANDSAT/LC08/C01/T1/LC08_044034_%s",
  c("20140318", "20140403","20140419","20140505")
)
display_l8ts <- list()
for (l8 in l8_ts) {
  ee_l8 <- ee$Image(l8)
  display_l8ts[[l8]] <- Map$addLayer(ee_l8)
}
Map$centerObject(ee_l8)
Reduce('+', display_l8ts)

相反，map()是服务器函数和客户端功能在传递给 map() 的函数中不起作用。例如：

   错误— 此代码不起作用！

出现差异的原因是分析的比例是由代码编辑器缩放级别设置的。通过调用reproject()您可以设置计算的比例而不是地图显示。reproject()出于本文档中描述的原因，请谨慎使用。

5.过滤和选择（）第一

通常，在对集合执行任何其他操作之前，按时间、位置和/或元数据过滤输入集合。在选择性较少的过滤器之前应用更多选择性过滤器。空间和/或时间过滤器通常更具选择性。例如，请注意select()和filter()之前应用map()：

images <- ee$ImageCollection("COPERNICUS/S2_SR")
sf <- ee$Geometry$Point(c(-122.463, 37.768))
# 减少图像邻域的昂贵功能。
reduceFunction <- function(image) {
  image$reduceNeighborhood(
    reducer = ee$Reducer$mean(),
    kernel = ee$Kernel$square(4)
  )
}
bands <- list("B4", "B3", "B2")
# Select and filter first!
reasonableComputation <- images$select(bands)$
  filterBounds(sf)$
  filterDate("2018-01-01", "2019-02-01")$
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", 1))$
  aside(ee_print)$ # Useful for debugging.
  map(reduceFunction)$
  reduce('mean')$
  rename(bands)
viz <- list(bands = bands, min = 0, max = 10000)
Map$addLayer(reasonableComputation, viz, "resonableComputation")

6. 使用updateMask()而不是mask()这个问题之前再JS中讲过，更新过后的掩膜更加准确

updateMask()和之间的区别在于mask()前者and()对参数（新掩膜）和现有图像掩膜进行逻辑处理，而mask()只是用参数替换图像掩膜。后者的危险在于您可能会无意中取消屏蔽像素。在此示例中，目标是屏蔽小于或等于 300 米高程的像素。正如您所看到的（缩小），使用mask()会导致很多像素被掩盖，这些像素不属于感兴趣的图像：

l8sr <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR")
sf <- ee$Geometry$Point(c(-122.40554461769182, 37.786807309873716))
aw3d30 <- ee$Image("JAXA/ALOS/AW3D30_V1_1")
Map$centerObject(sf, 7)
image <- l8sr$filterBounds(sf)$filterDate("2019-06-01", "2019-12-31")$first()
vis <- list(bands = c("B4", "B3", "B2"), min = 0, max = 3000)
Map$addLayer(image, vis, "image", FALSE)
mask <- aw3d30$select("AVE")$gt(300)
Map$addLayer(mask, {}, 'mask', FALSE)
# 两者的最直接区别
badMask <- image$mask(mask)
Map$addLayer(badMask, vis, "badMask")
goodMask <- image.updateMask(mask)
Map$addLayer(goodMask, vis, "goodMask", FALSE)

7.组合reducer

如果您需要来自单个输入（例如图像区域）的多个统计信息（例如均值和标准差），则组合减速器会更有效。例如，要获取图像统计信息，请按如下方式组合 reducer：

image <- ee$Image('COPERNICUS/S2/20150821T111616_20160314T094808_T30UWU')
# 通过组合减速器获得每个波段的均值和标准差。
stats <- image$reduceRegion(
  reducer = ee$Reducer$mean()$combine(
    reducer2 = ee$Reducer$stdDev(),
    sharedInputs = TRUE
  ),
  geometry = ee$Geometry$Rectangle(c(-2.15, 48.55, -1.83, 48.72)),
  scale = 10,
  bestEffort = TRUE # Use maxPixels if you care about scale.
)
print(stats$getInfo())
# 将平均值和 SD 提取到图像。
meansImage <- stats$toImage()$select('.*_mean')
sdsImage <- stats$toImage()$select('.*_stdDev')

link <- '86836482971a35a5e735a17e93c23272'
task <- ee$batch$Export$table$toDrive(
  collection = ee$FeatureCollection(ee$Feature(NULL, stats)),
  description = paste0("exported_stats_demo_", link),
  fileFormat = "CSV"
)
# Using rgee I/O
task <- ee_table_to_drive(
  collection = ee$FeatureCollection(ee$Feature(NULL, stats)),
  description = paste0("exported_stats_demo_", link),
  fileFormat = "CSV"
)
task$start()
ee_monitoring(task)
exported_stats <- ee_drive_to_local(task = task,dsn = "exported_stats.csv")
read.csv(exported_stats)

请注意，减少另一个区域内一个图像的像素的其他方法包括reduceConnectedComponents()和/或分组 reducers。

15.使用fastDistanceTransform()进行邻域操作

对于某些卷积运算，fastDistanceTransform()可能比reduceNeighborhood()或更有效convolve()。例如，要对二进制输入进行腐蚀和/或膨胀：

aw3d30 <- ee$Image("JAXA/ALOS/AW3D30_V1_1")
# 从高程阈值制作一个简单的二元层。
mask <- aw3d30$select("AVE")$gt(300)
Map$setCenter(-122.0703, 37.3872, 11)
Map$addLayer(mask, {}, "mask")
# 以像素为单位的距离。
distance <- mask$fastDistanceTransform()$sqrt()
# 膨胀的距离阈值（三个像素）。
dilation <- distance$lt(3)
Map$addLayer(dilation, {}, "dilation")
# 对腐蚀进行相反的操作。
notDistance <- mask$Not()$fastDistanceTransform()$sqrt()
erosion <- notDistance$gt(3)
Map$addLayer(erosion, {}, 'erosion')

16. 使用reduceNeighborhood() 中的优化

如果您需要执行卷积并且不能使用fastDistanceTransform()，请使用 中的优化reduceNeighborhood()。

l8raw <- ee$ImageCollection('LANDSAT/LC08/C01/T1_RT')
composite <- ee$Algorithms$Landsat$simpleComposite(l8raw)
bands <- c('B4', 'B3', 'B2')
optimizedConvolution <- composite$select(bands)$reduceNeighborhood(
  reducer = ee$Reducer$mean(),
  kernel = ee$Kernel$square(3),
  optimization = "boxcar" # Suitable optimization for mean.
)$rename(bands)
viz <- list(bands = bands, min = 0, max = 72)
Map$setCenter(-122.0703, 37.3872, 11)
Map$addLayer(composite, viz, "composite") +
Map$addLayer(optimizedConvolution, viz, "optimizedConvolution")

17. 不要采样超过你需要的数据

抵制不必要地增加训练数据集大小的冲动。尽管在某些情况下增加训练数据量是一种有效的机器学习策略，但它也会增加计算成本，而不会相应提高准确性。（要了解何时增加训练数据集大小，请参阅此参考资料）。以下示例演示了请求过多训练数据如何导致可怕的“计算值太大”错误：

   不好——不要采样太多数据！

l8raw <- ee$ImageCollection('LANDSAT/LC08/C01/T1_RT')
composite <- ee$Algorithms$Landsat$simpleComposite(l8raw)
labels <- ee$FeatureCollection('projects/google/demo_landcover_labels')
# 不！没有必要。不要这样做：
labels <- labels$map(function(f){f$buffer(100000, 1000)})
bands <- c('B2', 'B3', 'B4', 'B5', 'B6', 'B7')
training <- composite$select(bands)$sampleRegions(
  collection = labels,
  properties = list("landcover"),
  scale = 30
)
classifier <- ee$Classifier$smileCart()$train(
  features = training,
  classProperty = "landcover",
  inputProperties = bands
)
print(classifier$explain()) # Computed value is too large

更好的方法是从适量的数据开始并调整分类器的超参数以确定是否可以达到所需的准确度：

   好——调整超参数。

l8raw <- ee$ImageCollection("LANDSAT/LC08/C01/T1_RT")
composite <- ee$Algorithms$Landsat$simpleComposite(l8raw)
labels <- ee$FeatureCollection("projects/google/demo_landcover_labels")
# 稍微增加数据，可能会引入噪音。
labels <- labels$map(function(f) {f$buffer(100, 10)})
bands <- c('B2', 'B3', 'B4', 'B5', 'B6', 'B7')
data <- composite$select(bands)$sampleRegions(
  collection = labels,
  properties = list("landcover"),
  scale = 30
)
# 添加一列名为“随机”的统一随机数。
data <- data$randomColumn()
# 划分为训练和测试。
training <- data$filter(ee$Filter$lt("random", 0.5))
testing <- data$filter(ee$Filter$gte("random", 0.5))
# 调整 minLeafPopulation 参数。
minLeafPops <- ee$List$sequence(1, 10)
accuracies <- minLeafPops$map(
  ee_utils_pyfunc(
    function(p) {
      classifier <- ee$Classifier$smileCart(minLeafPopulation = p)$
        train(
          features = training,
          classProperty = "landcover",
          inputProperties = bands
        )
      
      testing$
        classify(classifier)$
        errorMatrix("landcover", "classification")$
        accuracy()
    }
  )
)
minLeafPopulation_array <- accuracies$getInfo()
plot(
  x = minLeafPopulation_array,
  type = "b", 
  col = "blue",
  lwd = 2,
  ylab = "accuracy",
  xlim = c(0,10),
  xlab = "value",
  main = "Hyperparameter tunning (minLeafPopulation)"
)

在这个例子中，分类器已经非常准确，所以没有太多的调整要做。您可能希望选择minLeafPopulation仍然具有所需精度的最小树（即最大）。

18.导出中间结果

假设您的目标是从相对复杂的计算图像中取样。通常Export对图像更有效toAsset()，加载导出的图像，然后采样。例如：

image <- ee$Image('UMD/hansen/global_forest_change_2018_v1_6')
geometry <- ee$Geometry$Polygon(
  coords = list(
    c(-76.64069800085349, 5.511777325802095),
    c(-76.64069800085349, -20.483938229362376),
    c(-35.15632300085349, -20.483938229362376),
    c(-35.15632300085349, 5.511777325802095)
  ),
  proj =  "EPSG:4326",
  geodesic =  FALSE
)
testRegion <- ee$Geometry$Polygon(
  coords = list(
    c(-48.86726050085349, -3.0475996402515717),
    c(-48.86726050085349, -3.9248707849303295),
    c(-47.46101050085349, -3.9248707849303295),
    c(-47.46101050085349, -3.0475996402515717)
  ),
  proj = "EPSG:4326",
  geodesic = FALSE
)
# 2016 年森林损失，对样本进行分层。
loss <- image$select("lossyear")
loss16 <- loss$eq(16)$rename("loss16")
# 云掩膜功能。
maskL8sr <- function(image) {
  cloudShadowBitMask <- bitwShiftL(1, 3)
  cloudsBitMask <- bitwShiftL(1, 5)
  qa <- image$select('pixel_qa')
  mask <- qa$bitwiseAnd(cloudShadowBitMask)$eq(0)$
    And(qa$bitwiseAnd(cloudsBitMask)$eq(0))
  
  image$updateMask(mask)$
    divide(10000)$
    select("B[0-9]*")$
    copyProperties(image, list("system:time_start"))
}
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR")$map(maskL8sr)
# 创建两个年度无云复合材料。
composite1 <- collection$filterDate('2015-01-01', '2015-12-31')$median()
composite2 <- collection$filterDate('2017-01-01', '2017-12-31')$median()
# 我们想要这个堆栈的分层样本。
stack <- composite1$addBands(composite2)$float() # Export the smallest size possible.
#导出图像。由于导出已完成，因此对该块进行了注释。
# link <- "0b8023b0af6c1b0ac7b5be649b54db06"
# desc <- paste0(ee_get_assethome(), "/Logistic_regression_stack_", link)
# 
# #ee_image_info(stack)
# task <- ee_image_to_asset(
#   image = stack,
#   description = link,
#   assetId = desc,
#   region = geometry,
#   scale = 100,
#   maxPixels = 1e10
# )
  
# 加载导出的图像。
exportedStack <- ee$Image(
  "projects/google/Logistic_regression_stack_0b8023b0af6c1b0ac7b5be649b54db06"
)
# 先取一个很小的样本，进行调试。
testSample <- exportedStack$addBands(loss16)$stratifiedSample(
  numPoints = 1,
  classBand = "loss16",
  region = testRegion,
  scale = 30,
  geometries = TRUE
)
print(testSample$getInfo()) # Check this in the console.
# 取一个大样本。
sample <- exportedStack$addBands(loss16)$stratifiedSample(
  numPoints = 10000,
  classBand = "loss16",
  region = geometry,
  scale = 30
)
# 导出大样本...

s2 <- ee$ImageCollection("COPERNICUS/S2")$
  filterBounds(ee$Geometry$Point(c(-2.0205, 48.647)))
l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR")
joined <- ee$Join$saveAll("landsat")$apply(
  primary = s2,
  secondary = l8,
  condition = ee$Filter$And(
    ee$Filter$maxDifference(
      difference = 1000 * 60 * 60 * 24, # One day in milliseconds
      leftField = "system:time_start",
      rightField = "system:time_start"
    ),
    ee$Filter$intersects(
      leftField = ".geo",
      rightField = ".geo"
    )
  )
)
print(joined$first()$getInfo())

s2 <- ee$ImageCollection("COPERNICUS/S2")$
  filterBounds(ee$Geometry$Point(c(-2.0205, 48.647)))
l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR")
mappedFilter <- s2$map(function(image) {
  date <- image$date()
  landsat <- l8$
    filterBounds(image$geometry())$
    filterDate(date$advance(-1, "day"), date$advance(1, "day"))
    # Return the input image with matching scenes in a property.
  image$set(
    list(
      landsat = landsat,
      size = landsat$size()
    )
  )
})$filter(ee$Filter$gt("size", 0))
print(mappedFilter$first()$getInfo())

# 选择国家也就是矢量边界
countriesTable <- ee$FeatureCollection("USDOS/LSIB_SIMPLE/2017")
# 选择影像
mod13a1 <- ee$ImageCollection("MODIS/006/MOD13A1")
# MODIS vegetation indices (always use the most recent version).
band <- "NDVI"
imagery <- mod13a1$select(band)
# 区域归约操作1: reduceRegions()
testTable <- countriesTable$limit(1) # Do this outside map()s and loops.
data <- imagery$map(function(image) {
  image$reduceRegions(
    collection = testTable,
    reducer = ee$Reducer$mean(),
    scale = 500
  )$map(function(f) {
    f$set(
      list(
        time = image$date()$millis(),
        date = image$date()$format()
      )
    )
  })
})$flatten()
print(data$first()$getInfo())
# 区域归约操作2: mapped reduceRegion()
data <- countriesTable$map(function(feature) {
  imagery$map(
    function(image) {
      ee$Feature(
        feature$geometry()$centroid(100),
        image$reduceRegion(
          reducer = ee$Reducer$mean(),
          geometry = feature$geometry(),
          scale = 500
        )
      )$set(
        list(
          time = image$date()$millis(),
          date = image$date()$format()
        )
      )$copyProperties(feature)
    }
  )
})$flatten()
print(data$first()$getInfo())
# 区域归约操作3: for-loop (WATCH OUT!)
size <- countriesTable$size()
print(size$getInfo()) # 312
countriesList <- countriesTable$toList(1) # Adjust size.
data <- ee$FeatureCollection(list()) # Empty table.
for (j in (seq_len(countriesList$length()$getInfo()) - 1)) {
  feature <- ee$Feature(countriesList$get(j))
  # 转换 ImageCollection > FeatureCollection
  fc <- ee$FeatureCollection(
    imagery$map(
      function(image) {
        ee$Feature(
          feature$geometry()$centroid(100),
          image$reduceRegion(
            reducer = ee$Reducer$mean(),
            geometry = feature$geometry(),
            scale = 500
          )
        )$set(
          list(
            time = image$date()$millis(),
            date = image$date()$format()
          )
        )$copyProperties(feature)
      }
    )
  )
  data <- data$merge(fc)
}
print(data$first()$getInfo())

sentinel2 <- ee$ImageCollection("COPERNICUS/S2")
sf <- ee$Geometry$Point(c(-122.47555371521855, 37.76884708376152))
s2 <- sentinel2$
  filterBounds(sf)$
  filterDate("2018-01-01", "2019-12-31")
withDoys <- s2$map(function(image) {
  ndvi <- image$normalizedDifference(c("B4", "B8"))$rename("ndvi")
  date <- image$date()
  doy <- date$getRelative("day", "year")
  time <- image$metadata("system:time_start")
  doyImage <- ee$Image(doy)$
    rename("doy")$
    int()
  
  ndvi$
    addBands(doyImage)$
    addBands(time)$
    clip(image$geometry()) # Appropriate use of clip.
})
array <- withDoys$toArray()
timeAxis <- 0
bandAxis <- 1
dedup <- function(array) {
  time <- array$arraySlice(bandAxis, -1)
  sorted <- array$arraySort(time)
  doy <- sorted$arraySlice(bandAxis, -2, -1)
  left <- doy$arraySlice(timeAxis, 1)
  right <- doy$arraySlice(timeAxis, 0, -1)
  mask <- ee$Image(ee$Array(list(list(1))))$
    arrayCat(left$neq(right), timeAxis)
  array$arrayMask(mask)
}
deduped <- dedup(array)
# 检查这些输出以确认已删除重复项。
print(array$reduceRegion("first", sf, 10)$getInfo())
print(deduped$reduceRegion("first", sf, 10)$getInfo())