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

基于 GEE 利用 Landsat4、5、7、8、9 数据计算 MNDWI 指数实现 1990—2024 年研究区水体变化分析

article 2025/4/2 3:59:01

目录

1 前言

2 代码解析

2.1 数据加载

2.2 时间与区域过滤

2.3 云掩膜处理

2.4 影像集合处理

2.5 指数计算与河流提取

2.6 可视化与导出

3 完整代码

4 运行结果

1 前言

随着遥感技术的飞速发展，Google Earth Engine（GEE）作为一款强大的云端地理空间数据处理平台，为研究人员提供了前所未有的便利。通过GEE，我们可以轻松访问并处理Landsat系列卫星影像，广泛应用于土地使用变化监测、环境分析和水资源管理等领域。本文将详细解析一段基于Landsat 4至9号卫星影像的河流提取代码，带你深入理解如何利用GEE实现多时相河流分布的提取与可视化。

Landsat系列卫星由美国国家航空航天局（NASA）和美国地质调查局（USGS）联合发射，自1972年以来持续提供地球表面的高分辨率影像数据。本代码利用了Landsat 4、5、7、8和9的二级产品（C02/T1_L2），这些数据经过大气校正和地形校正，适用于定量分析。代码的目标是通过计算改进的归一化水体指数（MNDWI）和归一化建筑指数（NDBI），结合云掩膜处理，从多时相影像中提取河流分布并导出结果。

代码结构与功能概述：

  1. 数据加载：加载Landsat 4至9的影像集。
  2. 时间与区域过滤：定义年份列表并实现影像的时间和空间过滤。
  3. 云掩膜处理：针对不同传感器设计云检测算法。
  4. 指数计算与河流提取：计算MNDWI和NDBI，基于阈值提取河流。
  5. 结果可视化与导出：将结果添加到地图并导出为GeoTIFF文件。

2 代码解析

2.1 数据加载

var l4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2"),
    l5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2"),
    l7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2"),
    l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2"),
    l9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2");

这段代码定义了五个影像集变量，分别对应Landsat 4、5、7、8和9的二级产品。这些数据集包含表面反射率（SR）和质量评估（QA）波段，适合用于后续的云检测和指数计算。值得注意的是，Landsat 4和5使用TM传感器，Landsat 7使用ETM+传感器，而Landsat 8和9使用OLI/TIRS传感器，波段命名和特性有所不同，因此需要分别处理。

2.2 时间与区域过滤

var yearList = [1990, 1995, 2000, 2005, 2010, 2015, 2020, 2024];

function filterCol(col, roi, date) {
  return col.filterDate(date[0], date[1]).filterBounds(roi);
}
  • yearList 定义了分析的年份，从1990年至2024年，每五年一个间隔（2024年为最新年份）。
  • filterCol 是一个通用过滤函数，接收影像集（col）、感兴趣区域（roi）和时间范围（date），通过 filterDate 和 filterBounds 方法分别按时间和空间过滤影像。这种模块化设计提高了代码复用性。

2.3 云掩膜处理

由于云覆盖会干扰地表特征的识别，代码为不同传感器设计了两套云掩膜函数：

Landsat 4、5、7（TM/ETM+传感器）云掩膜函数：

function cloudMaskTm(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));
  
  return image.select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'], ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'])
    .updateMask(mask).multiply(0.0000275).add(-0.2);
}
  • QA_PIXEL 是质量评估波段，包含云、云阴影等信息。
  • bitwiseAnd 使用位运算检查特定位（1：云填充，3：云，4：云阴影），eq(0) 表示该位为0（无云）。
  • 波段重命名：将原始波段名称标准化为 'B2', 'B3', 'B4', 'B5', 'B6', 'B7'，便于后续计算。
  • 反射率转换：乘以0.0000275并加-0.2，将DN值转换为表面反射率。

Landsat 8、9（OLI传感器）云掩膜函数：

function cloudMaskOli(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 2).eq(0))
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));
  
  return image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'], ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'])
    .updateMask(mask).multiply(0.0000275).add(-0.2);
}

与 cloudMaskTm 类似，但增加了对位2（稀薄云）的检查，因为OLI传感器的QA波段定义有所不同。

2.4 影像集合处理

function landsat457(roi, date) {
  var col = filterCol(l4, roi, date).merge(filterCol(l5, roi, date));
  return col.map(cloudMaskTm).median().clip(roi);
}

function landsat89(roi, date) {
  var col = filterCol(l8, roi, date).merge(filterCol(l9, roi, date));
  return col.map(cloudMaskOli).median().clip(roi);
}
  • landsat457 和 landsat89 分别处理早期（Landsat 4、5）和晚期（Landsat 8、9）影像。
  • merge 将多个影像集合合并。
  • map 应用云掩膜函数，**median** 计算中值影像以减少噪声，**clip** 裁剪到指定区域。

2.5 指数计算与河流提取

yearList.forEach(function(year) {
  var start = ee.Date.fromYMD(year - 1, 1, 1);
  var end = ee.Date.fromYMD(year + 1, 12, 31);
  var date = [start, end];

  var landsat = year < 2014 ? landsat457 : landsat89;
  var image = landsat(roi, date);

  var mndwi = image.expression('(GREEN - SWIR) / (GREEN + SWIR)', {
    GREEN: image.select('B3'),
    SWIR: image.select('B6')
  }).rename('MNDWI');

  var ndbi = image.expression('(SWIR - NIR) / (SWIR + NIR)', {
    SWIR: image.select('B6'),
    NIR: image.select('B5')
  }).rename('NDBI');

  var river = mndwi.gt(0.1).and(ndbi.lt(0)).selfMask().rename('river_' + year).toUint8();
  • 时间窗口：以目标年份为中心，前后各一年，确保数据充足。
  • 传感器选择：2014年为分界线（Landsat 8于2013年发射），早期使用 landsat457，晚期使用 landsat89。
  • MNDWI：改进归一化水体指数，利用绿波段（B3）和短波红外（SWIR，B6）增强水体特征。
  • NDBI：归一化建筑指数，利用SWIR（B6）和近红外（NIR，B5）区分建筑。
  • 河流提取：mndwi > 0.1 筛选水体，ndbi < 0 排除建筑，selfMask 只保留满足条件的像素。

2.6 可视化与导出

  Map.addLayer(river, {palette: 'blue'}, 'River_' + year, false);

  Export.image.toDrive({
    image: river,
    description: 'River_' + year,
    fileNamePrefix: 'River_' + year,
    scale: 30,
    region: roi,
    fileFormat: 'GeoTIFF'
  });
});
  • Map.addLayer：将河流图层添加到GEE地图，蓝色调显示，默认不显示（false）。
  • Export.image.toDrive：将结果导出到Google Drive，分辨率为30米，格式为GeoTIFF。

3 完整代码

var l4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2"),
    l5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2"),
    l7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2"),
    l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2"),
    l9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2");

var yearList = [1990, 1995, 2000, 2005, 2010, 2015, 2020,2024];

function filterCol(col, roi, date) {
  return col.filterDate(date[0], date[1]).filterBounds(roi);
}

function landsat457(roi, date) {
  var col = filterCol(l4, roi, date).merge(filterCol(l5, roi, date));
  return col.map(cloudMaskTm).median().clip(roi);
}

function landsat89(roi, date) {
  var col = filterCol(l8, roi, date).merge(filterCol(l9, roi, date));
  return col.map(cloudMaskOli).median().clip(roi);
}

function cloudMaskTm(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));
  
  return image.select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'], ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'])
    .updateMask(mask).multiply(0.0000275).add(-0.2);
}

function cloudMaskOli(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 2).eq(0))
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));
  
  return image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'], ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'])
    .updateMask(mask).multiply(0.0000275).add(-0.2);
}

yearList.forEach(function(year) {
  var start = ee.Date.fromYMD(year - 1, 1, 1);
  var end = ee.Date.fromYMD(year + 1, 12, 31);
  var date = [start, end];

  var landsat = year < 2014 ? landsat457 : landsat89;
  var image = landsat(roi, date);

  var mndwi = image.expression('(GREEN - SWIR) / (GREEN + SWIR)', {
    GREEN: image.select('B3'),
    SWIR: image.select('B6')
  }).rename('MNDWI');

  var ndbi = image.expression('(SWIR - NIR) / (SWIR + NIR)', {
    SWIR: image.select('B6'),
    NIR: image.select('B5')
  }).rename('NDBI');

  var river = mndwi.gt(0.1).and(ndbi.lt(0)).selfMask().rename('river_' + year).toUint8();
  Map.addLayer(river, {palette: 'blue'}, 'River_' + year, false);

  Export.image.toDrive({
    image: river,
    description: 'River_' + year,
    fileNamePrefix: 'River_' + year,
    scale: 30,
    region: roi,
    fileFormat: 'GeoTIFF'
  });
});

4 运行结果

点击RUN即可下载数据

http://www.kler.cn/a/623259.html

相关文章：

  • Spring Boot 3.4.3 基于 Caffeine 实现本地缓存
  • Linux基础指令（一）
  • golang 的io与os包中的常用方法
  • VITA 模型解读，实时交互式多模态大模型的 pioneering 之作
  • HarmonyOs学习 实验六：tabs标签与Swiper轮播图页面设计
  • 2023年3月全国计算机等级考试真题（二级C语言）
  • 【算法】并查集基础讲解
  • TCP协议与wireshark抓包分析
  • 现代优雅杂志海报徽标设计手写英文字体安装包 Attomes – Brush Handwritten Font
  • 【Prompt实战】邮件意图分类助手
  • git | 版本切换的相关指令
  • 深度学习入门（二）：从感知机到神经网络
  • （三）物理设备
  • 创作领域“＜em ＞一＜/em＞＜em＞分＜/em＞＜em＞快＜/em＞＜em＞3＜/em＞＜em＞官＜/em＞＜em＞网
  • 关于参加CSP-J/S认证需符合年龄条件的公告（2025年起）
  • 漏洞挖掘---灵当CRM客户管理系统getOrderList SQL注入漏洞
  • 保存预测图像时出现的文件名错误
  • Kubernetes 存储
  • NQA 网络质量分析协议
  • uniapp uni-swipe-action滑动内容排版改造