Google Earth Engine（GEE）——R 语言 Google 地球引擎20个基本案例分析
基本 rgee - 最佳实践
本文档描述了旨在最大化复杂或昂贵的地球引擎计算成功机会的编码实践。
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() 的函数中不起作用。例如：
错误— 此代码不起作用！
-
table <- ee$FeatureCollection('USDOS/LSIB_SIMPLE/2017')
-
# Error:
-
foobar <- table$map(function(f) {
-
print(f); # Can't use a client function here.
-
# Can't Export, either.
-
})
好- 使用map() set().
-
table <- ee$FeatureCollection('USDOS/LSIB_SIMPLE/2017')
-
# 对集合的每个元素做一些事情。
-
withMoreProperties = table$map(function(f) {
-
# 设置属性。
-
f$set("area_sq_meters", f$area())
-
})
-
print(withMoreProperties$first()$get("area_sq_meters")$getInfo())
您还可以filter()基于计算或现有属性和print()结果的集合。请注意，您无法打印包含超过 5000 个元素的集合。如果您收到“Collection query aborted after accumulating over 5000 elements”错误，filter()或者limit()在打印之前收集。
2. 避免不必要地转换为列表
Earth Engine 中的集合使用优化进行处理，这些优化通过将集合转换为 List或Array类型而被破坏。除非您需要随机访问集合元素（即您需要获取集合的第 i 个元素），否则请在集合上使用过滤器来访问单个集合元素。以下示例说明了类型转换（不推荐）和过滤（推荐）以访问集合中的元素之间的区别：
坏- 不要不必要地转换为列表！
-
table <- ee$FeatureCollection('USDOS/LSIB_SIMPLE/2017');
-
# 不要这么做
-
list <- table$toList(table$size())
-
print(list$get(13)$getInfo()) # 超出用户内存限制。
请注意，您可以通过将集合不必要地转换为列表来轻松触发错误。更安全的方法是使用filter()：
好- 使用filter().
print(table$filter(ee$Filter$eq('country_na', 'Niger'))$first()$getInfo())
请注意，您应该在分析中尽早使用过滤器。
3. 避免ee.Algorithms.If()
不要ee.Algorithms.If()用于实现分支逻辑，尤其是在映射函数中。如以下示例所示，ee.Algorithms.If()可能会占用大量内存，不推荐使用：
坏- 不要使用If()：
-
table <- ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
-
# 不要这样做！
-
veryBad = table$map(function(f) {
-
ee$Algorithms$If(
-
condition = ee$String(f$get('country_na'))$compareTo('Chad')$gt(0),
-
trueCase = f, #做点什么。
-
falseCase = NULL # 做点别的。
-
)
-
}, TRUE)
-
print(veryBad$getInfo()) # User memory limit exceeded.
-
# If() 可以评估 true 和 false 两种情况。
请注意，的第二个参数map()是TRUE。这意味着映射的函数可能会返回空值，并且它们将被删除到结果集合中。这可能很有用（没有If()），但这里最简单的解决方案是使用过滤器：
好- 使用filter().
print(table$filter(ee$Filter$eq('country_na', 'Chad')))
如本教程所示，使用过滤器的函数式编程方法是将一种逻辑应用于集合的某些元素并将另一种逻辑应用于集合的其他元素的正确方法。
4. 避免reproject()
reproject除非绝对必要，否则不要使用。您可能想要使用 reproject() 的一个原因是强制Map显示计算以特定比例发生，以便您可以在所需的分析比例下检查结果。在下一个示例中，计算热点像素的补丁并计算每个补丁中的像素数。运行示例并单击其中一个补丁。请注意，重新投影的数据和未重新投影的数据之间的像素数不同。
-
l8sr <- ee$ImageCollection("LANDSAT/LC08/C01/T1_SR")
-
sf <- ee$Geometry$Point(c(-122.405, 37.786))
-
Map$centerObject(sf, 13)
-
# 重新投影的原因 - 计算像素并以交互方式探索。
-
image <- l8sr$filterBounds(sf)$
-
filterDate("2019-06-01", "2019-12-31")$
-
first()
-
Map$addLayer(image, list(bands = "B10", min = 2800, max = 3100), "image")
-
hotspots <- image$select("B10")$
-
gt(3100)$
-
selfMask()$
-
rename("hotspots")
-
objectSize <- hotspots$connectedPixelCount(256)
-
# 小心重新规划！不要缩小重新投影的数据。
-
reprojected <- objectSize$reproject(hotspots$projection())
-
Map$addLayer(objectSize, list(min = 1, max = 256), "Size No Reproject", FALSE) +
-
Map$addLayer(reprojected, list(min = 1, max = 256), "Size Reproject", FALSE)
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')
在这个例子中，请注意均值归约器与标准偏差归约器相结合，并且sharedInputs能够单次通过输入像素。在输出字典中，reducer 的名称附加到带名称。要获得均值和 SD 图像（例如对输入图像进行归一化），您可以将值转换为图像并使用正则表达式分别提取均值和 SD，如示例中所示。
8. 使用导出
对于在代码编辑器中导致“超出用户内存限制”或“计算超时”错误的计算，使用Export. 这是因为在批处理系统（导出运行的地方）中运行时，超时时间更长，并且允许的内存占用量更大。（您可能想先尝试其他方法，如调试文档中所述）。继续前面的例子，假设字典返回了一个错误。您可以通过执行以下操作来获得结果：
-
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)
请注意，该链接已嵌入到资产名称中，以实现可重复性。另请注意，如果您想导出toAsset，您将需要提供一个几何图形，它可以是任何东西，例如图像质心，它很小且计算成本低。（即，如果您不需要，请不要使用复杂的几何图形）。
有关使用Export解决Computation timed out 和Too many concurrent aggregations 的示例，请参阅调试页面。有关一般导出的详细信息，请参阅此文档。
-
table <- ee$FeatureCollection('USDOS/LSIB_SIMPLE/2017')
-
l8sr <- ee$ImageCollection('LANDSAT/LC08/C01/T1_SR')
-
chad <- table$filter(ee$Filter$eq('country_na', 'Chad'))$first()
-
# 最好避免不必要的裁剪
-
unnecessaryClip <- l8sr$
-
select('B4')$ # Good.
-
filterBounds(chad$geometry())$ # Good.
-
filterDate('2019-01-01', '2019-12-31')$ # Good.
-
map(function(image) {
-
image$clip(chad$geometry()) # NO! Bad! Not necessary.
-
})$
-
median()$
-
reduceRegion(
-
reducer = ee$Reducer$mean(),
-
geometry = chad$geometry(),
-
scale = 30,
-
maxPixels = 1e10
-
)
-
print(unnecessaryClip$getInfo())
可以完全跳过剪切输入图像，因为在reduceRegion()调用中指定了区域：
好- 在输出上指定区域。
-
noClipNeeded <- l8sr$
-
select('B4')$ # Good.
-
filterBounds(chad$geometry())$ # Good.
-
filterDate('2019-01-01', '2019-12-31')$ # Good.
-
median()$
-
reduceRegion(
-
reducer = ee$Reducer$mean(),
-
geometry = chad$geometry(), # Geometry is specified here.
-
scale = 30,
-
maxPixels = 1e10
-
)
-
print(noClipNeeded$getInfo())
如果此计算超时，Export则如本例所示。
10.如果需要剪辑复杂的集合，使用clipToCollection()
如果您确实需要剪辑某些内容，并且要用于剪辑的几何图形位于集合中，请使用clipToCollection()：
-
ecoregions <- ee$FeatureCollection('RESOLVE/ECOREGIONS/2017')
-
image <- ee$Image('JAXA/ALOS/AW3D30_V1_1')
-
complexCollection <- ecoregions$
-
filter(
-
ee$Filter$eq(
-
'BIOME_NAME',
-
'Tropical & Subtropical Moist Broadleaf Forests'
-
)
-
)
-
Map$addLayer(complexCollection, {}, 'complexCollection')
-
clippedTheRightWay <- image$select('AVE')$
-
clipToCollection(complexCollection)
-
Map$addLayer(clippedTheRightWay, {}, 'clippedTheRightWay', FALSE)
12. 使用非零的errorMargin
对于可能昂贵的几何运算，在给定计算精度的情况下，尽可能使用最大的误差容限。误差幅度指定在几何操作期间（例如在重新投影期间）允许的最大允许误差（以米为单位）。指定较小的误差幅度可能会导致需要对几何图形（带坐标）进行密集化，这可能会占用大量内存。为您的计算指定尽可能大的误差范围是一种很好的做法：
-
ecoregions <- ee$FeatureCollection("RESOLVE/ECOREGIONS/2017")
-
complexCollection <- ecoregions$limit(10)
-
Map$centerObject(complexCollection)
-
Map$addLayer(complexCollection)
-
expensiveOps <- complexCollection$map(function(f) {
-
f$buffer(10000, 200)$bounds(200)
-
})
-
Map$addLayer(expensiveOps, {}, 'expensiveOps')
13. 不要在reduceToVectors() 中使用的规模
如果要将栅格转换为矢量，请使用适当的比例。指定一个非常小的比例会大大增加计算成本。将比例设置得尽可能高，以获得所需的精度。例如，要获取代表全球陆地的多边形：
-
etopo <- ee$Image('NOAA/NGDC/ETOPO1')
-
# 大致的陆地边界。
-
bounds <- etopo$select(0)$gt(-100)
-
# 非测地线多边形。
-
almostGlobal <- ee$Geometry$Polygon(
-
coords = list(
-
c(-180, -80),
-
c(180, -80),
-
c(180, 80),
-
c(-180, 80),
-
c(-180, -80)
-
),
-
proj = "EPSG:4326",
-
geodesic = FALSE
-
)
-
Map$addLayer(almostGlobal, {}, "almostGlobal")
-
vectors <- bounds$selfMask()$reduceToVectors(
-
reducer = ee$Reducer$countEvery(),
-
geometry = almostGlobal,
-
# 将比例设置为给定的最大可能
-
# 计算所需的精度。
-
scale = 50000
-
)
-
Map$addLayer(vectors, {}, "vectors")
在前面的示例中，请注意在全局缩减中使用非测地线多边形。
14. 不要将reduceToVectors()与reduceRegions() 一起使用
不要使用FeatureCollection返回的reduceToVectors()作为 的输入reduceRegions()。相反，在调用之前添加要减少的波段reduceToVectors()：
-
etopo <- ee$Image('NOAA/NGDC/ETOPO1')
-
mod11a1 <- ee$ImageCollection('MODIS/006/MOD11A1')
-
# 大致的陆地边界。
-
bounds <- etopo$select(0)$gt(-100)
-
# 非测地线多边形。
-
almostGlobal <- ee$Geometry$Polygon(
-
coords = list(c(-180, -80), c(180, -80), c(180, 80), c(-180, 80), c(-180, -80)),
-
proj = "EPSG:4326",
-
geodesic = FALSE
-
)
-
lst <- mod11a1$first()$select(0)
-
means <- bounds$selfMask()$addBands(lst)$reduceToVectors(
-
reducer = ee$Reducer$mean(),
-
geometry = almostGlobal,
-
scale = 1000,
-
maxPixels = 1e10
-
)
-
print(means$limit(10)$getInfo())
请注意，减少另一个区域内一个图像的像素的其他方法包括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
-
)
-
# 导出大样本...
在此示例中，请注意图像导出为浮点数。除非绝对必要，否则不要以双精度导出。
导出完成后，重新加载资产并继续从中采样。请注意，首先在非常小的测试区域上运行非常小的样本，以进行调试。当证明成功时，获取更大的样本并将其导出。如此大的样本通常需要出口。不要期望这些样本在print()没有先导出它们的情况下可以交互（例如通过）或可用（例如作为分类器的输入）。
19.加入与地图过滤器
假设您想根据时间、位置或某些元数据属性加入集合。通常，这是通过连接最有效地完成的。以下示例在 Landsat 8 和 Sentinel-2 集合之间进行时空连接：
-
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())
尽管您应该首先尝试连接（Export如果需要），但有时 afilter()内的 amap()也可能有效，尤其是对于非常大的集合。
-
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())
20. reduceRegion() vs. reduceRegions() vs. for-loop
reduceRegions()使用非常大或复杂的 FeatureCollection输入进行调用可能会导致可怕的“计算值太大”错误。一个可能的解决方案是映射reduceRegion() 在FeatureCollection代替。另一个潜在的解决方案是使用 (gasp) for 循环。尽管所描述的这种强烈的地球引擎气馁这里，这里和这里，reduceRegion()可以在一个for循环来实现，以执行大量减少。
假设你的目的是获得的像素中每个特征的平均（或任何统计量）在一个FeatureCollection用于在每个图像ImageCollection。以下示例比较了前面描述的三种方法：
-
# 选择国家也就是矢量边界
-
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())
请注意first()，出于调试目的，会打印每个集合中的事物。您不应该期望以交互方式提供完整的结果：您需要Export. 另请注意，应极其谨慎地使用 for 循环，并且仅作为最后的手段。最后，for 循环需要手动获取输入集合的大小并将其硬编码到适当的位置。如果其中任何一个听起来不清楚，请不要使用 for 循环。
21. 及时对邻居使用前向差分
假设您有一个时间排序ImageCollection（即时间序列），并且您想将每个图像与前一个（或下一个）图像进行比较。与其iterate()为此目的而使用，使用基于数组的前向差分可能更有效。以下示例使用此方法对 Sentinel-2 集合进行重复数据删除，其中重复项定义为一年中同一天的图像：
-
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())
文章来源: blog.csdn.net，作者：此星光明2021年博客之星云计算Top3，版权归原作者所有，如需转载，请联系作者。
原文链接：blog.csdn.net/qq_31988139/article/details/119968371
- 点赞
- 收藏
- 关注作者
评论（0）