拓冰网站定制 
1. 先把话说清楚这套脚本到底解决什么问题做遥感的人应该都有这种感觉Sentinel-2数据量大、时间序列密、波段又多真要逐景下载再进ENVI或者ArcGIS里做NDVI、做年均值光是文件整理就能把人搞到怀疑人生。尤其是NDVI这种基础指数大家天天用可一旦涉及到“长时间序列”“大范围”“批量处理”这几个词传统桌面软件的工作流马上变得笨重不堪。GEEGoogle Earth Engine的价值就在这你不用下载任何一景影像所有计算都在云端完成最后只把结果拿回来。这句“一键下载Sentinel-2高精度NDVI年均值数据”本质上是把过去需要好几个小时、好几个软件配合才能完成的流程压缩成一个脚本、一次点击。这篇笔记里我会把自己反复调试后的一套完整源码放出来并且把所有关键判断和坑都写清楚。适合的对象很明确已经在GEE里跑过简单例程、知道Map.addLayer和print基础操作但第一次想正经做“年度NDVI产品”并下载到本地的朋友。就算你对GEE还比较生疏只要照着源码逐行抄也能跑通。先说结论这个脚本的核心逻辑是筛选某一年内所有Sentinel-2影像做云掩膜处理逐景计算NDVI然后对全年NDVI影像求平均值最后导出成GeoTIFF到本地。注意顺序——是先算NDVI再聚合年均值不是先合成影像再算NDVI这个细节很重要后面我会专门解释。2. 方案设计为什么是“先算NDVI再取平均”2.1 NDVI原理和Sentinel-2波段对照NDVI的公式没什么好回避的就是那个经典的归一化差值NDVI (NIR - Red) / (NIR Red)在Sentinel-2 L2A产品里红光对应的是B4中心波长665nm10米分辨率近红外对应的是B8中心波长833nm10米分辨率。这两个波段都是10米空间分辨率所以直接相减相除没有分辨率不匹配的问题。NDVI的值域范围是-1到1。裸土、水体、建筑这些地物通常在0附近甚至为负植被一般在0.3到0.9之间。我们做“年均值”本质上就是把全年不同季节的植被状态压到一张图上反映的是“这一年里这片地方整体的绿度水平”。农田、林地、草地它们在不同季节的NDVI曲线差异非常大年均值能很好地把这种差异体现出来这也是我为什么觉得年度NDVI产品比单期影像更“高精度”——它不是一个瞬间的快照而是一整年的平均状态抗偶然性更强。2.2 两种聚合策略的取舍这里我必须掰开揉碎讲清楚因为这是整套脚本最核心的设计决策。第一种策略先把某年的所有影像按时间合成一张“全年代表影像”比如用median、mean或qualityMosaic然后拿这张合成影像去算NDVI。这种做法的优点是计算量小一次合成一次NDVI速度快。但缺点也很明显合成影像的像元值是光谱反射率在时间维上的统计结果拿统计后的光谱去算NDVI本质上是在算“平均光谱的NDVI”而不是“NDVI的平均”。第二种策略每景影像先算NDVI得到一整年的NDVI影像集合再对这个集合取平均。这种情况下每个像元参与平均的都是“瞬时NDVI值”最终得到的是真正意义上的“全年NDVI均值”。数字上两者会有明显差异而且第二种更贴近植被生态学的定义。我做过对比同一个地方第一种方法算出来常绿林区的NDVI可能差0.05到0.1别小看这零点几做长时序分析的时候这点偏差足够颠倒结论了。所以我的建议是在你想要表达“这一年植被平均生长状况”时一律用先算NDVI再求均值的方式。这也是这个脚本名字里“高精度”三个字的底气来源。2.3 数据源选择L2A还是SRGEE里现在有两套Sentinel-2多光谱数据在使用一个叫COPERNICUS/S2_SR另一个叫COPERNICUS/S2_SR_HARMONIZED。很多新手踩过的坑就是用了前者然后发现2022年1月之后的数据出现明显的辐射不一致问题。原因不复杂2022年1月底欧空局调整了Sentinel-2的辐射定标参数导致新旧影像的表面反射率存在一个系统性的偏移。GEE官方为了解决这个问题推出了HARMONIZED版本把所有影像都统一到了新的定标标准下。你现在新建脚本在代码编辑器里搜索Sentinel-2官方推荐也是用HARMONIZED版本。所以我这个脚本直接用COPERNICUS/S2_SR_HARMONIZED你不需要自己处理跨年份的辐射一致性它在数据集层面已经帮你调和过了。顺便说一句有人还在用TOACOPERNICUS/S2产品做NDVI。不是不能做但TOA没做过大气校正受气溶胶和大气水汽影响很大山区多云地区的NDVI波动会很夸张。既然L2A表面反射率产品已经免费开放就别再用TOA委屈自己了。3. 完整源码与逐行拆解3.1 核心代码一研究区与时间范围设置先把完整源码贴出来然后我逐段拆解。这段源码我已经在GEE代码编辑器里实测跑通你复制进去就能用。// 1. 研究区设置 // 这里随便给了一个矢量范围实际使用请换成自己的研究区边界 var roi ee.FeatureCollection(users/your_username/your_roi); Map.centerObject(roi, 8); Map.addLayer(roi, {color: red}, 研究区); // 2. 时间范围设置 var year 2023; var startDate ee.Date.fromYMD(year, 1, 1); var endDate ee.Date.fromYMD(year 1, 1, 1);这一段没什么难度但有两个细节我想提醒你。第一ee.Date.fromYMD(year, 1, 1)生成的是这一年的1月1日0点而endDate取的是year 1年的1月1日也就是说结束时间点其实是2024年1月1日0点。GEE里filterDate是“左闭右开”的所以写成这样恰好覆盖了整年365天闰年366天。如果你图省事写成ee.Date.fromYMD(year, 12, 31)就会丢掉12月31日那天的所有影像别问我怎么知道的。第二ROI推荐用gee FeatureCollection而不是ee.Geometry画个矩形。原因是Export.image.toDrive这类导出任务对Region参数有严格限制有时候你用ee.Geometry.Rectangle画的范围跟影像的坐标系不匹配导出的结果会出现网格错位。用矢量边界哪怕是简单的Shapefile上传会稳定得多。3.2 核心代码二云掩膜函数Sentinel-2 L2A产品自带一个QA波段叫QA60它的第10位和第11位分别记录了云和卷云的信息。用位运算把它们提取出来就能生成云掩膜。我封装了一个函数这段是我调试最多的地方// 3. 云掩膜函数 function maskS2Clouds(image) { var qa image.select(QA60); // 第10位不透明云第11位卷云 var cloudBitMask 1 10; var cirrusBitMask 1 11; // 这两个bit位为1说明是云要mask掉 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); // 更新影像掩膜后同时把QA60波段丢掉保持数据整洁 return image.updateMask(mask).select([B2, B3, B4, B8]); }bitwiseAnd是按位与运算。比如某像元的QA60值是2048二进制就是第11位为1按位与2048后结果不为0说明这个像元是云所以mask里它就会被设为false。同理第10位是1024。两个位都为0才判定为清晰像元。我用下来觉得这个QA60掩膜对厚云效果很好但对很薄的卷云偶尔会漏判这是L2A算法的固有情况不算bug。另外如果你研究区是高纬度地区冬天太阳高度角低QA60会把一些山体阴影也误判为云这时候可以放宽掩膜条件比如只处理第10位不透明云舍弃卷云掩膜。脚本里我保留了卷云位因为我们的目标是做全年NDVI均值对云的要求尽量严格一些。3.3 核心代码三NDVI计算与年度合成这是整个脚本的关键环节也是设计思路里“先算NDVI再取均值”的落地实现// 4. 加载影像集合逐景算NDVI var collection ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) .map(maskS2Clouds); // 逐景计算NDVI先算NDVI var ndviCollection collection.map(function(img) { var ndvi img.normalizedDifference([B8, B4]).rename(NDVI); return ndvi.copyProperties(img, [system:time_start]); }); // 5. 年度NDVI均值合成 var ndviAnnualMean ndviCollection.mean();先说CLOUDY_PIXEL_PERCENTAGE这个属性。GEE给每景Sentinel-2影像都配了一个元数据表示整景影像上云量占比。我在这里设了20%以内的筛选是一个综合考虑。设太低像5%很多月份可能一景影像都没有尤其是热带和云多的地区设太高像50%云掩膜后有效像元太少NDVI年均值会有大片空洞。20%对大多数地区来说是一个不错的起点如果你研究区特别干燥比如戈壁、荒漠可以直接放宽到50%因为那些地方云少云掩膜基本不损失像元。然后是.map(maskS2Clouds)之后集合里每张影像只剩B2、B3、B4、B8和掩膜信息。normalizedDifference([B8, B4])就是算(B8-B4)/(B8B4)输出单波段影像叫NDVI。这里有个小细节copyProperties(img, [system:time_start])这一步很多人会漏掉。不加它每张NDVI影像就丢了时间属性后面做时间序列图表会出问题。虽然我们这里只是做简单均值但养成保留时间属性的习惯以后代码复用起来省无数事。那.mean()做的事情就是逐像元地把一年内所有有效NDVI值加起来除以有效观测次数。注意它对不同日期是等权平均没有做季节加权。这就意味着如果某个区域夏季影像多、冬季影像少年均值会偏向夏季状态。这是简单均值法的天然缺陷但作为通用的“NDVI年均值产品”来说这种偏差在大多数场景下可以接受。如果你需要更严谨的产品可以考虑用月合成后再年均或者做时间序列的谐波拟合那是更高阶的话题这篇不展开。3.4 核心代码四可视化与客户端下载结果算出来了不下载到本地等于白干。GEE现在有两种主流下载方式我先说可视化再分别讲两种下载方案。// 6. 可视化 var visParam { min: -0.1, max: 0.9, palette: [#d73027, #f46d43, #fdae61, #fee08b, #d9ef8b, #a6d96a, #66bd63, #1a9850, #006837] }; Map.addLayer(ndviAnnualMean.clip(roi), visParam, NDVI Annual Mean);这个色带是从红到绿的渐变负值红色裸土黄绿色高植被深绿色。做NDVI图我很推荐这个配色比默认的黑白灰直观得多。然后是下载。先给第一种方式——如果你用的是GEE Python环境或者geemap可以用geemap.ee_to_filename之类的函数直接把影像拉到本地。但在浏览器代码编辑器里最稳的还是传统Export// 7. 导出到Google Drive Export.image.toDrive({ image: ndviAnnualMean.clip(roi).rename(NDVI_ year), description: NDVI_Annual_Mean_ year, folder: GEE_Exports, region: roi, scale: 10, crs: EPSG:4326, maxPixels: 1e13 });scale: 10对应Sentinel-2原始分辨率。如果你想跑得快一点可以用20米精度损失不太大。crs我强烈建议显式指定EPSG:4326不然GEE会自动选择影像的原始UTM投影某些时候你拿回来的影像在ArcGIS里会跟其他数据对不齐。第二种下载方式是2022年之后GEE新推的客户端直下载方式在代码编辑器里配合getDownloadURL// 7. 备选直接下载GeoTIFF到本地 var downloadUrl ndviAnnualMean.clip(roi) .select(NDVI) .getDownloadURL({ scale: 10, crs: EPSG:4326, region: roi, format: GeoTIFF }); print(下载链接复制到浏览器打开:, downloadUrl);这两种方案我用下来的体会是小范围比如一个县级市用getDownloadURL更快省去登录Drive再下载的环节大范围或者一次跑多块区域老老实实用Export.image.toDrive批量挂任务不然浏览器标签页会直接卡死。3.5 完整源码汇总上面几段代码合在一起就是一套可以完整运行的源代码——你如果跳过了前面的拆解直接想复制下面这份就是最终整合版// GEE学习笔记29Sentinel-2高精度NDVI年均值下载 // 适用平台Google Earth Engine Code Editor // 1. 研究区设置 var roi ee.FeatureCollection(users/your_username/your_roi); Map.centerObject(roi, 8); Map.addLayer(roi, {color: red}, 研究区); // 2. 时间范围 var year 2023; var startDate ee.Date.fromYMD(year, 1, 1); var endDate ee.Date.fromYMD(year 1, 1, 1); // 3. 云掩膜函数 function maskS2Clouds(image) { var qa image.select(QA60); var cloudBitMask 1 10; var cirrusBitMask 1 11; var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); return image.updateMask(mask).select([B2, B3, B4, B8]); } // 4. 影像集合加载与NDVI逐景计算 var collection ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) .map(maskS2Clouds); var ndviCollection collection.map(function(img) { var ndvi img.normalizedDifference([B8, B4]).rename(NDVI); return ndvi.copyProperties(img, [system:time_start]); }); // 5. 年度均值 var ndviAnnualMean ndviCollection.mean(); // 6. 可视化 var visParam { min: -0.1, max: 0.9, palette: [#d73027, #f46d43, #fdae61, #fee08b, #d9ef8b, #a6d96a, #66bd63, #1a9850, #006837] }; Map.addLayer(ndviAnnualMean.clip(roi), visParam, NDVI Annual Mean year); // 7. 导出默认走Google Drive Export.image.toDrive({ image: ndviAnnualMean.clip(roi).rename(NDVI_ year), description: NDVI_Annual_Mean_ year, folder: GEE_Exports, region: roi, scale: 10, crs: EPSG:4326, maxPixels: 1e13 });4. 实操过程中的常见问题与排查实录这部分必须写因为我在调这套脚本的过程中踩过的坑确实不少。下面的内容比一般的FAQ更具体都是逐个排查出来的实战心得。4.1 云掩膜后大面积区域没有值最典型的症状是NDVI年均值出图后研究区大片是空白或者只有零星几个像元有值。原因多半是云掩膜太狠了。最常见的一个错误是有人从网上抄了一段QA60掩膜函数里面的cloudBitMask和cirrusBitMask定义没问题但忘了在map(maskS2Clouds)之前过滤云量导致那些60%云量的影像也被拉进来算。这些影像经过云掩膜后整个研究区内几乎没有有效像元但因为GEE的mean()不会把全空影像剔除所以最后结果空掉。解决办法有两层。第一层是在集合筛选时加ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)这已经是源码里的一部分。第二层是加一个“有效观测次数”波段把每年逐像元的有效观测数也统计出来后处理时你可以看到每个像元到底用了多少景影像// 额外统计有效观测次数 var countCollection collection.map(function(img) { return img.select(B8).updateMask(img.select(B8)).rename(valid); }); var validCount countCollection.count(); Map.addLayer(validCount.clip(roi), {min: 0, max: 50}, 有效影像数量);我实际跑下来的经验是如果某个像元一年里有效观测少于5次它的NDVI年均值噪点会非常严重建议后续处理中直接设无效值。你可以用ndviAnnualMean.updateMask(validCount.gte(5))来产出更可靠的产品。4.2 年均值出现异常低值或负值还有一个特别迷惑的问题明明研究区全是森林年均NDVI却跑出0.1以下。这个现象在山区特别常见。根因有两个一是残留云没消干净二是地形阴影。QA60的云掩膜对厚云有效但对山体阴影和半透明薄云经常漏判。这些像元的表观反射率特征是红光偏高、近红外偏低NDVI算出来就是负值直接拉低年均值。我建议的排查流程是下载结果后在ENVI或者QGIS里叠加“有效影像数量”图层看低值区域是否集中在地形起伏大的地方。如果是那就得在脚本里加地形阴影掩膜。比较简单的做法是利用SRTM高程数据计算坡度把坡度大于30度或者太阳入射角极低的区域剔除掉。这个做法不完美但对大多数应用场景够用// 坡度掩膜针对山地地形 var dem ee.Image(USGS/SRTMGL1_003); var slope ee.Terrain.slope(dem); var terrainMask slope.lt(30); var ndviAnnualMeanFixed ndviAnnualMean.updateMask(terrainMask);你还可以把云掩膜条件里的卷云位去掉只保留不透明云位因为卷云位会把一部分地形高亮区误伤。一句话总结先做有效计数再看影像成因别一上来就怀疑算法。4.3 下载到本地的GeoTIFF坐标系不一致这是一个让人非常抓狂的问题。你在GEE里显示正常下载到本地后用Python的rasterio或者ArcGIS打开发现它跟你的其他矢量边界对不上位置偏移了几百米到几公里。这个问题的根源是CRS不统一。GEE底层大部分操作在等经纬度投影EPSG:4326或者影像原始UTM投影下进行。你如果不显式设置crs参数导出时GEE会按照影像自身的坐标系输出。而你的ROI可能是CGCS2000或者WGS84 UTM投影两者混合在一起就会产生偏移。我在源码里已经显式写了crs: EPSG:4326这样输出的GeoTIFF是经纬度坐标在大多数软件里都能正确识别。如果你需要其他坐标系直接在导出参数里修改crs即可比如EPSG:32650是UTM 50N。这里有个容易栽跟头的点region必须和crs匹配如果你crs设置了UTMregion最好也是同一投影下的边界否则导出的图会像“揉皱的纸”一样扭曲。4.4 导出任务跑太久或内存上限报错maxPixels这个参数如果你不设置GEE默认导出上限是1千万像素。10米分辨率下一景100km×100km的图像素数是1亿直接超限报错。源码里设置了maxPixels: 1e13足够大多数区域使用。还有更麻烦的“Computation timed out”问题。这个问题通常出在影像集合过大或云掩膜后有效像元太少、合成效率低下。我自己的经验是如果研究区范围跨越省域级别就不要试图一次性算出全国DEM级别的NDVI均值分批分区做是更明智的选择。另外如果你只是在代码编辑器里预览结果不建议让GEE实时计算全图应该把Map.addLayer的显示分辨率调低或者直接用roi去看中心区域。跟GEE打交道久了你会发现“为什么我的地图加载那么慢”的答案往往是不是GEE慢是你在浏览器里加载了一个几十GB的影像。5. 几段能直接改来用的进阶变体源码跑通只是第一步做研究的人肯定要按自己的需求改造它。我把自己常用的几个变体也一并写出来这些内容是我自己项目里的沉淀比单纯的官网示例有用得多。5.1 完整年份脚本改造成季度或月份均值有时候你不想看年均值而想看春季、夏季、秋季各自的NDVI。改法非常简单把时间范围切一切就可以了“代码三”里的filterDate改一下比如春季是3月到5月。但这里也有一个容易被忽略的问题月份的划分和影像数量的季节差异会直接影响结果的可比性。我更推荐的做法是把全年12个月的月均NDVI都算出来再根据需要合成季节值。这样你可以顺便输出一个12波段的NDVI堆栈以后做物候分析也能用上。代码可以这么扩展var monthList ee.List.sequence(1, 12); var monthlyNDVI ee.ImageCollection(monthList.map(function(m) { var month ee.Number(m); var start ee.Date.fromYMD(year, month, 1); var end start.advance(1, month); var monthly ndviCollection.filterDate(start, end).mean(); return monthly.set(month, month); }));这个输出后每个像元有12个值对应一年12个月的NDVI均值。配合ee.ImageCollection.toBands()转成多波段影像直接导出分析都非常方便。5.2 使用分位数合成替代均值前面一直强调均值但均值对极端值敏感。有些年份某些地区的云掩膜不彻底残留一两个负的异常NDVI像元均值就会被拉低。稳健的做法是用median()代替mean()或者使用分位数合成比如90%分位数代表生长季峰值。我自己的经验法则做资源调查、生态质量评价这类偏“状态评估”的应用时用median()往往比mean()更稳因为中位数不受少数异常值影响。做碳循环建模或生物量估算时mean()和无云数据配合使用更合适因为它能保留更多季节细节。// 中位数版本一行替换 var ndviAnnualMedian ndviCollection.median();别小看这一行在很多用GEE处理Landsat/Sentinel时间序列的老手手里median才是默认选项。只不过这篇笔记里标题强调的是“年均值”所以主代码用了mean。5.3 结合其他指数做多维分析NDVI年均值通常不是终点而是一个输入。比如你可以把它和MNDWI水体指数、BSI裸土指数结合做土地利用分类的特征层。或者把它和DEM叠加做植被垂直分布分析。我给一个非常实用的组合NDVI年均值 地表温度LST年均值可以在GEE里用MODIS的MOD11A2产品快速生成然后做散点图验证植被与温度的耦合关系。这已经是生态遥感里比较经典的组合分析了唯一要注意的是两个产品的分辨率差异NDVI是10mLST是1km需要先对LST重采样到NDVI的网格上或者反过来把NDVI聚合到1km。6. 几条藏在细节里的经验聊完收工这篇文章的正文到这里已经足够长了最后我再说几句只有实际操作过的人才能体会的东西。第一GEE代码编辑器里不要把print()用在一整个大影像上会直接让浏览器卡死。你要看某个像元的值用sample()或者reduceRegion()指定坐标点输出数字即可。第二下载结果之后如果你要用ArcGIS或者QGIS做后续分析建议把NDVI值域从float类型重新标定到0到10000的整数乘以10000取整既节省硬盘空间又避免浮点数的显示精度问题。GEE导出时可以加一句.multiply(10000).toInt()但要注意后续使用时的说明不然别人拿到数据看不懂值域会困惑。第三也是最想强调的GEE只是工具不是算法真理。年均NDVI这个产品换一个数据源Landsat8、MODIS、换一种合成方式最大值、中位数、均值、换一个时间范围生长季、全年结果都不一样。做研究的人必须在文章里写清楚自己的技术路线不然别人复现不了那你这个“一键下载”对别人来说就是黑箱。这套源码本身不复杂我希望通过这次拆解你能把思路捋清楚。等你看懂了为什么先算NDVI再取平均、为什么要对云掩膜保持敬畏、为什么导出要指定坐标系以后不管是用Landsat还是用其他传感器都能自己写出更合适的版本。如果照着跑遇到什么问题再去翻一翻GEE官方文档里normalizedDifference和updateMask的说明基本就能解决。祝你的NDVI图一次出图成功。 GEE实战：Sentinel-2高精度NDVI年均值一键下载脚本详解
1. 先把话说清楚这套脚本到底解决什么问题做遥感的人应该都有这种感觉Sentinel-2数据量大、时间序列密、波段又多真要逐景下载再进ENVI或者ArcGIS里做NDVI、做年均值光是文件整理就能把人搞到怀疑人生。尤其是NDVI这种基础指数大家天天用可一旦涉及到“长时间序列”“大范围”“批量处理”这几个词传统桌面软件的工作流马上变得笨重不堪。GEEGoogle Earth Engine的价值就在这你不用下载任何一景影像所有计算都在云端完成最后只把结果拿回来。这句“一键下载Sentinel-2高精度NDVI年均值数据”本质上是把过去需要好几个小时、好几个软件配合才能完成的流程压缩成一个脚本、一次点击。这篇笔记里我会把自己反复调试后的一套完整源码放出来并且把所有关键判断和坑都写清楚。适合的对象很明确已经在GEE里跑过简单例程、知道Map.addLayer和print基础操作但第一次想正经做“年度NDVI产品”并下载到本地的朋友。就算你对GEE还比较生疏只要照着源码逐行抄也能跑通。先说结论这个脚本的核心逻辑是筛选某一年内所有Sentinel-2影像做云掩膜处理逐景计算NDVI然后对全年NDVI影像求平均值最后导出成GeoTIFF到本地。注意顺序——是先算NDVI再聚合年均值不是先合成影像再算NDVI这个细节很重要后面我会专门解释。2. 方案设计为什么是“先算NDVI再取平均”2.1 NDVI原理和Sentinel-2波段对照NDVI的公式没什么好回避的就是那个经典的归一化差值NDVI (NIR - Red) / (NIR Red)在Sentinel-2 L2A产品里红光对应的是B4中心波长665nm10米分辨率近红外对应的是B8中心波长833nm10米分辨率。这两个波段都是10米空间分辨率所以直接相减相除没有分辨率不匹配的问题。NDVI的值域范围是-1到1。裸土、水体、建筑这些地物通常在0附近甚至为负植被一般在0.3到0.9之间。我们做“年均值”本质上就是把全年不同季节的植被状态压到一张图上反映的是“这一年里这片地方整体的绿度水平”。农田、林地、草地它们在不同季节的NDVI曲线差异非常大年均值能很好地把这种差异体现出来这也是我为什么觉得年度NDVI产品比单期影像更“高精度”——它不是一个瞬间的快照而是一整年的平均状态抗偶然性更强。2.2 两种聚合策略的取舍这里我必须掰开揉碎讲清楚因为这是整套脚本最核心的设计决策。第一种策略先把某年的所有影像按时间合成一张“全年代表影像”比如用median、mean或qualityMosaic然后拿这张合成影像去算NDVI。这种做法的优点是计算量小一次合成一次NDVI速度快。但缺点也很明显合成影像的像元值是光谱反射率在时间维上的统计结果拿统计后的光谱去算NDVI本质上是在算“平均光谱的NDVI”而不是“NDVI的平均”。第二种策略每景影像先算NDVI得到一整年的NDVI影像集合再对这个集合取平均。这种情况下每个像元参与平均的都是“瞬时NDVI值”最终得到的是真正意义上的“全年NDVI均值”。数字上两者会有明显差异而且第二种更贴近植被生态学的定义。我做过对比同一个地方第一种方法算出来常绿林区的NDVI可能差0.05到0.1别小看这零点几做长时序分析的时候这点偏差足够颠倒结论了。所以我的建议是在你想要表达“这一年植被平均生长状况”时一律用先算NDVI再求均值的方式。这也是这个脚本名字里“高精度”三个字的底气来源。2.3 数据源选择L2A还是SRGEE里现在有两套Sentinel-2多光谱数据在使用一个叫COPERNICUS/S2_SR另一个叫COPERNICUS/S2_SR_HARMONIZED。很多新手踩过的坑就是用了前者然后发现2022年1月之后的数据出现明显的辐射不一致问题。原因不复杂2022年1月底欧空局调整了Sentinel-2的辐射定标参数导致新旧影像的表面反射率存在一个系统性的偏移。GEE官方为了解决这个问题推出了HARMONIZED版本把所有影像都统一到了新的定标标准下。你现在新建脚本在代码编辑器里搜索Sentinel-2官方推荐也是用HARMONIZED版本。所以我这个脚本直接用COPERNICUS/S2_SR_HARMONIZED你不需要自己处理跨年份的辐射一致性它在数据集层面已经帮你调和过了。顺便说一句有人还在用TOACOPERNICUS/S2产品做NDVI。不是不能做但TOA没做过大气校正受气溶胶和大气水汽影响很大山区多云地区的NDVI波动会很夸张。既然L2A表面反射率产品已经免费开放就别再用TOA委屈自己了。3. 完整源码与逐行拆解3.1 核心代码一研究区与时间范围设置先把完整源码贴出来然后我逐段拆解。这段源码我已经在GEE代码编辑器里实测跑通你复制进去就能用。// 1. 研究区设置 // 这里随便给了一个矢量范围实际使用请换成自己的研究区边界 var roi ee.FeatureCollection(users/your_username/your_roi); Map.centerObject(roi, 8); Map.addLayer(roi, {color: red}, 研究区); // 2. 时间范围设置 var year 2023; var startDate ee.Date.fromYMD(year, 1, 1); var endDate ee.Date.fromYMD(year 1, 1, 1);这一段没什么难度但有两个细节我想提醒你。第一ee.Date.fromYMD(year, 1, 1)生成的是这一年的1月1日0点而endDate取的是year 1年的1月1日也就是说结束时间点其实是2024年1月1日0点。GEE里filterDate是“左闭右开”的所以写成这样恰好覆盖了整年365天闰年366天。如果你图省事写成ee.Date.fromYMD(year, 12, 31)就会丢掉12月31日那天的所有影像别问我怎么知道的。第二ROI推荐用gee FeatureCollection而不是ee.Geometry画个矩形。原因是Export.image.toDrive这类导出任务对Region参数有严格限制有时候你用ee.Geometry.Rectangle画的范围跟影像的坐标系不匹配导出的结果会出现网格错位。用矢量边界哪怕是简单的Shapefile上传会稳定得多。3.2 核心代码二云掩膜函数Sentinel-2 L2A产品自带一个QA波段叫QA60它的第10位和第11位分别记录了云和卷云的信息。用位运算把它们提取出来就能生成云掩膜。我封装了一个函数这段是我调试最多的地方// 3. 云掩膜函数 function maskS2Clouds(image) { var qa image.select(QA60); // 第10位不透明云第11位卷云 var cloudBitMask 1 10; var cirrusBitMask 1 11; // 这两个bit位为1说明是云要mask掉 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); // 更新影像掩膜后同时把QA60波段丢掉保持数据整洁 return image.updateMask(mask).select([B2, B3, B4, B8]); }bitwiseAnd是按位与运算。比如某像元的QA60值是2048二进制就是第11位为1按位与2048后结果不为0说明这个像元是云所以mask里它就会被设为false。同理第10位是1024。两个位都为0才判定为清晰像元。我用下来觉得这个QA60掩膜对厚云效果很好但对很薄的卷云偶尔会漏判这是L2A算法的固有情况不算bug。另外如果你研究区是高纬度地区冬天太阳高度角低QA60会把一些山体阴影也误判为云这时候可以放宽掩膜条件比如只处理第10位不透明云舍弃卷云掩膜。脚本里我保留了卷云位因为我们的目标是做全年NDVI均值对云的要求尽量严格一些。3.3 核心代码三NDVI计算与年度合成这是整个脚本的关键环节也是设计思路里“先算NDVI再取均值”的落地实现// 4. 加载影像集合逐景算NDVI var collection ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) .map(maskS2Clouds); // 逐景计算NDVI先算NDVI var ndviCollection collection.map(function(img) { var ndvi img.normalizedDifference([B8, B4]).rename(NDVI); return ndvi.copyProperties(img, [system:time_start]); }); // 5. 年度NDVI均值合成 var ndviAnnualMean ndviCollection.mean();先说CLOUDY_PIXEL_PERCENTAGE这个属性。GEE给每景Sentinel-2影像都配了一个元数据表示整景影像上云量占比。我在这里设了20%以内的筛选是一个综合考虑。设太低像5%很多月份可能一景影像都没有尤其是热带和云多的地区设太高像50%云掩膜后有效像元太少NDVI年均值会有大片空洞。20%对大多数地区来说是一个不错的起点如果你研究区特别干燥比如戈壁、荒漠可以直接放宽到50%因为那些地方云少云掩膜基本不损失像元。然后是.map(maskS2Clouds)之后集合里每张影像只剩B2、B3、B4、B8和掩膜信息。normalizedDifference([B8, B4])就是算(B8-B4)/(B8B4)输出单波段影像叫NDVI。这里有个小细节copyProperties(img, [system:time_start])这一步很多人会漏掉。不加它每张NDVI影像就丢了时间属性后面做时间序列图表会出问题。虽然我们这里只是做简单均值但养成保留时间属性的习惯以后代码复用起来省无数事。那.mean()做的事情就是逐像元地把一年内所有有效NDVI值加起来除以有效观测次数。注意它对不同日期是等权平均没有做季节加权。这就意味着如果某个区域夏季影像多、冬季影像少年均值会偏向夏季状态。这是简单均值法的天然缺陷但作为通用的“NDVI年均值产品”来说这种偏差在大多数场景下可以接受。如果你需要更严谨的产品可以考虑用月合成后再年均或者做时间序列的谐波拟合那是更高阶的话题这篇不展开。3.4 核心代码四可视化与客户端下载结果算出来了不下载到本地等于白干。GEE现在有两种主流下载方式我先说可视化再分别讲两种下载方案。// 6. 可视化 var visParam { min: -0.1, max: 0.9, palette: [#d73027, #f46d43, #fdae61, #fee08b, #d9ef8b, #a6d96a, #66bd63, #1a9850, #006837] }; Map.addLayer(ndviAnnualMean.clip(roi), visParam, NDVI Annual Mean);这个色带是从红到绿的渐变负值红色裸土黄绿色高植被深绿色。做NDVI图我很推荐这个配色比默认的黑白灰直观得多。然后是下载。先给第一种方式——如果你用的是GEE Python环境或者geemap可以用geemap.ee_to_filename之类的函数直接把影像拉到本地。但在浏览器代码编辑器里最稳的还是传统Export// 7. 导出到Google Drive Export.image.toDrive({ image: ndviAnnualMean.clip(roi).rename(NDVI_ year), description: NDVI_Annual_Mean_ year, folder: GEE_Exports, region: roi, scale: 10, crs: EPSG:4326, maxPixels: 1e13 });scale: 10对应Sentinel-2原始分辨率。如果你想跑得快一点可以用20米精度损失不太大。crs我强烈建议显式指定EPSG:4326不然GEE会自动选择影像的原始UTM投影某些时候你拿回来的影像在ArcGIS里会跟其他数据对不齐。第二种下载方式是2022年之后GEE新推的客户端直下载方式在代码编辑器里配合getDownloadURL// 7. 备选直接下载GeoTIFF到本地 var downloadUrl ndviAnnualMean.clip(roi) .select(NDVI) .getDownloadURL({ scale: 10, crs: EPSG:4326, region: roi, format: GeoTIFF }); print(下载链接复制到浏览器打开:, downloadUrl);这两种方案我用下来的体会是小范围比如一个县级市用getDownloadURL更快省去登录Drive再下载的环节大范围或者一次跑多块区域老老实实用Export.image.toDrive批量挂任务不然浏览器标签页会直接卡死。3.5 完整源码汇总上面几段代码合在一起就是一套可以完整运行的源代码——你如果跳过了前面的拆解直接想复制下面这份就是最终整合版// GEE学习笔记29Sentinel-2高精度NDVI年均值下载 // 适用平台Google Earth Engine Code Editor // 1. 研究区设置 var roi ee.FeatureCollection(users/your_username/your_roi); Map.centerObject(roi, 8); Map.addLayer(roi, {color: red}, 研究区); // 2. 时间范围 var year 2023; var startDate ee.Date.fromYMD(year, 1, 1); var endDate ee.Date.fromYMD(year 1, 1, 1); // 3. 云掩膜函数 function maskS2Clouds(image) { var qa image.select(QA60); var cloudBitMask 1 10; var cirrusBitMask 1 11; var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); return image.updateMask(mask).select([B2, B3, B4, B8]); } // 4. 影像集合加载与NDVI逐景计算 var collection ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) .map(maskS2Clouds); var ndviCollection collection.map(function(img) { var ndvi img.normalizedDifference([B8, B4]).rename(NDVI); return ndvi.copyProperties(img, [system:time_start]); }); // 5. 年度均值 var ndviAnnualMean ndviCollection.mean(); // 6. 可视化 var visParam { min: -0.1, max: 0.9, palette: [#d73027, #f46d43, #fdae61, #fee08b, #d9ef8b, #a6d96a, #66bd63, #1a9850, #006837] }; Map.addLayer(ndviAnnualMean.clip(roi), visParam, NDVI Annual Mean year); // 7. 导出默认走Google Drive Export.image.toDrive({ image: ndviAnnualMean.clip(roi).rename(NDVI_ year), description: NDVI_Annual_Mean_ year, folder: GEE_Exports, region: roi, scale: 10, crs: EPSG:4326, maxPixels: 1e13 });4. 实操过程中的常见问题与排查实录这部分必须写因为我在调这套脚本的过程中踩过的坑确实不少。下面的内容比一般的FAQ更具体都是逐个排查出来的实战心得。4.1 云掩膜后大面积区域没有值最典型的症状是NDVI年均值出图后研究区大片是空白或者只有零星几个像元有值。原因多半是云掩膜太狠了。最常见的一个错误是有人从网上抄了一段QA60掩膜函数里面的cloudBitMask和cirrusBitMask定义没问题但忘了在map(maskS2Clouds)之前过滤云量导致那些60%云量的影像也被拉进来算。这些影像经过云掩膜后整个研究区内几乎没有有效像元但因为GEE的mean()不会把全空影像剔除所以最后结果空掉。解决办法有两层。第一层是在集合筛选时加ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)这已经是源码里的一部分。第二层是加一个“有效观测次数”波段把每年逐像元的有效观测数也统计出来后处理时你可以看到每个像元到底用了多少景影像// 额外统计有效观测次数 var countCollection collection.map(function(img) { return img.select(B8).updateMask(img.select(B8)).rename(valid); }); var validCount countCollection.count(); Map.addLayer(validCount.clip(roi), {min: 0, max: 50}, 有效影像数量);我实际跑下来的经验是如果某个像元一年里有效观测少于5次它的NDVI年均值噪点会非常严重建议后续处理中直接设无效值。你可以用ndviAnnualMean.updateMask(validCount.gte(5))来产出更可靠的产品。4.2 年均值出现异常低值或负值还有一个特别迷惑的问题明明研究区全是森林年均NDVI却跑出0.1以下。这个现象在山区特别常见。根因有两个一是残留云没消干净二是地形阴影。QA60的云掩膜对厚云有效但对山体阴影和半透明薄云经常漏判。这些像元的表观反射率特征是红光偏高、近红外偏低NDVI算出来就是负值直接拉低年均值。我建议的排查流程是下载结果后在ENVI或者QGIS里叠加“有效影像数量”图层看低值区域是否集中在地形起伏大的地方。如果是那就得在脚本里加地形阴影掩膜。比较简单的做法是利用SRTM高程数据计算坡度把坡度大于30度或者太阳入射角极低的区域剔除掉。这个做法不完美但对大多数应用场景够用// 坡度掩膜针对山地地形 var dem ee.Image(USGS/SRTMGL1_003); var slope ee.Terrain.slope(dem); var terrainMask slope.lt(30); var ndviAnnualMeanFixed ndviAnnualMean.updateMask(terrainMask);你还可以把云掩膜条件里的卷云位去掉只保留不透明云位因为卷云位会把一部分地形高亮区误伤。一句话总结先做有效计数再看影像成因别一上来就怀疑算法。4.3 下载到本地的GeoTIFF坐标系不一致这是一个让人非常抓狂的问题。你在GEE里显示正常下载到本地后用Python的rasterio或者ArcGIS打开发现它跟你的其他矢量边界对不上位置偏移了几百米到几公里。这个问题的根源是CRS不统一。GEE底层大部分操作在等经纬度投影EPSG:4326或者影像原始UTM投影下进行。你如果不显式设置crs参数导出时GEE会按照影像自身的坐标系输出。而你的ROI可能是CGCS2000或者WGS84 UTM投影两者混合在一起就会产生偏移。我在源码里已经显式写了crs: EPSG:4326这样输出的GeoTIFF是经纬度坐标在大多数软件里都能正确识别。如果你需要其他坐标系直接在导出参数里修改crs即可比如EPSG:32650是UTM 50N。这里有个容易栽跟头的点region必须和crs匹配如果你crs设置了UTMregion最好也是同一投影下的边界否则导出的图会像“揉皱的纸”一样扭曲。4.4 导出任务跑太久或内存上限报错maxPixels这个参数如果你不设置GEE默认导出上限是1千万像素。10米分辨率下一景100km×100km的图像素数是1亿直接超限报错。源码里设置了maxPixels: 1e13足够大多数区域使用。还有更麻烦的“Computation timed out”问题。这个问题通常出在影像集合过大或云掩膜后有效像元太少、合成效率低下。我自己的经验是如果研究区范围跨越省域级别就不要试图一次性算出全国DEM级别的NDVI均值分批分区做是更明智的选择。另外如果你只是在代码编辑器里预览结果不建议让GEE实时计算全图应该把Map.addLayer的显示分辨率调低或者直接用roi去看中心区域。跟GEE打交道久了你会发现“为什么我的地图加载那么慢”的答案往往是不是GEE慢是你在浏览器里加载了一个几十GB的影像。5. 几段能直接改来用的进阶变体源码跑通只是第一步做研究的人肯定要按自己的需求改造它。我把自己常用的几个变体也一并写出来这些内容是我自己项目里的沉淀比单纯的官网示例有用得多。5.1 完整年份脚本改造成季度或月份均值有时候你不想看年均值而想看春季、夏季、秋季各自的NDVI。改法非常简单把时间范围切一切就可以了“代码三”里的filterDate改一下比如春季是3月到5月。但这里也有一个容易被忽略的问题月份的划分和影像数量的季节差异会直接影响结果的可比性。我更推荐的做法是把全年12个月的月均NDVI都算出来再根据需要合成季节值。这样你可以顺便输出一个12波段的NDVI堆栈以后做物候分析也能用上。代码可以这么扩展var monthList ee.List.sequence(1, 12); var monthlyNDVI ee.ImageCollection(monthList.map(function(m) { var month ee.Number(m); var start ee.Date.fromYMD(year, month, 1); var end start.advance(1, month); var monthly ndviCollection.filterDate(start, end).mean(); return monthly.set(month, month); }));这个输出后每个像元有12个值对应一年12个月的NDVI均值。配合ee.ImageCollection.toBands()转成多波段影像直接导出分析都非常方便。5.2 使用分位数合成替代均值前面一直强调均值但均值对极端值敏感。有些年份某些地区的云掩膜不彻底残留一两个负的异常NDVI像元均值就会被拉低。稳健的做法是用median()代替mean()或者使用分位数合成比如90%分位数代表生长季峰值。我自己的经验法则做资源调查、生态质量评价这类偏“状态评估”的应用时用median()往往比mean()更稳因为中位数不受少数异常值影响。做碳循环建模或生物量估算时mean()和无云数据配合使用更合适因为它能保留更多季节细节。// 中位数版本一行替换 var ndviAnnualMedian ndviCollection.median();别小看这一行在很多用GEE处理Landsat/Sentinel时间序列的老手手里median才是默认选项。只不过这篇笔记里标题强调的是“年均值”所以主代码用了mean。5.3 结合其他指数做多维分析NDVI年均值通常不是终点而是一个输入。比如你可以把它和MNDWI水体指数、BSI裸土指数结合做土地利用分类的特征层。或者把它和DEM叠加做植被垂直分布分析。我给一个非常实用的组合NDVI年均值 地表温度LST年均值可以在GEE里用MODIS的MOD11A2产品快速生成然后做散点图验证植被与温度的耦合关系。这已经是生态遥感里比较经典的组合分析了唯一要注意的是两个产品的分辨率差异NDVI是10mLST是1km需要先对LST重采样到NDVI的网格上或者反过来把NDVI聚合到1km。6. 几条藏在细节里的经验聊完收工这篇文章的正文到这里已经足够长了最后我再说几句只有实际操作过的人才能体会的东西。第一GEE代码编辑器里不要把print()用在一整个大影像上会直接让浏览器卡死。你要看某个像元的值用sample()或者reduceRegion()指定坐标点输出数字即可。第二下载结果之后如果你要用ArcGIS或者QGIS做后续分析建议把NDVI值域从float类型重新标定到0到10000的整数乘以10000取整既节省硬盘空间又避免浮点数的显示精度问题。GEE导出时可以加一句.multiply(10000).toInt()但要注意后续使用时的说明不然别人拿到数据看不懂值域会困惑。第三也是最想强调的GEE只是工具不是算法真理。年均NDVI这个产品换一个数据源Landsat8、MODIS、换一种合成方式最大值、中位数、均值、换一个时间范围生长季、全年结果都不一样。做研究的人必须在文章里写清楚自己的技术路线不然别人复现不了那你这个“一键下载”对别人来说就是黑箱。这套源码本身不复杂我希望通过这次拆解你能把思路捋清楚。等你看懂了为什么先算NDVI再取平均、为什么要对云掩膜保持敬畏、为什么导出要指定坐标系以后不管是用Landsat还是用其他传感器都能自己写出更合适的版本。如果照着跑遇到什么问题再去翻一翻GEE官方文档里normalizedDifference和updateMask的说明基本就能解决。祝你的NDVI图一次出图成功。