| GEE：对Sentinel | 您所在的位置：网站首页 › 直方图计算最大值与最小值的 › GEE：对Sentinel |
GEE：对Sentinel
|
作者：CSDN @ _养乐多_ 本文介绍了通过Google Earth Engine平台，并使用哨兵数据提取水体掩膜的方法和代码。通过裁剪和去除云等处理步骤，最终得到具有水体掩膜的影像，并进行可视化和导出。这种方法基于归一化水体指数（NDWI）和OTSU阈值计算技术，无需复杂的图像处理算法，适用于快速获取水体信息的需求。  图1 彩色合成
 图2 水体掩膜
文章目录
一、OTSU阈值方法原理二、代码详解三、代码链接四、OTSU算法详解
一、OTSU阈值方法原理
OTSU阈值方法是一种常用的图像分割算法，用于将图像分成背景和目标两个部分。在提取水体的应用中，可以使用OTSU阈值方法来自动分割水体和非水体区域。 OTSU阈值方法基于图像的灰度直方图，通过寻找最佳的阈值来实现分割。以下是OTSU阈值方法的原理步骤： 计算图像的灰度直方图：首先，将图像转换为灰度图像，并统计每个灰度级别的像素数量。这样可以得到图像的灰度直方图。 计算像素的总数：统计图像中的总像素数。 计算类间方差：对于每个可能的阈值T，将图像分成两个类别：低于阈值T的像素类别为背景，高于阈值T的像素类别为目标（水体）。然后，计算两个类别的权重、均值和方差，并使用这些值来计算类间方差。类间方差是一个衡量分割质量的指标，当类间方差最大时，分割效果最佳。 寻找最大类间方差：通过遍历所有可能的阈值T，计算类间方差，并找到使类间方差最大的阈值T*。 应用阈值分割：将图像的每个像素与阈值T*进行比较，将像素分为水体和非水体两个部分。 OTSU阈值方法的原理是基于类间方差最大化的思想。通过最大化类间方差，可以实现最佳的分割效果，将水体和非水体区域有效地分开。这种方法不需要预先设定阈值，而是通过自动计算得到最优的阈值，因此具有较好的适应性和稳定性。 在应用中，可以将OTSU阈值方法应用于水体提取任务，通过对输入图像进行OTSU阈值分割，得到水体和非水体区域的二值图像，从而实现水体的提取和分割。 二、代码详解 var roi = table Map.centerObject(roi,10) // 可视化参数 var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01', '011D01', '011301']; //可视化参数，按843波段合成 var rgbVis = { min: 0.0, max: 0.35, bands: ['B8', 'B4', 'B3'], }; //按矢量边界裁剪 function roiClip(image){ return image.clip(roi) } //S2_SR去云 function remove_cloud(image) { var qa = image.select('QA60') var cloudBitMask = 1 var counts = ee.Array(ee.Dictionary(histogram).get('histogram')); var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans')); var size = means.length().get([0]); var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]); var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]); var mean = sum.divide(total); var indices = ee.List.sequence(1, size); var bss = indices.map(function(i) { var aCounts = counts.slice(0, 0, i); var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]); var aMeans = means.slice(0, 0, i); var aMean = aMeans.multiply(aCounts) .reduce(ee.Reducer.sum(), [0]).get([0]) .divide(aCount); var bCount = total.subtract(aCount); var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount); return aCount.multiply(aMean.subtract(mean).pow(2)).add( bCount.multiply(bMean.subtract(mean).pow(2))); }); return means.sort(bss).get([-1]); } //计算水体掩膜 function calWaterMask(image){ var MNDWI_S2 = image.select("NDWI"); var histogram = MNDWI_S2.reduceRegion({ reducer: ee.Reducer.histogram(), geometry: roi, scale: 30, maxPixels: 1e13, tileScale: 8 }); var threshold_S2 = otsu(histogram.get("NDWI")); var mask = MNDWI_S2.gte(threshold_S2); return image.addBands(mask.rename('mask')); } //运行水体掩模函数 newImgS2 = calNDWI(newImgS2) newImgS2 = calWaterMask(newImgS2).select('mask') print(newImgS2) Map.addLayer(newImgS2, {min:0, max:1, palette:['#DDDDDD', '#0099FF']}, 'S2_Water_body_mask'); //下载非水体概率影像 Export.image.toDrive({ image:newImgS2, description: 'S2_water_mask', scale:30, region:roi, fileFormat: 'GeoTIFF', maxPixels:1e13, });这段代码的目的是对Sentinel-2遥感影像进行处理，提取水体掩膜，并将结果可视化和导出。 首先，定义了一个变量roi，表示感兴趣区域（Region of Interest），这是一个表格（table）对象。然后使用Map.centerObject(roi,10)将该区域设置为地图的中心，并指定缩放级别为10。 接下来，定义了一个调色板（palette）数组，用于可视化图像时的颜色映射。 然后，定义了一个可视化参数rgbVis，设置了最小值（min）和最大值（max），以及要显示的波段（bands）。 之后，定义了一个函数roiClip，用于将图像按照矢量边界进行裁剪。 接着是一个函数remove_cloud，用于去除图像中的云和雾。函数内部使用了位掩码（bitmask）的方式来选择需要保留的像素，然后将图像的像素值除以10000并选择特定波段，最后将处理后的图像的时间属性复制到结果中。 然后定义了起始日期和结束日期，用于选择要处理的影像时间范围。 接下来，使用ee.ImageCollection函数选择了名为"COPERNICUS/S2_SR"的Sentinel-2影像集合，并通过一系列操作对影像进行筛选、裁剪和去除云，最后将结果映射到ROI区域，并存储在变量S2中。 之后，使用S2.max()函数选取ROI区域内同一天的影像，并将结果存储在变量newImgS2中。 接着，定义了一个名为calNDWI的函数，用于计算归一化水体指数（NDWI）。 然后，定义了一个名为otsu的函数，用于计算阈值，该阈值用于将图像的像素分为水体和非水体。 接下来，定义了一个名为calWaterMask的函数，该函数通过计算NDWI的直方图，并使用OTSU方法计算阈值来生成水体掩膜。 然后，分别调用calNDWI和calWaterMask函数对newImgS2进行处理，并将结果可视化。 最后，使用Export.image.toDrive函数将非水体概率影像导出到Google Drive中。 三、代码链接https://code.earthengine.google.com/f145e92abf2d0196aba61655745bad9c?noload=true 四、OTSU算法详解 // OTSU计算阈值 // 此函数接受一个直方图作为输入，并使用OTSU算法计算出最佳阈值作为图像分割的阈值。 function otsu(histogram) { // 从直方图中获取像素计数和桶均值数组 var counts = ee.Array(ee.Dictionary(histogram).get('histogram')); var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans')); // 获取桶均值数组的长度，也就是直方图的桶数目 var size = means.length().get([0]); // 计算像素总数 var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]); // 计算桶均值乘以像素计数的总和 var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]); // 计算全局均值 var mean = sum.divide(total); // 创建从1到size的索引列表 var indices = ee.List.sequence(1, size); // 计算类间方差（between-class variance） var bss = indices.map(function(i) { // 将counts数组切片，获取前i个桶的像素计数 var aCounts = counts.slice(0, 0, i); // 计算前i个桶的像素计数之和 var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]); // 将means数组切片，获取前i个桶的桶均值 var aMeans = means.slice(0, 0, i); // 计算前i个桶的桶均值乘以像素计数的总和，然后除以像素计数之和 var aMean = aMeans.multiply(aCounts) .reduce(ee.Reducer.sum(), [0]).get([0]) .divide(aCount); // 计算后i个桶的像素计数之和 var bCount = total.subtract(aCount); // 计算后i个桶的桶均值，公式为：(总和 - 前i个桶的像素计数之和 * 前i个桶的桶均值) / 后i个桶的像素计数之和 var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount); // 计算类间方差 return aCount.multiply(aMean.subtract(mean).pow(2)).add( bCount.multiply(bMean.subtract(mean).pow(2))); }); // 将类间方差数组按照从小到大排序，获取最后一个（最大值），即最佳阈值 return means.sort(bss).get([-1]); }OTSU算法是一种自适应的图像阈值分割方法，它通过最大化类间方差来选择最佳阈值。 具体步骤如下： 计算图像的直方图，即每个像素值的像素计数。计算总像素数和每个像素值的占比。对每个像素值，计算类内方差，即以该像素值作为阈值进行分割时，前景和背景之间的方差。选择使类内方差最小的像素值作为阈值。利用该阈值将图像分割为前景和背景两部分。这段代码的作用是实现OTSU算法的第4步，即计算类间方差并选择最佳阈值。 |
【本文地址】