GEE C24 区域统计
导读：
1.以点为中心缓冲成方形或者圆形；
2.编写和应用带有可选参数的函数；
3.学习什么是区域统计并使用reducer;
4.将计算结果导出到表中;
5.将属性从一个图像复制到另一个图像。
一、原理
1.1从给定区域内的栅格计算统计数据称为区域统计。
1.2在第22章和23章中使用ee.Image.ReduceRegions计算分区统计。本章提出了一个更通用的方式——自定义函数——计算区域统计，适用单幅的影像或者影像集合。（ee.Image 和ee.ImageCollection objects.）这种通用的方式除了比较灵活之外，也不容易出现“计算值太大”的错误。（当使用ReduceRegions处理比较复杂或者数据较大的影像集合时通常会出现“计算值太大”的问题）
二、函数
2.1 bufferPoints()参数
function bufferPoints(radius, bounds) {
return function(pt) {
pt = ee.Feature(pt);
return bounds ? pt.buffer(radius).bounds() :
pt.buffer(radius);
};
}
- var ptsTopo = pts.map(bufferPoints(45, false)); //缓冲半径为45，缓冲的形状为圆形
- var ptsTopo1 = pts.map(bufferPoints(45, true));//缓冲半径为45，缓冲的形状为正方形
2.2 zonalStats（）参数
function zonalStats(ic, fc, params) { // 初始化内部参数字典。 var _params = { reducer: ee.Reducer.mean(), scale: null, crs: null, bands: null, bandsRename: null, imgProps: null, imgPropsRename: null, datetimeName: 'datetime', datetimeFormat: 'YYYY-MM-dd HH:mm:ss' }; // 用提供的参数替换初始化的参数。 if (params) { for (var param in params) { _params[param] = params[param] || _params[param]; } } // 根据图像代表设置默认参数. var imgRep = ic.first(); var nonSystemImgProps = ee.Feature(null) .copyProperties(imgRep).propertyNames(); if (!_params.bands) _params.bands = imgRep.bandNames(); if (!_params.bandsRename) _params.bandsRename = _params.bands; if (!_params.imgProps) _params.imgProps = nonSystemImgProps; if (!_params.imgPropsRename) _params.imgPropsRename = _params .imgProps; // 对影像集合迭代使用reduceRegions函数. var results = ic.map(function(img) { // 选择波段 (可选重命名)，设置日期时间和时间戳属性。. img = ee.Image(img.select(_params.bands, _params .bandsRename)) // Add datetime and timestamp features. .set(_params.datetimeName, img.date().format( _params.datetimeFormat)) .set('timestamp', img.get('system:time_start')); // Define final image property dictionary to set in output features. var propsFrom = ee.List(_params.imgProps) .cat(ee.List([_params.datetimeName, 'timestamp'])); var propsTo = ee.List(_params.imgPropsRename) .cat(ee.List([_params.datetimeName, 'timestamp'])); var imgProps = img.toDictionary(propsFrom).rename( propsFrom, propsTo); // Subset points that intersect the given image. var fcSub = fc.filterBounds(img.geometry()); // Reduce the image by regions. return img.reduceRegions({ collection: fcSub, reducer: _params.reducer, scale: _params.scale, crs: _params.crs }) // Add metadata to each feature. .map(function(f) { return f.set(imgProps); }); // Converts the feature collection of feature collections to a single //feature collection. }).flatten(); return results; }
代码解释：
这段代码定义了一个名为 zonalStats 的函数，用于对指定的影像集合（ic）内的每个影像在给定的要素集合（fc）内进行区域统计分析。该函数使用了多个参数，允许用户自定义统计分析的具体细节。下面是对这个函数的详细解释：
函数参数
-
ic：要分析的影像集合（ee.ImageCollection）。fc：要在其中进行统计分析的要素集合（ee.FeatureCollection），通常代表一组地理区域。params：一个包含多个可选参数的对象，用于控制分析的细节。
参数字典 _params 的初始化与自定义
-
_params：一个内部参数字典，初始化包含默认的统计分析参数，例如使用ee.Reducer.mean()作为默认简化器（reducer），意味着默认计算区域的平均值。- 如果提供了
params参数，代码会使用params中的值来更新_params中的相应值。
设置默认参数
-
- 通过影像集合的第一个影像（
imgRep）来设置一些默认参数，例如波段名（bands）和非系统影像属性（nonSystemImgProps）。
- 通过影像集合的第一个影像（
影像集合的迭代处理
-
- 对
ic中的每个影像进行迭代，为每个影像执行以下操作：- 选择并可选地重命名波段。
- 设置日期时间和时间戳属性。
- 定义最终输出要素（features）中要设置的影像属性字典。
- 对
减少区域操作
-
- 对每个影像，使用
reduceRegions函数将影像数据按照fc中定义的区域进行聚合。聚合操作根据_params中的简化器、比例尺（scale）和坐标参考系统（CRS）来执行。 - 每个区域（feature）的统计结果被添加上影像的属性。
- 对每个影像，使用
返回结果
-
- 最后，
map函数的调用结果（一个特征集合的集合）被flatten，转换成一个单一的特征集合。这个结果集合包含了原始fc中每个区域针对ic中每个影像计算得到的统计结果。
- 最后，
使用场景
这个 zonalStats 函数可以用于多种场景，如计算给定区域内的平均温度、降水量、植被指数等，或者进行时间序列分析，以监测某个区域随时间的环境变化。通过提供不同的简化器，用户可以计算总和、最大值、最小值等统计量。此外，该函数的灵活性允许用户根据需求调整比例尺、坐标参考系统以及进行统计分析的波段。
三、创建点要素集合
var pts = ee.FeatureCollection([
ee.Feature(ee.Geometry.Point([-118.6010, 37.0777]), { plot_id: 1 }),
ee.Feature(ee.Geometry.Point([-118.5896, 37.0778]), { plot_id: 2 }),
ee.Feature(ee.Geometry.Point([-118.5842, 37.0805]), { plot_id: 3 }),
ee.Feature(ee.Geometry.Point([-118.5994, 37.0936]), { plot_id: 4 }),
ee.Feature(ee.Geometry.Point([-118.5861, 37.0567]), { plot_id: 5 })
]);
print('Points of interest', pts);
Map.addLayer(pts,{plette:'red'});
Map.centerObject(pts,9);
四、邻域统计示例
4.1 地形变量 Topographic Variables
//==========================section 1=============================================
var pts = ee.FeatureCollection([
ee.Feature(ee.Geometry.Point([-118.6010, 37.0777]), { plot_id: 1 }),
ee.Feature(ee.Geometry.Point([-118.5896, 37.0778]), { plot_id: 2 }),
ee.Feature(ee.Geometry.Point([-118.5842, 37.0805]), { plot_id: 3 }),
ee.Feature(ee.Geometry.Point([-118.5994, 37.0936]), { plot_id: 4 }),
ee.Feature(ee.Geometry.Point([-118.5861, 37.0567]), { plot_id: 5 })
]);
print('Points of interest', pts);
//Map.addLayer(pts);
//Map.centerObject(pts,9);
//==========================section 2 定义函数 bufferPoints和zonalStats =============================================
function bufferPoints(radius, bounds) {
return function(pt) {
pt = ee.Feature(pt);
return bounds ? pt.buffer(radius).bounds() :
pt.buffer(radius);
};
}
function zonalStats(ic, fc, params) {
// 初始化内部参数字典。
var _params = {
reducer: ee.Reducer.mean(),
scale: null,
crs: null,
bands: null,
bandsRename: null,
imgProps: null,
imgPropsRename: null,
datetimeName: 'datetime',
datetimeFormat: 'YYYY-MM-dd HH:mm:ss'
};
// 用提供的参数替换初始化的参数。
if (params) {
for (var param in params) {
_params[param] = params[param] || _params[param];
}
}
// 根据图像代表设置默认参数.
var imgRep = ic.first();
var nonSystemImgProps = ee.Feature(null)
.copyProperties(imgRep).propertyNames();
if (!_params.bands) _params.bands = imgRep.bandNames();
if (!_params.bandsRename) _params.bandsRename = _params.bands;
if (!_params.imgProps) _params.imgProps = nonSystemImgProps;
if (!_params.imgPropsRename) _params.imgPropsRename = _params
.imgProps;
// 对影像集合迭代使用reduceRegions函数.
var results = ic.map(function(img) {
// 选择波段 (可选重命名)，设置日期时间和时间戳属性。.
img = ee.Image(img.select(_params.bands, _params
.bandsRename))
// Add datetime and timestamp features.
.set(_params.datetimeName, img.date().format(
_params.datetimeFormat))
.set('timestamp',
img.get('system:time_start'));
// Define final image property dictionary to set in output features.
var propsFrom = ee.List(_params.imgProps)
.cat(ee.List([_params.datetimeName,
'timestamp']));
var propsTo = ee.List(_params.imgPropsRename)
.cat(ee.List([_params.datetimeName,
'timestamp']));
var imgProps = img.toDictionary(propsFrom).rename(
propsFrom, propsTo);
// Subset points that intersect the given image.
var fcSub = fc.filterBounds(img.geometry());
// Reduce the image by regions.
return img.reduceRegions({
collection: fcSub,
reducer: _params.reducer,
scale: _params.scale,
crs: _params.crs
})
// Add metadata to each feature.
.map(function(f) {
return f.set(imgProps);
});
// Converts the feature collection of feature collections to a single
//feature collection.
}).flatten();
return results;
}
//========================== section 3 调用 bufferPoints =====================================
// Buffer the points.
var ptsTopo = pts.map(bufferPoints(45, false));
print('ptsTopo', ptsTopo);
//var ptsTopo1 = pts.map(bufferPoints(45, true));
//print('ptsTopo1', ptsTopo1);
//Map.addLayer(ptsTopo);
//Map.addLayer(ptsTopo1);
//==========================section 4 调用zonalStats=============================================
//!!!zonalStats函数中的ic参数是影像集合，如果只有单张影像的话，要放到影像集合中。
//此外，影像集合中的每一幅影像必须有‘system:time_start’ 属性，一般影像都会有，没有的话，添加上.
// Import the MERIT global elevation dataset.
var elev = ee.Image('MERIT/DEM/v1_0_3');
// 根据DEM计算坡度.
var slope = ee.Terrain.slope(elev);
// Concatenate elevation and slope as two bands of an image.
// Computed images do not have a 'system:time_start' property; add one based on when the data were collected.
var topo = ee.Image.cat(elev, slope)
.set('system:time_start', ee.Date('2000-01-01').millis());
// Wrap the single image in an ImageCollection for use in the zonalStats function.
var topoCol = ee.ImageCollection([topo]);
// Define parameters for the zonalStats function.
var params = {
bands: [0, 1],
bandsRename: ['elevation', 'slope']
};
// Extract zonal statistics per point per image.
var ptsTopoStats = zonalStats(topoCol, ptsTopo, params);
print('Topo zonal stats table', ptsTopoStats);
// ================================section 5 Display the layers on the map. ===============================================
Map.setCenter(-118.5957, 37.0775, 13);
Map.addLayer(topoCol.select(0), {
min: 2400,
max: 4200
}, 'Elevation');
Map.addLayer(topoCol.select(1), {
min: 0,
max: 60
}, 'Slope');
Map.addLayer(pts, {
color: 'purple'
}, 'Points');
Map.addLayer(ptsTopo, {
color: 'yellow'
}, 'Points w/ buffer');
代码解释：
.set('system:time_start', ee.Date('2000-01-01').millis())：这部分代码使用 .set() 方法为上面创建的组合影像设置一个属性，具体是设置 system:time_start 属性的值。ee.Date('2000-01-01').millis() 生成了一个代表2000年1月1日的 ee.Date 对象，并通过 .millis() 方法将这个日期转换为毫秒值（自1970年1月1日以来的毫秒数）。这意味着新的组合影像被赋予了一个时间戳，标记为2000年1月1日。
续：
//===============================第二个例子 MODIS Time Series======================================== var ptsModis = pts.map(bufferPoints(50, true)); var modisCol = ee.ImageCollection('MODIS/006/MOD09A1') .filterDate('2015-01-01', '2020-01-01') .filter(ee.Filter.calendarRange(183, 245, 'DAY_OF_YEAR')); //.filter(ee.Filter.calendarRange(183, 245, 'DAY_OF_YEAR'))： //这个方法进一步过滤影像集合，只保留每年的第183天至第245天之间的影像。 //这个日期范围大约对应每年的7月初到9月初，常用于分析特定季节的地表特征，如植被生长情况。 // Define parameters for the zonalStats function. var params = { reducer: ee.Reducer.median(), scale: 500, crs: 'EPSG:5070', bands: ['sur_refl_b01', 'sur_refl_b02', 'sur_refl_b06'], bandsRename: ['modis_red', 'modis_nir', 'modis_swir'], datetimeName: 'date', datetimeFormat: 'YYYY-MM-dd' }; // Extract zonal statistics per point per image. var ptsModisStats = zonalStats(modisCol, ptsModis, params); print('Limited MODIS zonal stats table', ptsModisStats.limit(50)); //===============================第3个例子 Landsat Time Series ======================================== //--------------- Mask clouds from images and apply scaling factors.-------------------------- function maskScale(img) { var qaMask = img.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); var saturationMask = img.select('QA_RADSAT').eq(0); // Apply the scaling factors to the appropriate bands. var getFactorImg = function(factorNames) { var factorList = img.toDictionary().select(factorNames) .values(); return ee.Image.constant(factorList); }; var scaleImg = getFactorImg(['REFLECTANCE_MULT_BAND_.']); var offsetImg = getFactorImg(['REFLECTANCE_ADD_BAND_.']); var scaled = img.select('SR_B.').multiply(scaleImg).add( offsetImg); // Replace the original bands with the scaled ones and apply the masks. return img.addBands(scaled, null, true) .updateMask(qaMask) .updateMask(saturationMask); } // ------------------------Selects and renames bands of interest for Landsat OLI-------------------. function renameOli(img) { return img.select( ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'], ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2']); } // ----------------------Selects and renames bands of interest for TM/ETM+.--------------------- function renameEtm(img) { return img.select( ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'], ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2']); } // Prepares (cloud masks and renames) OLI images. function prepOli(img) { img = maskScale(img); img = renameOli(img); return img; } // Prepares (cloud masks and renames) TM/ETM+ images. function prepEtm(img) { img = maskScale(img); img = renameEtm(img); return img; } var ptsLandsat = pts.map(bufferPoints(15, true)); var oliCol = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterBounds(ptsLandsat) .map(prepOli); var etmCol = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2') .filterBounds(ptsLandsat) .map(prepEtm); var tmCol = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2') .filterBounds(ptsLandsat) .map(prepEtm); var landsatCol = oliCol.merge(etmCol).merge(tmCol); //---------------------------Calculate Zonal Statistics---------------------------------------- // Define parameters for the zonalStats function. var params = { reducer: ee.Reducer.max(), scale: 30, crs: 'EPSG:5070', bands: ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'], bandsRename: ['ls_blue', 'ls_green', 'ls_red', 'ls_nir', 'ls_swir1', 'ls_swir2' ], imgProps: ['SENSOR_ID', 'SPACECRAFT_ID'], imgPropsRename: ['img_id', 'satellite'], datetimeName: 'date', datetimeFormat: 'YYYY-MM-dd' }; // Extract zonal statistics per point per image. // Filter out observations where image pixels were all masked. var ptsLandsatStats = zonalStats(landsatCol, ptsLandsat, params) .filter(ee.Filter.notNull(params.bandsRename)); print('Limited Landsat zonal stats table', ptsLandsatStats.limit(50)); Export.table.toAsset({ collection: ptsLandsatStats, description: 'EEFA_export_Landsat_to_points', assetId: 'EEFA_export_values_to_points' }); Export.table.toDrive({ collection: ptsLandsatStats, folder: 'EEFA_outputs', // this will create a new folder if it doesn't exist description: 'EEFA_export_values_to_points', fileFormat: 'CSV' });
五、其他
5.1 Weighted Versus Unweighted Region Reduction
5.2 Copy Properties to Computed Images
copyProperties（）
5.3 Understanding Which Pixels Are Included in Polygon Statistics
分类:
2-GEE
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 全程不用写代码，我用AI程序员写了一个飞机大战
· MongoDB 8.0这个新功能碉堡了，比商业数据库还牛
· 记一次.NET内存居高不下排查解决与启示
· DeepSeek 开源周回顾「GitHub 热点速览」
· 白话解读 Dapr 1.15：你的「微服务管家」又秀新绝活了