当前位置：首页 > news >正文

NDVI计算

news 2026/5/16 21:33:01

1.研究区

// 1. 加载研究区并可视化 var roi = table; Map.centerObject(roi, 8); Map.addLayer(roi, {'color': 'yellow'}, '研究区'); // 2. 优化云掩膜函数 function maskS2clouds(image) { var roiMask = ee.Image.constant(1).clip(roi).mask(); // 细分云掩码 var opaqueCloudMask = image.select('MSK_CLASSI_OPAQUE').lte(0); var cirrusMask = image.select('MSK_CLASSI_CIRRUS').lte(0); var snowIceMask = image.select('MSK_CLASSI_SNOW_ICE').eq(0); var cloudMask = opaqueCloudMask.and(cirrusMask).and(snowIceMask); // SCL分类掩码 var scl = image.select('SCL'); var sclMask = scl.neq(1).and(scl.neq(3)) .and(scl.neq(7)).and(scl.neq(8)).and(scl.neq(9) .and(scl.neq(10)).and(scl.neq(11))); var finalMask = cloudMask.and(sclMask); return image.updateMask(finalMask) .divide(10000) .copyProperties(image, ['system:time_start']) .set('mask', finalMask); } // 3. 加载固定时间段影像集合 var Sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") .filterDate('2025-07-01', '2025-08-31') .filterBounds(roi) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) .map(maskS2clouds) .select(['B2', 'B3', 'B4', 'B8']); print('可用影像数量', Sentinel2.size()); if (Sentinel2.size().getInfo() < 1) { print('警告：没有可用影像，请调整云量阈值或时间范围'); } // 4. 直接生成单幅合成影像（不进行填充） // 使用中值合成，对云敏感区域效果较好 var medianImage = Sentinel2.median().clip(roi); // 可视化单幅影像 Map.addLayer(Sentinel2, {bands: ['B4','B3','B2'], min:0, max:0.3}, '原始影像集合'); Map.addLayer(medianImage, {bands: ['B4','B3','B2'], min:0, max:0.3}, '中值合成影像'); // 5. 导出单幅影像数据至Drive Export.image.toDrive({ image: medianImage.select(['B4', 'B8']), description: 'shiyan_NDVI', folder: 'shiyan_NDVI', fileNamePrefix: 'median_image', region: roi, scale: 10, maxPixels: 1e13, fileFormat: 'GeoTIFF', crs: 'EPSG:4326' });

2.NDVI

2.1 直接计算NDVI

由于研究时段仅为2025年的7月至8月，加上云掩膜函数，导致原始的遥感影像存在一些空缺区域，直接在原始单幅影像上计算NDVI，NDVI影像也会继承这些空缺。

// 1. 加载研究区并可视化 var roi = table; Map.centerObject(roi, 8); Map.addLayer(roi, {'color': 'yellow'}, '研究区'); // 2. 优化云掩膜函数 function maskS2clouds(image) { var roiMask = ee.Image.constant(1).clip(roi).mask(); // 细分云掩码 var opaqueCloudMask = image.select('MSK_CLASSI_OPAQUE').lte(0); var cirrusMask = image.select('MSK_CLASSI_CIRRUS').lte(0); var snowIceMask = image.select('MSK_CLASSI_SNOW_ICE').eq(0); var cloudMask = opaqueCloudMask.and(cirrusMask).and(snowIceMask); // SCL分类掩码 var scl = image.select('SCL'); var sclMask = scl.neq(1).and(scl.neq(3)) .and(scl.neq(7)).and(scl.neq(8)).and(scl.neq(9) .and(scl.neq(10)).and(scl.neq(11))); var finalMask = cloudMask.and(sclMask); return image.updateMask(finalMask) .divide(10000) .copyProperties(image, ['system:time_start']) .set('mask', finalMask); } // 3. 加载固定时间段影像集合 var Sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") .filterDate('2025-07-01', '2025-08-31') .filterBounds(roi) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) .map(maskS2clouds) .select(['B2', 'B3', 'B4', 'B8']); print('可用影像数量', Sentinel2.size()); if (Sentinel2.size().getInfo() < 1) { print('警告：没有可用影像，请调整云量阈值或时间范围'); } // 4. 生成单幅合成影像 var medianImage = Sentinel2.median().clip(roi); // 计算NDVI值（仅保留内置方法） var ndvi = medianImage.normalizedDifference(['B8', 'B4']).rename('NDVI'); // NDVI可视化参数设置（优化冗余颜色） var ndviVis = { min: -1, max: 1, palette: ['blue', 'green'] // 颜色说明： // 空白 - 空值区域（保持地图底图显示） // 从blue到green渐变： // NDVI ≤ 0 → 蓝色（水体/裸地等） // NDVI > 0 → 绿色（植被区域） }; // 添加NDVI图层到地图 Map.addLayer(ndvi, ndviVis, 'NDVI'); // 导出NDVI数据至Drive Export.image.toDrive({ image: ndvi, description: 'shiyan_NDVI', folder: 'shiyan_NDVI', fileNamePrefix: 'ndvi_image', region: roi, scale: 10, maxPixels: 1e13, fileFormat: 'GeoTIFF', crs: 'EPSG:4326' });

2.2 NDVI插补

改进了原始遥感影像的获取方法，先获取25年6至9月的影像，以7、8月份的影像为主，对其进行云掩膜和中值合成，空缺区域由6、9月份的影像进行插补。这样就生成了一张无空缺的遥感影像图。

再在无空缺的遥感影像图的基础上运用公式计算NDVI，生成NDVI影像图。

/** * Sentinel-2 NDVI计算与插值处理代码（定向导出版） * 功能：生成6月、7-8月、9月三幅影像，仅导出插值后的7-8月NDVI影像 * 版本：1.9（仅导出插值后的7-8月结果） * 最后更新：2025-09-30 */ // ==================== 1. 研究区设置 ==================== var roi = table; // 研究区矢量数据（需在GEE中预先定义） Map.centerObject(roi, 7); // 地图中心定位到研究区，缩放级别7 // 可视化研究区边界 var styling = {color: "red", fillColor: "00000000"}; // 红色边界，透明填充 Map.addLayer(roi.style(styling), {}, "研究区边界"); // ==================== 2. 辅助函数定义 ==================== /** * 影像归一化函数 * 将影像波段值归一化到[0,1]范围 * @param {ee.Image} img - 输入影像 * @return {ee.Image} 归一化后的影像 */ var img_normalize = function(img) { // 计算研究区内各波段的最大最小值 var minMax = img.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 10, // 采用10m分辨率（Sentinel-2原生分辨率） maxPixels: 10e13, // 支持大区域计算 tileScale: 16 // 提高计算并行度，避免内存超限 }); var year = img.get('year'); // 获取影像年份属性 // 对每个波段进行归一化处理 var normalize = ee.ImageCollection.fromImages( img.bandNames().map(function(name) { name = ee.String(name); var band = img.select(name); // 线性归一化到[0,1] return band.unitScale( ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))) ); }) ).toBands().rename(img.bandNames()); return normalize; }; /** * Sentinel-2云掩膜函数（基于SCL波段） * 替代已弃用的QA60波段，保留特定地物类型 * SCL分类说明： * 0 - 无数据 | 1 - 饱和或缺陷像素 | 2 - 暗部像素 * 3 - 云影 | 4 - 植被 | 5 - 非植被 | 6 - 水 * 7 - 未分类 | 8 - 云（低概率） | 9 - 云（中概率） * 10 - 云（高概率） | 11 - 薄卷云 * @param {ee.Image} image - 输入Sentinel-2影像 * @return {ee.Image} 掩膜处理后的影像 */ function maskS2clouds(image) { // 选择SCL波段（场景分类图） var scl = image.select('SCL'); // 定义保留的地物类型：2（暗部像素）、4（植被）、5（非植被）、6（水） // 排除云、云影、卷云、无数据等干扰类型 var mask = scl.eq(2) .or(scl.eq(4)) .or(scl.eq(5)) .or(scl.eq(6)); // 应用掩码并缩放反射率值（Sentinel-2数据原始值需除以10000） return image.updateMask(mask) .divide(10000) .copyProperties(image, ['system:time_start']); // 保留时间属性 } // ==================== 3. 数据集加载与处理 ==================== /** * 加载Sentinel-2数据集 * 数据源：COPERNICUS/S2_SR_HARMONIZED（经过大气校正的地表反射率数据） * 包含SCL波段用于场景分类 */ var imageCollection = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") .filterBounds(roi) .select(['B4', 'B8', 'SCL']); // 仅选择计算NDVI所需波段和SCL波段 // 生成目标月份序列（共3组：6月、7-8月、9月） var targetMonths = ee.List([ {start: 6, end: 6, name: "6月"}, // 6月（插值数据源） {start: 7, end: 8, name: "7-8月"}, // 7-8月（主要研究期） {start: 9, end: 9, name: "9月"} // 9月（插值数据源） ]); // ==================== 4. 目标月份NDVI合成 ==================== /** * 生成目标月份NDVI中值合成影像 * 对每个目标月份/月份组，计算中值，提高数据稳定性 */ var composites = ee.ImageCollection.fromImages(targetMonths.map(function(group) { group = ee.Dictionary(group); var startMonth = group.getNumber('start'); var endMonth = group.getNumber('end'); var groupName = group.getString('name'); var startYear = ee.Number(2025); // 计算该月份组的代表性日期（用于时间属性） var representativeDate = ee.Date.fromYMD(startYear, startMonth, 15) .advance(ee.Number(endMonth).subtract(startMonth), 'month'); // 筛选过去2年同期的影像（提高数据可用性） var filtered = imageCollection.filter(ee.Filter.calendarRange({ start: startYear, // 起始年份：当前年份 end: startYear, // 结束年份：当前年份 field: 'year' })).filter(ee.Filter.calendarRange({ start: startMonth, // 筛选起始月份 end: endMonth, // 筛选结束月份 field: 'month' })) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)); // 筛选云量<20%的影像 // 云掩膜处理并计算中值合成 var composite = filtered.map(maskS2clouds).median().clip(roi); // 计算NDVI（Sentinel-2：近红外B8 - 红波段B4） return composite.normalizedDifference(['B8', 'B4']).rename('NDVI') .set('period', groupName) .set('system:time_start', representativeDate.millis()) // 标记主要研究期（7-8月） .set('is_primary_period', ee.String(groupName).equals("7-8月")); })); print("目标月份NDVI合成影像集合", composites); // ==================== 5. 影像集合堆叠函数 ==================== /** * 将影像集合堆叠为多波段影像 * 便于后续插值处理和结果导出 */ var stackCollection = function(collection) { // 创建空影像作为初始值 var first = ee.Image(collection.first()).select([]); // 定义波段追加函数 var appendBands = function(image, previous) { return ee.Image(previous).addBands(image); }; // 迭代追加所有波段 return ee.Image(collection.iterate(appendBands, first)); }; // 堆叠插值前的目标月份影像 var compos = stackCollection(composites); print('插值前的多波段影像', compos); // ==================== 6. 插值处理 ==================== /** * 对7-8月影像使用6月和9月影像进行插值 * 使用修复后的ImageCollection.mean()方法 */ var replacedVals = composites.map(function(image) { var currentPeriod = image.get('period'); var currentDate = ee.Date(image.get('system:time_start')); // 仅对7-8月影像进行插值处理 return ee.Algorithms.If( ee.String(currentPeriod).equals("7-8月"), // 1. 获取6月和9月的影像，转换为ImageCollection再计算平均值 ee.ImageCollection.fromImages([ composites.filter(ee.Filter.eq('period', '6月')).first(), composites.filter(ee.Filter.eq('period', '9月')).first() ]).mean() // 2. 仅保留插值数据源中原始影像为空的区域 .updateMask(image.mask().not()) // 3. 将原始影像和插值结果合并 .addBands(image) // 4. 关键步骤：优先使用原始影像的值，空值区域使用插值结果 .reduce(ee.Reducer.firstNonNull()) .rename('NDVI') .set('period', currentPeriod) .set('system:time_start', currentDate.millis()), // 其他月份影像：直接返回原始影像 image.set('period', currentPeriod) ); }); // 将结果转换为ImageCollection var replacedValsCol = ee.ImageCollection(replacedVals); // 堆叠插值后的目标月份影像 var stacked = stackCollection(replacedValsCol); print('插值后的多波段影像', stacked); // ==================== 7. 结果可视化 ==================== // 简化的NDVI可视化参数：小于0为蓝色，大于0为绿色 var Vis = { min: -1, max: 1.0, palette: [ 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01', '011D01', '011301' ] }; // 仅显示7-8月插值前后的影像（索引1对应7-8月） Map.addLayer(compos.select(1), Vis, '插值前（7-8月）'); Map.addLayer(stacked.select(1), Vis, 'NDVI（7-8月）'); // ==================== 8. 结果导出（仅导出插值后的7-8月影像） ==================== // 找到7-8月在目标月份列表中的索引 var julyAugustIndex = targetMonths.map(function(group) { return ee.String(ee.Dictionary(group).get('name')).equals("7-8月"); }).indexOf(true).getInfo(); // 仅导出插值后的7-8月NDVI影像 var exportBand = stacked.bandNames().get(julyAugustIndex); exportBand.evaluate(function(band) { var periodName = "7至8月"; Export.image.toDrive({ image: stacked.select(band), description:'shiyan_NDVI', folder: 'shiyan_NDVI', fileNamePrefix: 'NDVI_image', crs: "EPSG:4326", scale: 10, region: roi, maxPixels: 1e13, fileFormat: 'GeoTIFF', }); });

2.3 NDVI栅格

gee平台对单次上传的文件有大小限制，若生成的栅格较大，gee会将其分割后再进行导出。

需要下载分割后的栅格，在arcgispro中使用“镶嵌”工具将其合并为完整影像，再使用“按掩膜提取”工具进行提取，从而得到所需要的研究区的影像。

gee导出的栅格的空间参考是WGS 1984地理坐标系，需要使用“投影栅格”工具将其转化为“WGS 1984 UTM Zone 49N”投影坐标系。

3.FVC

3.1 NDVI分位数

NDVI值用于量化地表植被的覆盖情况和生长状态。其值域为-1到1之间，具有明确的物理意义：

-1：地面覆盖为云、水、雪等，这些地表类型对可见光波段具有高反射性。

0：地表有岩石或裸土等，此时近红外波段（NIR）和红光波段（R）的反射率大致相等。

正值：有植被覆盖，且NDVI值随植被覆盖度的增大而增大。

计算FVC需要NDVImin和NDVImax（即NDVIsoil和 NDVIveg）

NDVIsoil：裸地或无植被覆盖的像元NDVI值（5%）；

NDVIveg：全部被植被覆盖的像元NDVI值（95%）。

查看方法如下：

在栅格图层的“符号系统”窗格中，将绘制类型从“拉伸”更改为“分类”。在方法中选择“分位数”，并将类的数量设置为 20。应用设置后，ArcGIS Pro 会将NDVI值分成20组，每组包含大约5%的像元。

列表中第1个类别的上限值就非常接近您要找的 5%分位数（NDVIsoil）。因为它表示有5%的像元NDVI值低于或等于这个数值。
相应地，第19个类别的上限值则非常接近 95%分位数（NDVIveg），表示有95%的像元NDVI值低于或等于这个数

结果：

NDVImin = NDVIsoil = 0.364706

NDVImax = NDVIveg = 0.921569

3.2 FVC计算

// ==================== 1. 研究区设置 ==================== var roi = table; Map.centerObject(roi, 7); var styling = {color: "red", fillColor: "00000000"}; Map.addLayer(roi.style(styling), {}, "研究区边界"); // ==================== 2. 辅助函数定义 ==================== // 影像归一化函数（保持10m分辨率） var img_normalize = function(img) { var minMax = img.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 10, // 保持原始分辨率 maxPixels: 1e13, tileScale: 16, bestEffort: true // 允许在必要时微调 }); var year = img.get('year'); var normalize = ee.ImageCollection.fromImages( img.bandNames().map(function(name) { name = ee.String(name); var band = img.select(name); return band.unitScale( ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))) ); }) ).toBands().rename(img.bandNames()); return normalize; }; // 云掩膜函数 function maskS2clouds(image) { var scl = image.select('SCL'); var mask = scl.eq(2).or(scl.eq(4)).or(scl.eq(5)).or(scl.eq(6)); return image.updateMask(mask).divide(10000) .copyProperties(image, ['system:time_start']); } // ==================== 3. 数据集加载与处理 ==================== var imageCollection = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") .filterBounds(roi) .select(['B4', 'B8', 'SCL']); var targetMonths = ee.List([ {start: 6, end: 6, name: "6月"}, {start: 7, end: 8, name: "7-8月"}, {start: 9, end: 9, name: "9月"} ]); // ==================== 4. 目标月份NDVI合成 ==================== var composites = ee.ImageCollection.fromImages(targetMonths.map(function(group) { group = ee.Dictionary(group); var startMonth = group.getNumber('start'); var endMonth = group.getNumber('end'); var groupName = group.getString('name'); var startYear = ee.Number(2025); var representativeDate = ee.Date.fromYMD(startYear, startMonth, 15) .advance(ee.Number(endMonth).subtract(startMonth), 'month'); var filtered = imageCollection.filter(ee.Filter.calendarRange({ start: startYear.subtract(1), end: startYear, field: 'year' })).filter(ee.Filter.calendarRange({ start: startMonth, end: endMonth, field: 'month' })) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)); var composite = filtered.map(maskS2clouds).median().clip(roi); return composite.normalizedDifference(['B8', 'B4']).rename('NDVI') .set('period', groupName) .set('system:time_start', representativeDate.millis()) .set('is_primary_period', ee.String(groupName).equals("7-8月")); })); // ==================== 5. 影像集合堆叠函数 ==================== var stackCollection = function(collection) { var first = ee.Image(collection.first()).select([]); var appendBands = function(image, previous) { return ee.Image(previous).addBands(image); }; return ee.Image(collection.iterate(appendBands, first)); }; var compos = stackCollection(composites); // ==================== 6. 插值处理 ==================== var replacedVals = composites.map(function(image) { var currentPeriod = image.get('period'); var currentDate = ee.Date(image.get('system:time_start')); return ee.Algorithms.If( ee.String(currentPeriod).equals("7-8月"), ee.ImageCollection.fromImages([ composites.filter(ee.Filter.eq('period', '6月')).first(), composites.filter(ee.Filter.eq('period', '9月')).first() ]).mean() .updateMask(image.mask().not()) .addBands(image) .reduce(ee.Reducer.firstNonNull()) .rename('NDVI') .set('period', currentPeriod) .set('system:time_start', currentDate.millis()), image.set('period', currentPeriod) ); }); var replacedValsCol = ee.ImageCollection(replacedVals); var stacked = stackCollection(replacedValsCol); //print('插值后的多波段影像', stacked); // ==================== 8. NDVI结果可视化 ==================== var Vis = { min: -1, max: 1.0, palette: [ 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01', '011D01', '011301' ] }; // 从堆叠影像中提取插值后的7-8月NDVI（索引为1） var NDVI_2025 = stacked.select(1).rename('B1'); Map.addLayer(NDVI_2025, Vis, 'NDVI_2025'); // ==================== 9. FVC计算及可视化 ==================== // 已知的端元值 var NDVIsoil = ee.Number(0.443137); // 5%分位数（裸土NDVI） var NDVIveg = ee.Number(0.921569); // 95%分位数（全植被NDVI） // 获取NDVI影像（波段名为B1） var ndvi = NDVI_2025.select('B1').clip(roi); // 复现ENVI的FVC计算逻辑 var fvc = ee.Image(0) .where(ndvi.lt(NDVIsoil), 0) // 处理纯土壤像元：NDVI < NDVIsoil 时，FVC=0 .where(ndvi.gt(NDVIveg), 1) // 处理纯植被像元：NDVI > NDVIveg 时，FVC=1 .where(ndvi.gte(NDVIsoil).and(ndvi.lte(NDVIveg)), // 处理混合像元：介于两者之间时使用像元二分模型 ndvi.subtract(NDVIsoil).divide(NDVIveg.subtract(NDVIsoil))) .rename('FVC') .clip(roi); // 可视化FVC结果 var fvcVis = { min: 0, max: 1, palette: ['white', 'lightgreen', 'green', 'darkgreen'] }; Map.addLayer(fvc, fvcVis, 'FVC_image'); // 导出FVC影像到Google Drive Export.image.toDrive({ image: fvc, // 要导出的影像 description: 'shiyan_FVC', // 导出任务的描述（会成为文件名的一部分） folder: 'shiyan_NDVI', // Google Drive中的目标文件夹（可自定义） fileNamePrefix: 'FVC_image', // 导出文件的前缀名 region: roi, // 导出范围（研究区） scale: 10, // 分辨率（与原始NDVI保持一致，10米） maxPixels: 1e13, // 允许导出的最大像元数 fileFormat: 'GeoTIFF', // 导出格式（GeoTIFF支持空间信息） crs: 'EPSG:4326' // 坐标参考系（WGS84，可选，默认也是这个） });