- 当Earth Engine用户从教程转向开发自己的处理脚本时,他们会遇到可怕的错误消息,“计算超时”或“用户内存限制超出”。
本节主要内容:
- 了解GEE资源使用的限制条件。
- 熟悉扩展地球引擎操作的多种策略。
- 管理大型项目和多阶段工作流程。
- 识别何时使用Python API可能有利于执行大量任务。
Earth Engine: Under the Hood (GEE :幕后故事)
GEE是一个并行的分布式系统,当你提交任务时,它会将你的查询碎片分解到不同的处理器上,以更有效地完成它们。
Topic 1: Scaling Across Time
1.1 Scaling up with Earth Engine Operators: Annual Daily Climate Data
作为一个例子,我们将提取美国县多边形的降水、最高温度和最低温度的日时间序列。我们将使用 GRIDMET 气候再分析产品(Abatzoglou 2013),该产品提供1979年至今美国连续地区的每日4000米分辨率网格化气象数据。为了节省实习时间,我们将把重点放在美国中部的印第安纳州、伊利诺伊州和爱荷华州,这三个州共包括293个县。
使用 ee.Image.reduceRegions operator
//==========================================================================
//step1:定义featurecall集合,ImageCollection集合和时间序列。
//==========================================================================
// Load county dataset.
// Filter counties in Indiana, Illinois, and Iowa by state FIPS code.
// Select only the unique ID column for simplicity.
var countiesAll = ee.FeatureCollection('TIGER/2018/Counties');
var states = ['17', '18', '19'];
var uniqueID = 'GEOID';
var featColl = countiesAll.filter(ee.Filter.inList('STATEFP', states))
.select(uniqueID);
print(featColl.size());
print(featColl.limit(1));
// Visualize target features (create Figure F6.2.1).
Map.centerObject(featColl, 5);
Map.addLayer(featColl);
// specify years of interest
var startYear = 2020;
var endYear = 2020;
// climate dataset info
var imageCollectionName = 'IDAHO_EPSCOR/GRIDMET';
var bandsWanted = ['pr', 'tmmn', 'tmmx'];
var scale = 4000;
// Load and format climate data.
var startDate = startYear + '-01-01';
var endYear_adj = endYear + 1;
var endDate = endYear_adj + '-01-01';
var imageCollection = ee.ImageCollection(imageCollectionName)
.select(bandsWanted)
.filterBounds(featColl)
.filterDate(startDate, endDate);
// get values at features
var sampledFeatures = imageCollection.map(function(image) {
return image.reduceRegions({
collection: featColl,
reducer: ee.Reducer.mean(),
scale: scale
}).filter(ee.Filter.notNull(
bandsWanted)) // drop rows with no data
.map(function(f) { // add date property
var time_start = image.get( 'system:time_start');
var dte = ee.Date(time_start).format( 'YYYYMMdd');
return f.set('date_ymd', dte);
});
}).flatten();
print(sampledFeatures.limit(1));
// export info
var exportFolder = 'GEE_scalingUp';
var filename = 'Gridmet_counties_IN_IL_IA_' + scale + 'm_' + startYear + '-' + endYear;
// prepare export: specify properties/columns to include
var columnsWanted = [uniqueID].concat(['date_ymd'], bandsWanted);
print(columnsWanted);
Export.table.toDrive({
collection: sampledFeatures,
description: filename,
folder: exportFolder,
fileFormat: 'CSV',
selectors: columnsWanted
});
1.2 Scaling Across Time by Batching: Get 20 Y ears of Daily Climate Data
通过批处理跨时间缩放: 获得20年的每日气候数据
//==========================================================================
//step1:定义featurecall集合,ImageCollection集合和时间序列。
//==========================================================================
// Load county dataset.
// Filter counties in Indiana, Illinois, and Iowa by state FIPS code.
// Select only the unique ID column for simplicity.
var countiesAll = ee.FeatureCollection('TIGER/2018/Counties');
var states = ['17', '18', '19'];
var uniqueID = 'GEOID';
var featColl = countiesAll.filter(ee.Filter.inList('STATEFP', states)) .select(uniqueID);
print(featColl.size());
print(featColl.limit(1));
// Visualize target features (create Figure F6.2.1).
Map.centerObject(featColl, 5);
Map.addLayer(featColl);
// specify years of interest
//var startYear = 2020;
var startYear = 2001;
var endYear = 2020;
//===============================================================================================
// climate dataset info
var imageCollectionName = 'IDAHO_EPSCOR/GRIDMET';
var bandsWanted = ['pr', 'tmmn', 'tmmx'];
var scale = 4000;
//=========================================================================================
// export info
var exportFolder = 'GEE_scalingUp';
var filenameBase = 'Gridmet_counties_IN_IL_IA_' + scale + 'm_';
//==============================================================================================
for (var i = startYear; i <= endYear; i++) { // for each year....
// Load climate collection for that year.
var startDate = i + '-01-01';
var endYear_adj = i + 1;
var endDate = endYear_adj + '-01-01';
var imageCollection = ee.ImageCollection(imageCollectionName)
.select(bandsWanted)
.filterBounds(featColl)
.filterDate(startDate, endDate);
var sampledFeatures = imageCollection.map(function(image) {
return image.reduceRegions({
collection: featColl,
reducer: ee.Reducer.mean(),
tileScale: 1,
scale: scale
}).filter(ee.Filter.notNull(bandsWanted)) // remove rows without data
.map(function(f) { // add date property
var time_start = image.get('system:time_start');
var dte = ee.Date(time_start).format('YYYYMMdd');
return f.set('date_ymd', dte);
});
}).flatten();
// prepare export: specify properties/columns to include
var columnsWanted = [uniqueID].concat(['date_ymd'], bandsWanted);
print(columnsWanted);
var filename = filenameBase + i;
Export.table.toDrive({
collection: sampledFeatures,
description: filename,
folder: exportFolder,
fileFormat: 'CSV',
selectors: columnsWanted
});
}
结果:
2.1 Scaling Across Space via Spatial Tiling
举个例子:基于S2 SR产品合成无云影像。
//===============================step1====================================
//定义感性兴趣区、时间范围、函数
var roi = ee.Geometry.Point([-122.33524518034544, 47.61356183942883]);
// Dates over which to create a median composite.
var start = ee.Date('2019-03-01');
var end = ee.Date('2019-09-01');
// Specify module with cloud mask functions.
var s2mask_tools = require('projects/gee-edu/book:Part F - Fundamentals/F6 - Advanced Topics/F6.2 Scaling Up/modules/s2cloudmask.js');
//===============================step2====================================
//load S2的三种影像产品
// Sentinel-2 surface reflectance data
//'COPERNICUS/S2_SR'已更新为'COPERNICUS/S2_SR_HARMONIZED'
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterDate(start, end)
.filterBounds(roi)
.select(['B2', 'B3', 'B4', 'B5']);
print("s2Sr" ,s2Sr);
// Sentinel-2 Level 1C data (top-of-atmosphere).
//'COPERNICUS/S2'已更新为'COPERNICUS/S2_HARMONIZED'
var s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
.filterBounds(roi)
.filterDate(start, end)
.select(['B7', 'B8', 'B8A', 'B10']);
print("s2" ,s2);
// Cloud probability dataset - used in cloud mask function
var s2c = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
.filterDate(start, end)
.filterBounds(roi);
print("s2c" ,s2c);
//应用云掩膜,将cloud probability dataset 合并到surface reflectance.
var withCloudProbability = s2mask_tools.indexJoin(s2Sr, s2c, 'cloud_probability');
print("withCloudProbability" ,withCloudProbability);
// Join the L1C data to get the bands needed for CDI.
var withS2L1C = s2mask_tools.indexJoin(withCloudProbability, s2, 'l1c');
//对合并的影像都循环使用一下云掩膜函数,将输出结果转换为影像集合。
var masked = ee.ImageCollection(withS2L1C.map(s2mask_tools.maskImage));
//===============================step3====================================
//生成和可视化中值合成影像
// Take the median, specifying a tileScale to avoid memory errors.
var median = masked.reduce(ee.Reducer.median(), 8);
// Display the results.
Map.centerObject(roi, 12);
Map.addLayer(roi);
var viz = {
bands: ['B4_median', 'B3_median', 'B2_median'],
min: 0,
max: 3000
};
Map.addLayer(median, viz, 'median');
结果:
2.2 Generate a Regional Composite Through Spatial Tiling
通过空间平铺生成区域合成影像
// Specify helper functions.
var s2mask_tools = require('projects/gee-edu/book:Part F - Fundamentals/
F6 - Advanced Topics/F6.2 Scaling Up/modules/s2cloudmask.js');
// Set the Region of Interest: Washington, USA.
var roi = ee.FeatureCollection('TIGER/2018/States')
.filter(ee.Filter.equals('NAME', 'Washington'));
// Specify grid size in projection, x and y units (based on projection).
var projection = 'EPSG:4326'; // WGS84 lat lon
var dx = 2.5;
var dy = 1.5;
// Dates over which to create a median composite.
var start = ee.Date('2019-03-01');
var end = ee.Date('2019-09-01');
// Make grid and visualize.
var proj = ee.Projection(projection).scale(dx, dy);
var grid = roi.geometry().coveringGrid(proj);
Map.addLayer(roi, {}, 'roi');
Map.addLayer(grid, {}, 'grid');
创建一个新的空ImageCollection资产,用作导出目的地
//==================================step1============================================
// Specify helper functions.
var s2mask_tools = require('projects/gee-edu/book:Part F - Fundamentals/F6 - Advanced Topics/F6.2 Scaling Up/modules/s2cloudmask.js');
// Set the Region of Interest: Washington, USA.
var roi = ee.FeatureCollection('TIGER/2018/States')
.filter(ee.Filter.equals('NAME', 'Washington'));
// Specify grid size in projection, x and y units (based on projection).
var projection = 'EPSG:4326'; // WGS84 lat lon
var dx = 2.5;
var dy = 1.5;
// Dates over which to create a median composite.
var start = ee.Date('2019-03-01');
var end = ee.Date('2019-09-01');
// Make grid and visualize.
var proj = ee.Projection(projection).scale(dx, dy);
var grid = roi.geometry().coveringGrid(proj);
Map.addLayer(roi, {}, 'roi');
Map.addLayer(grid, {}, 'grid');
// Export info.
var assetCollection = 'projects/ee-maziwei306/assets/S2compositeWA';
var imageBaseName = 'S2_median_';
// Get a list based on grid number.
var gridSize = grid.size().getInfo();
//应该谨慎使用getInfo,并且不要在for循环中使用它
var gridList = grid.toList(gridSize);
//==================================step2============================================
//通过循环批量生成每个网格的合成图像任务导出
// In each grid cell, export a composite
for (var i = 0; i < gridSize; i++) {
// Extract grid polygon and filter S2 datasets for this region.
var gridCell = ee.Feature(gridList.get(i)).geometry();
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterDate(start, end)
.filterBounds(gridCell)
.select(['B2', 'B3', 'B4', 'B5']);
var s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
.filterDate(start, end)
.filterBounds(gridCell)
.select(['B7', 'B8', 'B8A', 'B10']);
var s2c = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
.filterDate(start, end)
.filterBounds(gridCell);
// Apply the cloud mask.
var withCloudProbability = s2mask_tools.indexJoin(s2Sr, s2c, 'cloud_probability');
var withS2L1C = s2mask_tools.indexJoin(withCloudProbability, s2, 'l1c');
var masked = ee.ImageCollection(withS2L1C.map(s2mask_tools.maskImage));
// Generate a median composite and export.
var median = masked.reduce(ee.Reducer.median(), 8);
// Export.
var imagename = imageBaseName + 'tile' + i;
Export.image.toAsset({
image: median,
description: imagename,
assetId: assetCollection + '/' + imagename,
scale: 10,
region: gridCell,
maxPixels: 1e13
});
}
//==================================step3============================================
// 加载图像采集和拼接成单个图像
var assetCollection = 'projects/ee-maziwei306/assets/S2compositeWA';
var composite = ee.ImageCollection(assetCollection).mosaic();
// Display the results
var geometry = ee.Geometry.Point([-120.5873563817392, 47.39035206888694 ]);
Map.centerObject(geometry, 6);
var vizParams = {
bands: ['B4_median', 'B3_median', 'B2_median'],
min: 0,
max: 3000
};
Map.addLayer(composite, vizParams, 'median');
3.1 Multistep Workflows and Intermediate Assets
多步骤工作流和中间资产
通常,我们的目标需要几个处理步骤,不能在一个地球引擎计算链中完成。在这些情况下,最好的策略是将任务分解成单独的部分,这些部分被创建、存储在资产中,并跨多个脚本使用。每个顺序脚本创建一个中间输出,这个中间输出将成为下一个脚本的输入。
举个栗子:
考虑确定灌溉农业用地的土地利用分类任务:包括要用到卫星数据、综合天气信息、土壤信息 作物类型位置等。可以按以下几个步骤去实现:
- 生成年或月NDVi的卫星合成数据;
- 将气候数据处理成月或季度值;
- 从地面真值层生成随机点位置,用作特征训练数据集和准确性验证,并提取这些特征的复合值和天气值。
- 训练分类器并应用它,可能跨越多年;研究人员经常会实现多个分类器,并比较不同方法的性能。
- 执行分类后的后处理步骤,例如去除“斑点”。
- 使用每个行政边界的总面积来评估地面真相验证点和政府统计数据的准确性。
- 将结果导出为空间层、可视化或其他格式。