用 Earth Engine 绘制维也纳夏季无云影像哨兵2数据筛选与 NDVI 计算实战本文带你一步步解析一段 Earth Engine 代码从指定区域、筛选影像、去云、计算植被指数到可视化最终获取维也纳上空 2024 年夏季最晴朗的哨兵2合成影像。一、为什么要写这段代码遥感数据处理往往面临两个痛点数据量大哨兵2每天全球产生数 TB 数据和云层遮挡光学影像受天气影响严重。Google Earth EngineGEE提供了云端并行计算能力让我们能在浏览器中快速完成大规模影像的筛选、预处理和分析。今天我们要实现的目的是锁定维也纳经纬度 16.3738°E, 48.2082°N周边区域获取 2024 年 6 月 1 日至 9 月 1 日期间所有哨兵2 Level-2A 地表反射率影像剔除云量超过 20% 的影像为每景影像计算NDVI归一化植被指数并附上标签按云量从低到高排序最终合成一张“最清晰”的 mosaic 影像在交互地图上展示真彩色合成并打印出一些关键元数据。二、代码分步拆解1. 定义研究区域varviennaee.Geometry.Point([16.3738,48.2082]);我们创建一个点几何对象代表维也纳市中心。这个点将作为空间筛选的基准。2. 加载哨兵2影像集合varsentinel2ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED).filterBounds(vienna).filterDate(2024-06-01,2024-09-01).filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE,20))使用官方提供的COPERNICUS/S2_SR_HARMONIZED数据集该数据集已进行辐射定标和大气校正且将 Sentinel-2A 与 2B 的波段进行了统一“HARMONIZED”。空间筛选只保留覆盖维也纳点的影像。时间筛选限定 2024 年夏季北半球植被生长旺季。云量筛选CLOUDY_PIXEL_PERCENTAGE是每景影像自带的元数据字段这里要求小于 20%避免云层干扰后续分析。3. 计算 NDVI 并添加为波段.map(function(image){varndviimage.normalizedDifference([B8,B4]).rename(NDVI);returnimage.addBands(ndvi).set(summary,summer clear-sky candidate);})对集合中的每一景影像使用normalizedDifference方法计算 (B8 - B4) / (B8 B4)即 NDVI。其中 B8 为近红外波段NIRB4 为红波段。将新生成的 NDVI 波段添加到原影像中并设置一个自定义属性summary方便后续识别这批影像的来源。4. 排序并合成.sort(CLOUDY_PIXEL_PERCENTAGE);varfirstClearImagesentinel2.mosaic();按照云量百分比升序排列云量越少越靠前。调用mosaic()将集合中的所有影像合成一张。由于云量已排序mosaic 会按照集合顺序叠加最后显示在最上层的是云量最少的像素但实际 mosaic 的规则是后加入的覆盖先加入的所以排序后低云影像会覆盖高云影像从而得到整体最清晰的合成图。这里严格来说mosaic 并不一定是“第一张清晰影像”而是所有影像的镶嵌结果但结合排序效果相当于优先使用低云像素。5. 地图可视化Map.addLayer(firstClearImage.select([B4,B3,B2]),{bands:[B4,B3,B2],min:0,max:3000},Sentinel-2 true color,true,0.7);选择红B4、绿B3、蓝B2三个波段合成真彩色影像。拉伸范围设为 03000对应地表反射率缩放因子原始值为整数一般乘以 0.0001 得到反射率但这里直接显示。图层透明度设为 0.7方便与底图叠加观察。6. 打印调试信息print(Loaded collection,COPERNICUS/S2_SR_HARMONIZED);print(Filtered image count,sentinel2.size());print(Scene ids,sentinel2.limit(4).aggregate_array(system:index));print(Cloud percentages,sentinel2.limit(4).aggregate_array(CLOUDY_PIXEL_PERCENTAGE));print(First clear image,firstClearImage);print(Selected bands from first 3,sentinel2.limit(3).select([B4,B8,NDVI]));输出集合总数、前 4 景的影像 ID 及其云量、最终合成的影像对象、以及前 3 景的 B4/B8/NDVI 波段信息方便调试和验证。三、运行结果展示在 GEE 编辑器code.earthengine.google.com中执行上述代码后我们会看到地图中心显示维也纳及周边区域的真彩色合成影像色彩自然无明显云层遮挡。Console 输出影像总数例如 12 景具体数量取决于实际过境次数。前几景的 ID 和云量百分比可以看到云量最低的可能只有 2.3%。合成影像对象包含所有波段包括我们添加的 NDVI。前 3 景的 B4、B8 和 NDVI 数值数组。四、代码改进与思考云掩膜处理虽然我们过滤了整体云量 20% 的影像但每景内部仍有局部云。更严谨的做法是使用QA60波段生成云掩膜并替换为无云像素但代码中未做适合快速预览。时序合成如果希望得到整个夏季的“中值合成”或“最大值合成”如 NDVI 最大值可以改用median()或max()替代mosaic()这样能更好地消除残留云的影响。缩放因子显示真彩色时官方建议乘以 0.0001 转化为反射率0~1但代码直接用了 0~3000 的拉伸这在实际显示中也能获得不错的效果但数值意义需注意。元数据打印print语句在 GEE 中会触发计算如果集合很大建议先用limit()限制数量。五、总结通过这段简洁的 GEE 代码我们完成了从数据获取、预处理、指数计算到可视化的一整套遥感分析流程。无需下载任何数据仅需几行代码即可获得维也纳夏季高质量无云影像并附带植被指数信息。这体现了 GEE 在大地理空间数据处理中的巨大优势。如果你也对遥感数据处理感兴趣不妨在 GEE 平台上修改坐标、时间或筛选条件探索你所在城市或研究区域的影像变化。下一步我们可以尝试时间序列分析、土地覆盖分类或变化检测敬请期待附完整代码可直接复制到 GEE 编辑器运行// Earth Engine image collection example with a map layer.varviennaee.Geometry.Point([16.3738,48.2082]);varsentinel2ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED).filterBounds(vienna).filterDate(2024-06-01,2024-09-01).filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE,20)).map(function(image){varndviimage.normalizedDifference([B8,B4]).rename(NDVI);returnimage.addBands(ndvi).set(summary,summer clear-sky candidate);}).sort(CLOUDY_PIXEL_PERCENTAGE);varfirstClearImagesentinel2.mosaic();Map.addLayer(firstClearImage.select([B4,B3,B2]),{bands:[B4,B3,B2],min:0,max:3000},Sentinel-2 true color,true,0.7);print(Loaded collection,COPERNICUS/S2_SR_HARMONIZED);print(Filtered image count,sentinel2.size());print(Scene ids,sentinel2.limit(4).aggregate_array(system:index));print(Cloud percentages,sentinel2.limit(4).aggregate_array(CLOUDY_PIXEL_PERCENTAGE));print(First clear image,firstClearImage);print(Selected bands from first 3,sentinel2.limit(3).select([B4,B8,NDVI]));如果你有任何问题或希望看到更多 EE studio案例欢迎在评论区留言交流
用 Earth Engine 绘制维也纳夏季无云影像哨兵2数据筛选与 NDVI 计算实战本文带你一步步解析一段 Earth Engine 代码从指定区域、筛选影像、去云、计算植被指数到可视化最终获取维也纳上空 2024 年夏季最晴朗的哨兵2合成影像。一、为什么要写这段代码遥感数据处理往往面临两个痛点数据量大哨兵2每天全球产生数 TB 数据和云层遮挡光学影像受天气影响严重。Google Earth EngineGEE提供了云端并行计算能力让我们能在浏览器中快速完成大规模影像的筛选、预处理和分析。今天我们要实现的目的是锁定维也纳经纬度 16.3738°E, 48.2082°N周边区域获取 2024 年 6 月 1 日至 9 月 1 日期间所有哨兵2 Level-2A 地表反射率影像剔除云量超过 20% 的影像为每景影像计算NDVI归一化植被指数并附上标签按云量从低到高排序最终合成一张“最清晰”的 mosaic 影像在交互地图上展示真彩色合成并打印出一些关键元数据。二、代码分步拆解1. 定义研究区域varviennaee.Geometry.Point([16.3738,48.2082]);我们创建一个点几何对象代表维也纳市中心。这个点将作为空间筛选的基准。2. 加载哨兵2影像集合varsentinel2ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED).filterBounds(vienna).filterDate(2024-06-01,2024-09-01).filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE,20))使用官方提供的COPERNICUS/S2_SR_HARMONIZED数据集该数据集已进行辐射定标和大气校正且将 Sentinel-2A 与 2B 的波段进行了统一“HARMONIZED”。空间筛选只保留覆盖维也纳点的影像。时间筛选限定 2024 年夏季北半球植被生长旺季。云量筛选CLOUDY_PIXEL_PERCENTAGE是每景影像自带的元数据字段这里要求小于 20%避免云层干扰后续分析。3. 计算 NDVI 并添加为波段.map(function(image){varndviimage.normalizedDifference([B8,B4]).rename(NDVI);returnimage.addBands(ndvi).set(summary,summer clear-sky candidate);})对集合中的每一景影像使用normalizedDifference方法计算 (B8 - B4) / (B8 B4)即 NDVI。其中 B8 为近红外波段NIRB4 为红波段。将新生成的 NDVI 波段添加到原影像中并设置一个自定义属性summary方便后续识别这批影像的来源。4. 排序并合成.sort(CLOUDY_PIXEL_PERCENTAGE);varfirstClearImagesentinel2.mosaic();按照云量百分比升序排列云量越少越靠前。调用mosaic()将集合中的所有影像合成一张。由于云量已排序mosaic 会按照集合顺序叠加最后显示在最上层的是云量最少的像素但实际 mosaic 的规则是后加入的覆盖先加入的所以排序后低云影像会覆盖高云影像从而得到整体最清晰的合成图。这里严格来说mosaic 并不一定是“第一张清晰影像”而是所有影像的镶嵌结果但结合排序效果相当于优先使用低云像素。5. 地图可视化Map.addLayer(firstClearImage.select([B4,B3,B2]),{bands:[B4,B3,B2],min:0,max:3000},Sentinel-2 true color,true,0.7);选择红B4、绿B3、蓝B2三个波段合成真彩色影像。拉伸范围设为 03000对应地表反射率缩放因子原始值为整数一般乘以 0.0001 得到反射率但这里直接显示。图层透明度设为 0.7方便与底图叠加观察。6. 打印调试信息print(Loaded collection,COPERNICUS/S2_SR_HARMONIZED);print(Filtered image count,sentinel2.size());print(Scene ids,sentinel2.limit(4).aggregate_array(system:index));print(Cloud percentages,sentinel2.limit(4).aggregate_array(CLOUDY_PIXEL_PERCENTAGE));print(First clear image,firstClearImage);print(Selected bands from first 3,sentinel2.limit(3).select([B4,B8,NDVI]));输出集合总数、前 4 景的影像 ID 及其云量、最终合成的影像对象、以及前 3 景的 B4/B8/NDVI 波段信息方便调试和验证。三、运行结果展示在 GEE 编辑器code.earthengine.google.com中执行上述代码后我们会看到地图中心显示维也纳及周边区域的真彩色合成影像色彩自然无明显云层遮挡。Console 输出影像总数例如 12 景具体数量取决于实际过境次数。前几景的 ID 和云量百分比可以看到云量最低的可能只有 2.3%。合成影像对象包含所有波段包括我们添加的 NDVI。前 3 景的 B4、B8 和 NDVI 数值数组。四、代码改进与思考云掩膜处理虽然我们过滤了整体云量 20% 的影像但每景内部仍有局部云。更严谨的做法是使用QA60波段生成云掩膜并替换为无云像素但代码中未做适合快速预览。时序合成如果希望得到整个夏季的“中值合成”或“最大值合成”如 NDVI 最大值可以改用median()或max()替代mosaic()这样能更好地消除残留云的影响。缩放因子显示真彩色时官方建议乘以 0.0001 转化为反射率0~1但代码直接用了 0~3000 的拉伸这在实际显示中也能获得不错的效果但数值意义需注意。元数据打印print语句在 GEE 中会触发计算如果集合很大建议先用limit()限制数量。五、总结通过这段简洁的 GEE 代码我们完成了从数据获取、预处理、指数计算到可视化的一整套遥感分析流程。无需下载任何数据仅需几行代码即可获得维也纳夏季高质量无云影像并附带植被指数信息。这体现了 GEE 在大地理空间数据处理中的巨大优势。如果你也对遥感数据处理感兴趣不妨在 GEE 平台上修改坐标、时间或筛选条件探索你所在城市或研究区域的影像变化。下一步我们可以尝试时间序列分析、土地覆盖分类或变化检测敬请期待附完整代码可直接复制到 GEE 编辑器运行// Earth Engine image collection example with a map layer.varviennaee.Geometry.Point([16.3738,48.2082]);varsentinel2ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED).filterBounds(vienna).filterDate(2024-06-01,2024-09-01).filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE,20)).map(function(image){varndviimage.normalizedDifference([B8,B4]).rename(NDVI);returnimage.addBands(ndvi).set(summary,summer clear-sky candidate);}).sort(CLOUDY_PIXEL_PERCENTAGE);varfirstClearImagesentinel2.mosaic();Map.addLayer(firstClearImage.select([B4,B3,B2]),{bands:[B4,B3,B2],min:0,max:3000},Sentinel-2 true color,true,0.7);print(Loaded collection,COPERNICUS/S2_SR_HARMONIZED);print(Filtered image count,sentinel2.size());print(Scene ids,sentinel2.limit(4).aggregate_array(system:index));print(Cloud percentages,sentinel2.limit(4).aggregate_array(CLOUDY_PIXEL_PERCENTAGE));print(First clear image,firstClearImage);print(Selected bands from first 3,sentinel2.limit(3).select([B4,B8,NDVI]));如果你有任何问题或希望看到更多 EE studio案例欢迎在评论区留言交流
