## 1. 从零开始：为什么我们需要水体提取，以及NDWI和MNDWI是什么？ 大家好，我是老张，在遥感圈子里摸爬滚打了十几年，用Python处理过的卫星影像数据少说也有几个TB了。今天想和大家聊聊一个非常经典且实用的技术话题：如何用Python从卫星影像里，把水体（比如河流、湖泊、水库）给“抠”出来。这听起来是不是有点像Photoshop里的“魔棒”工具？没错，原理上还真有点类似，只不过我们是在代码的世界里，用数学公式和波段运算来实现的。 你可能想问，这事儿有啥用？用处可太大了。比如，监测洪涝灾害时，我们需要快速知道洪水淹了多大范围；评估一个地区的水资源状况，需要知道水库、湖泊的面积变化；甚至在城市规划里，也需要搞清楚城市水系的分布。传统的人工目视解译，效率低不说，还容易出错。而利用遥感技术，特别是像Sentinel-2、Landsat这些免费又高频的卫星数据，我们可以自动化、批量化地完成这项工作，又快又准。 那么，怎么“抠”呢？这就引出了我们今天要深入对比的两位“主角”：**NDWI**和**MNDWI**。它们都是“归一化差异水体指数”，名字听起来有点唬人，但原理其实很简单，就是利用水体在不同光谱波段（可以简单理解为不同颜色的光）下的反射特性不同，通过一个公式把水体给“凸显”出来。 **NDWI**，全称归一化差异水体指数，是McFeeters在1996年提出的。它的公式是 `(绿光波段 - 近红外波段) / (绿光波段 + 近红外波段)`。为什么这么算？因为水在绿光波段反射相对高一点（所以水看起来是绿的或蓝的），但在近红外波段吸收极强，反射率几乎为零。而植被和土壤在近红外反射很强。这么一减一除，水的值就会是正的（因为绿光>近红外），而植被和土壤的值会是负的（因为近红外>绿光），水体就被“亮”出来了。 听起来很完美对吧？但老方法总会遇到新问题。NDWI在城市区域就有点“水土不服”。因为城市里的建筑阴影、沥青路面，在近红外波段的反射也很低，容易被NDWI误判成水体。这就好比你的魔棒工具，不光选中了蓝色的湖，还把旁边灰色的水泥地也一块儿选中了。 于是，**MNDWI**，也就是改进的归一化差异水体指数，在2005年被徐涵秋教授提出来了。它把公式里的近红外波段换成了**短波红外波段**，变成了 `(绿光波段 - 短波红外波段) / (绿光波段 + 短波红外波段)`。这个改动非常巧妙，因为水体在短波红外的反射率比近红外还低，几乎可以忽略；而建筑物、土壤在短波红外的反射率却比在近红外要低得多。这样一来，水体和建筑物之间的“光谱反差”就被拉得更大了，在城市区域提取水体时，就能更好地避开建筑阴影的干扰。 所以，简单总结一下：**NDWI像是为广阔的自然水域（大江大湖）设计的“标准滤镜”**，而**MNDWI则是为城镇、城乡结合部这种复杂场景准备的“增强滤镜”**。但MNDWI也有个“硬伤”：它需要短波红外波段，而很多高分辨率商业卫星（比如高分二号）或者一些早期的数据，可能没有这个波段，这就限制了它的使用。 在接下来的内容里，我会手把手带你，用Python从读取一张卫星影像开始，一步步计算出这两种指数，直观地对比它们的效果，并分享我这些年积累下来的阈值选取“土办法”和后处理技巧，帮你避开我当年踩过的那些坑。无论你是遥感专业的学生，还是刚入行的数据分析师，相信这篇实战指南都能让你快速上手，做出漂亮又准确的水体提取结果。 ## 2. 实战第一步：搭建你的Python遥感处理环境与数据准备 工欲善其事，必先利其器。在开始写代码之前，我们得先把“厨房”收拾好。别担心，这个过程不复杂，我保证用最直白的方式带你搞定。 ### 2.1 安装必备的“武器库” 遥感数据处理离不开几个核心的Python库。我建议你直接使用Anaconda来管理环境，这样可以避免很多依赖冲突的烦心事。打开你的终端（Windows叫命令提示符或PowerShell，Mac/Linux叫Terminal），创建一个新的虚拟环境是个好习惯： ```bash conda create -n rs_water python=3.9 conda activate rs_water ``` 接下来，安装我们需要的库。这里我分两类，一类是数据处理和科学计算的“基石”，另一类是专门处理地理空间数据的“神器”。 **基石库**： ```bash pip install numpy matplotlib pillow jupyter ``` - `numpy`：不用多说，处理数组（我们的影像数据就是巨大的数字矩阵）的核心。 - `matplotlib` 和 `pillow`：用来显示和简单处理图像，让我们能直观地看到计算结果。 - `jupyter`：强烈推荐在Jupyter Notebook里跟着操作，可以边写代码边看图，交互体验极佳。 **地理空间神器库**： ```bash conda install -c conda-forge gdal rasterio geopandas ``` 这里我强烈推荐用`conda`从`conda-forge`这个频道安装。因为GDAL这个库（地理空间数据抽象库）是很多地理处理功能的底层依赖，用`pip`安装经常出幺蛾子，`conda`能更好地处理它的复杂依赖。`rasterio`是基于GDAL的、更Pythonic的栅格数据处理库，用起来比直接调GDAL的接口舒服多了。`geopandas`则是处理矢量数据（比如我们最后生成的水体边界多边形）的利器。 安装可能稍微需要点时间，喝杯咖啡等待一下。如果遇到网络问题，可以尝试更换国内的conda镜像源。 ### 2.2 理解并获取你的遥感数据 数据是分析的粮食。对于水体提取，我们通常需要多光谱影像，也就是包含多个波段（如红、绿、蓝、近红外等）的数据。免费的“午餐”有很多： - **Landsat系列**（美国）：历史悠久，数据丰富，包含短波红外波段，非常适合做MNDWI。可以从USGS EarthExplorer网站下载。 - **Sentinel-2系列**（欧盟）：更新、分辨率更高（10米/20米），重访周期短，也是绝佳选择。可以从欧空局的Copernicus Open Access Hub下载。 为了方便大家跟着练习，我准备了一小块Sentinel-2的示例数据。在实际项目中，你需要下载整景影像，然后根据你的研究区进行裁剪。下载下来的数据通常是一个包含多个波段文件的文件夹，或者是一个整合了所有波段的“堆叠”文件（比如GeoTIFF格式的`.tif`文件）。 在我们的代码里，我们将直接读取一个多波段的GeoTIFF文件。你需要知道每个波段对应的是什么。以Sentinel-2为例，通常波段顺序是：B1（海岸气溶胶）， B2（蓝）， B3（绿）， B4（红）， B8（近红外）， B11（短波红外1）等等。**计算NDWI我们需要绿光波段（B3）和近红外波段（B8）；计算MNDWI则需要绿光波段（B3）和短波红外波段（B11或B12）**。这个对应关系千万不能搞错，否则公式算出来就是一团糟。 假设你的数据文件叫 `my_sentinel2_image.tif`，把它放在一个你记得住的文件夹里，比如 `./data/`。接下来，我们就进入激动人心的代码环节。 ## 3. 核心算法实现：手把手编写NDWI与MNDWI计算代码 环境搭好，数据备齐，现在可以撸起袖子写代码了。我会把每一步都拆开讲清楚，你完全可以复制粘贴，但更希望你能理解每一行在做什么。 ### 3.1 读取影像数据：打开遥感图的“黑盒子” 卫星影像文件（GeoTIFF）不仅仅是一张图片，它里面还藏着地理坐标、投影信息等重要“元数据”。我们用 `rasterio` 来读取，它比原始的GDAL接口友好得多。 ```python import rasterio import numpy as np import matplotlib.pyplot as plt # 替换成你的影像文件路径 image_path = './data/my_sentinel2_image.tif' # 使用rasterio打开影像 with rasterio.open(image_path) as src: # 读取所有波段的数据，shape为 (波段数, 高度, 宽度) img_data = src.read() # 获取仿射变换参数（包含左上角坐标和像素分辨率） transform = src.transform # 获取投影信息 crs = src.crs # 获取图像尺寸 height, width = src.height, src.width print(f"影像形状（波段，高，宽）: {img_data.shape}") print(f"地理变换参数: {transform}") print(f"坐标系: {crs}") ``` 这段代码执行后，你就把影像数据加载到 `img_data` 这个三维的numpy数组里了。比如输出是 `(10, 10980, 10980)`，就表示有10个波段，每个波段是10980像素见方。 ### 3.2 可视化真彩色影像：先看看我们面对的是什么 在计算指数前，我们先看看原始影像长啥样，对地物有个直观认识。卫星数据原始值（DN值或反射率）范围很大，直接显示可能是一片黑或一片白，我们需要做个简单的拉伸。 ```python # 假设波段顺序中，第2、3、4个索引是蓝、绿、红波段（对应Sentinel-2的B2, B3, B4） # 注意：数组索引从0开始，所以是1，2，3 blue_band = img_data[1, :, :] green_band = img_data[2, :, :] red_band = img_data[3, :, :] # 将三个波段堆叠成RGB图像，形状变为（高，宽，3） rgb_image = np.stack([red_band, green_band, blue_band], axis=-1) # 定义一个简单的百分比拉伸函数，让图像看起来更清晰 def percent_stretch(band, lower_percent=2, upper_percent=98): lower = np.percentile(band, lower_percent) upper = np.percentile(band, upper_percent) band_stretched = np.clip(band, lower, upper) # 将超出范围的值截断 band_stretched = (band_stretched - lower) / (upper - lower) # 归一化到0-1 return band_stretched # 对每个波段进行拉伸 red_stretched = percent_stretch(rgb_image[:, :, 0]) green_stretched = percent_stretch(rgb_image[:, :, 1]) blue_stretched = percent_stretch(rgb_image[:, :, 2]) rgb_stretched = np.stack([red_stretched, green_stretched, blue_stretched], axis=-1) # 显示图像 plt.figure(figsize=(12, 10)) plt.imshow(rgb_stretched) plt.title('原始区域真彩色影像 (经过拉伸)') plt.axis('off') plt.show() ``` 现在，你应该能看到一幅清晰的卫星图了，河流、湖泊、城镇、农田都应该能分辨出来。 ### 3.3 计算NDWI：经典方法的代码实现 根据公式 NDWI = (Green - NIR) / (Green + NIR)，我们开始计算。这里要特别注意处理分母为零的情况，否则程序会报错。 ```python # 假设绿光波段是索引2（B3），近红外波段是索引7（B8） # 请根据你的数据实际情况调整索引！ green = img_data[2, :, :].astype(np.float32) # 转为浮点型以便计算 nir = img_data[7, :, :].astype(np.float32) # 初始化一个和影像一样大的数组，填充NaN（非数字）作为缺省值 ndwi = np.full_like(green, np.nan, dtype=np.float32) # 计算分母 denominator = green + nir # 只对分母不为零的像素进行计算 valid_mask = denominator != 0 ndwi[valid_mask] = (green[valid_mask] - nir[valid_mask]) / denominator[valid_mask] print(f"NDWI值范围: [{np.nanmin(ndwi):.3f}, {np.nanmax(ndwi):.3f}]") ``` 计算完成后，`ndwi` 数组里的值范围通常在 -1 到 1 之间。理论上，水体应该是正值，植被和裸土是负值。 ### 3.4 计算MNDWI：应对城镇场景的改进方法 MNDWI的计算过程几乎一模一样，只是把近红外波段换成短波红外波段。 ```python # 假设短波红外波段是索引10（B11） # 再次强调，根据你的数据调整索引！ swir = img_data[10, :, :].astype(np.float32) mndwi = np.full_like(green, np.nan, dtype=np.float32) denominator_m = green + swir valid_mask_m = denominator_m != 0 mndwi[valid_mask_m] = (green[valid_mask_m] - swir[valid_mask_m]) / denominator_m[valid_mask_m] print(f"MNDWI值范围: [{np.nanmin(mndwi):.3f}, {np.nanmax(mndwi):.3f}]") ``` 现在，我们手上已经有了 `ndwi` 和 `mndwi` 两个计算结果。光看数字没感觉，我们得把它们画出来，对比着看。 ## 4. 效果可视化与初步对比：谁才是真正的“水体侦探”？ 算完了不看看，等于白干。可视化不仅能帮我们理解结果，更是后续调整参数的关键依据。 ### 4.1 并排显示NDWI与MNDWI结果图 我们把原始影像、NDWI结果、MNDWI结果放在一起看。 ```python fig, axes = plt.subplots(1, 3, figsize=(18, 6)) # 第一幅：原始真彩色影像 axes[0].imshow(rgb_stretched) axes[0].set_title('原始真彩色影像') axes[0].axis('off') # 第二幅：NDWI结果 # 我们用‘viridis’色带，黄色/亮色代表高值（可能是水体），暗色代表低值 im_ndwi = axes[1].imshow(ndwi, cmap='viridis', vmin=-1, vmax=1) axes[1].set_title('NDWI 计算结果') axes[1].axis('off') plt.colorbar(im_ndwi, ax=axes[1], fraction=0.046, pad=0.04) # 添加颜色条 # 第三幅：MNDWI结果 im_mndwi = axes[2].imshow(mndwi, cmap='viridis', vmin=-1, vmax=1) axes[2].set_title('MNDWI 计算结果') axes[2].axis('off') plt.colorbar(im_mndwi, ax=axes[2], fraction=0.046, pad=0.04) plt.tight_layout() plt.show() ``` 这张对比图非常关键！你可能会立刻发现一些区别： - 在宽阔的自然河流或湖泊区域，两者可能都显示为亮色（高值），效果差不多。 - 但在城市区域，NDWI图上可能会出现一些零散的亮斑，这些很可能就是被误提取的建筑阴影或深色屋顶。而在MNDWI图上，这些误提的亮斑会少很多，水体（比如城市内的河道）的轮廓反而更清晰、更连续。 ### 4.2 利用直方图洞察数据分布，寻找阈值线索 光看图像还不够，我们需要更定量地分析。直方图能告诉我们指数值的分布情况，帮我们初步判断一个大概的阈值范围。 ```python fig, axes = plt.subplots(1, 2, figsize=(15, 5)) # NDWI直方图 axes[0].hist(ndwi[~np.isnan(ndwi)].flatten(), bins=200, range=(-1, 1), color='blue', alpha=0.7, edgecolor='black') axes[0].axvline(x=0, color='red', linestyle='--', linewidth=1, label='Zero Line') # 标记0值线 axes[0].set_xlabel('NDWI Value') axes[0].set_ylabel('Pixel Count') axes[0].set_title('NDWI Value Distribution Histogram') axes[0].legend() axes[0].grid(True, alpha=0.3) # MNDWI直方图 axes[1].hist(mndwi[~np.isnan(mndwi)].flatten(), bins=200, range=(-1, 1), color='green', alpha=0.7, edgecolor='black') axes[1].axvline(x=0, color='red', linestyle='--', linewidth=1, label='Zero Line') axes[1].set_xlabel('MNDWI Value') axes[1].set_ylabel('Pixel Count') axes[1].set_title('MNDWI Value Distribution Histogram') axes[1].legend() axes[1].grid(True, alpha=0.3) plt.tight_layout() plt.show() ``` 观察直方图，你会发现大部分像素都集中在负值区域（主要是植被和土壤），而在大于0的右侧，通常会有一个小的“鼓包”或者拖尾，这部分就很有可能是水体像素。**这个“鼓包”开始抬升的位置，比如0.1或0.2附近，就可以作为我们初始阈值的参考**。MNDWI的直方图可能比NDWI的“鼓包”更靠右，且与主峰分离得更开，这说明它区分水体和背景的能力可能更强。 ## 5. 阈值优化：找到区分水体与背景的“黄金分割线” 计算出的指数图是一个连续的值，我们需要一个阈值来做一个二值判断：大于阈值的是水体，小于等于的不是。阈值选多少？这是水体提取中最关键、也最需要经验的一步。没有放之四海而皆准的值，必须根据你的影像和区域来定。 ### 5.1 动态滑动观察法：我的“土办法” 我常用的一个笨办法但非常有效的方法，就是写个循环，动态生成一系列不同阈值下的二值图，像放幻灯片一样快速浏览，肉眼找到效果最好的那个。 ```python # 定义一个函数，用于生成和显示不同阈值下的二值掩膜 def preview_thresholds(index_array, index_name, thresholds): num_thresholds = len(thresholds) fig, axes = plt.subplots(1, num_thresholds, figsize=(num_thresholds*4, 4)) for i, thresh in enumerate(thresholds): # 创建二值掩膜：大于阈值为1（水体），否则为0 binary_mask = np.where(index_array > thresh, 1, 0) axes[i].imshow(binary_mask, cmap='gray') axes[i].set_title(f'{index_name}\nThreshold = {thresh:.2f}') axes[i].axis('off') plt.tight_layout() plt.show() # 为NDWI尝试一组阈值，比如从0到0.4，步长0.1 ndwi_thresholds_to_try = [0.0, 0.1, 0.2, 0.3, 0.4] preview_thresholds(ndwi, 'NDWI', ndwi_thresholds_to_try) # 为MNDWI尝试另一组阈值，因为它的值分布可能不同 mndwi_thresholds_to_try = [0.0, 0.1, 0.2, 0.3, 0.4] preview_thresholds(mndwi, 'MNDWI', mndwi_thresholds_to_try) ``` 通过快速浏览这组图片，你可以直观地看到： - **阈值太低（如0.0）**：很多非水体（如潮湿土壤、阴影）也被误认为是水体，提取范围过大，噪声极多。 - **阈值太高（如0.4）**：只有反射特性最强的核心水体部分被提取出来，一些边缘的、浑浊的或者有植被覆盖的水体（如湿地）会被漏掉，水体变得不连续。 - **合适的阈值（如0.1-0.3之间）**：能较好地平衡“误提”和“漏提”。你会发现NDWI和MNDWI的最佳阈值很可能不一样。通常，在城镇区域，MNDWI可以用比NDWI更高的阈值，依然能保持较好的提取效果，同时抑制更多噪声。 ### 5.2 基于统计的自动化阈值初选方法 除了肉眼观察，也可以借助一些统计方法获得初始阈值，再人工微调。这里介绍两种我常用的： **方法一：Otsu算法（大津法）** Otsu算法自动寻找一个阈值，使得分割后的前景（水体）和背景（非水体）的类间方差最大。对于双峰分布明显的直方图效果很好。 ```python from skimage.filters import threshold_otsu # 注意：Otsu需要输入一维数组，且忽略NaN ndwi_flat = ndwi[~np.isnan(ndwi)].flatten() if len(ndwi_flat) > 0: otsu_thresh_ndwi = threshold_otsu(ndwi_flat) print(f"Otsu算法推荐的NDWI阈值: {otsu_thresh_ndwi:.3f}") mndwi_flat = mndwi[~np.isnan(mndwi)].flatten() if len(mndwi_flat) > 0: otsu_thresh_mndwi = threshold_otsu(mndwi_flat) print(f"Otsu算法推荐的MNDWI阈值: {otsu_thresh_mndwi:.3f}") ``` **方法二：均值+标准差法** 这是一种更简单的经验方法。由于水体是图像中的“高亮”目标，我们可以认为水体像素是指数值大于“所有像素均值加上N倍标准差”的部分。N通常取1到3，需要试验。 ```python ndwi_mean, ndwi_std = np.nanmean(ndwi), np.nanstd(ndwi) mndwi_mean, mndwi_std = np.nanmean(mndwi), np.nanstd(mndwi) for n in [1, 1.5, 2]: ndwi_auto_thresh = ndwi_mean + n * ndwi_std mndwi_auto_thresh = mndwi_mean + n * mndwi_std print(f"均值+{n}倍标准差 -> NDWI阈值: {ndwi_auto_thresh:.3f}, MNDWI阈值: {mndwi_auto_thresh:.3f}") ``` 这些自动化方法给出的阈值可以作为你手动调整的起点，**但绝不能完全依赖**。一定要结合可视化结果，特别是与你目视判断的水体边界进行比对，来最终确定阈值。我个人的习惯是，先让自动化方法给个建议，然后用滑动观察法在建议值附近微调，比如尝试 `[建议值-0.05， 建议值， 建议值+0.05]` 这几个阈值，看哪个效果最顺眼。 ## 6. 后处理精修：让粗糙的“毛坯”变成精致的“成品” 即使用最优阈值得到了二值水体掩膜，它通常也还是个“毛坯房”，存在两个主要问题：1）内部可能有零星的非水体噪声点（椒盐噪声）；2）边界可能凹凸不平，存在很多基于单个像素的小图斑。这就需要后处理来精修。 ### 6.1 形态学滤波：平滑边界与填充空洞 形态学操作是图像处理中用于处理二值图形状的基本工具。我们主要用到“开运算”和“闭运算”。 - **开运算（先腐蚀后膨胀）**：可以消除细小的噪声点（孤立的白色点）。 - **闭运算（先膨胀后腐蚀）**：可以填充水体内部小的空洞，并连接邻近的细小水体。 ```python from scipy import ndimage # 假设我们确定了最终阈值，生成二值掩膜 final_ndwi_threshold = 0.15 binary_water_ndwi = np.where(ndwi > final_ndwi_threshold, 1, 0).astype(np.uint8) # 定义结构元素（核），这里用一个3x3的正方形 structure = np.ones((3, 3), dtype=np.uint8) # 开运算去除小噪声 opened = ndimage.binary_opening(binary_water_ndwi, structure=structure).astype(np.uint8) # 闭运算填充小空洞和平滑边界 closed = ndimage.binary_closing(opened, structure=structure).astype(np.uint8) # 可视化对比 fig, axes = plt.subplots(1, 3, figsize=(15, 5)) axes[0].imshow(binary_water_ndwi, cmap='gray') axes[0].set_title('原始二值掩膜') axes[0].axis('off') axes[1].imshow(opened, cmap='gray') axes[1].set_title('开运算后（去噪）') axes[1].axis('off') axes[2].imshow(closed, cmap='gray') axes[2].set_title('闭运算后（平滑）') axes[2].axis('off') plt.tight_layout() plt.show() ``` 你可以调整结构元素的大小（比如 `(5,5)`），来控制平滑和填充的力度。力度太大会改变水体的真实形状，需要谨慎。 ### 6.2 去除小图斑：基于连通域分析 有些噪声可能不是孤立的点，而是一小片区域。形态学滤波可能去不掉，或者会过度平滑。这时，基于连通域分析的方法就更精准。我们可以设定一个面积阈值，比如认为面积小于100个像素的连通区域都是噪声，将其剔除。 ```python from skimage import measure # 对后处理后的二值图像进行连通域标记 labeled_array, num_features = measure.label(closed, connectivity=2, return_num=True) print(f"找到 {num_features} 个连通区域") # 计算每个连通区域的属性，这里我们关心面积 props = measure.regionprops(labeled_array) areas = [prop.area for prop in props] # 设定面积阈值，小于此值的认为是小图斑（噪声） min_area_threshold = 100 # 创建一个和原图一样大的掩膜，初始为0 filtered_mask = np.zeros_like(closed, dtype=np.uint8) # 只保留面积大于阈值的区域 for prop in props: if prop.area >= min_area_threshold: # 将该区域对应的坐标在filtered_mask中设为1 filtered_mask[labeled_array == prop.label] = 1 # 可视化 fig, axes = plt.subplots(1, 2, figsize=(10, 5)) axes[0].imshow(closed, cmap='gray') axes[0].set_title('形态学滤波后') axes[0].axis('off') axes[1].imshow(filtered_mask, cmap='gray') axes[1].set_title(f'去除面积<{min_area_threshold}像素的小图斑后') axes[1].axis('off') plt.tight_layout() plt.show() ``` 经过这几步处理，你的水体掩膜应该干净、平滑多了，更接近真实水体的连续形状。 ## 7. 成果输出与应用：从栅格到矢量，完成闭环 得到漂亮的水体掩膜后，我们还需要把它保存下来，并转换成更通用的格式，方便在GIS软件（如QGIS, ArcGIS）中查看或进行进一步的空间分析。 ### 7.1 保存栅格结果 我们可以用 `rasterio` 将二值掩膜保存为新的GeoTIFF文件，同时保留原始影像的地理信息。 ```python # 使用原始影像的元数据作为模板 with rasterio.open(image_path) as src: meta = src.meta.copy() # 更新元数据以匹配我们的二值掩膜（单波段，数据类型为uint8） meta.update({ 'count': 1, 'dtype': 'uint8', 'nodata': 0 # 将0视为背景（非水体） }) # 保存NDWI和MNDWI的最终掩膜 output_ndwi_path = './data/water_mask_ndwi_final.tif' with rasterio.open(output_ndwi_path, 'w', **meta) as dst: dst.write(filtered_mask, 1) # 将我们的掩膜数组写入第一个波段 print(f"NDWI水体掩膜已保存至: {output_ndwi_path}") # 同样的流程可以保存MNDWI的结果 # ... (计算MNDWI的二值掩膜并后处理) # output_mndwi_path = './data/water_mask_mndwi_final.tif' # ... ``` ### 7.2 栅格转矢量：获取水体多边形边界 在GIS分析中，矢量多边形比栅格像元更方便计算面积、进行叠加分析等。我们可以将栅格掩膜转换成矢量面文件（如Shapefile）。 ```python import geopandas as gpd from rasterio import features # 为我们的掩膜生成多边形 shapes = features.shapes(filtered_mask.astype(np.uint8), transform=transform) # shapes生成的是 (多边形几何体, 像素值) 的迭代器 # 我们只关心值为1（水体）的多边形 polygons = [] values = [] for geom, val in shapes: if val == 1: # 只提取水体部分 polygons.append(geom) values.append(val) # 创建GeoDataFrame gdf = gpd.GeoDataFrame({'value': values, 'geometry': polygons}, crs=crs) # 保存为Shapefile output_shapefile_path = './data/water_polygons.shp' gdf.to_file(output_shapefile_path) print(f"水体矢量多边形已保存至: {output_shapefile_path}") # 快速查看一下矢量结果 fig, ax = plt.subplots(figsize=(10, 10)) # 在背景上显示原始RGB影像 ax.imshow(rgb_stretched, extent=[0, width, height, 0]) # 注意坐标范围 # 叠加显示水体多边形 gdf.boundary.plot(ax=ax, edgecolor='red', linewidth=1.5, facecolor='none') ax.set_title('最终提取的水体边界（叠加在真彩色影像上）') ax.axis('off') plt.show() ``` 现在，你不仅有了可视化的结果，还有了可以用于定量分析（如计算总面积、监测面积变化）和制图输出的标准地理数据产品了。整个从数据到信息的自动化流程就走通了。回顾一下，我们对比了NDWI和MNDWI在不同场景下的表现，探索了阈值优化的实用技巧，并运用后处理技术提升了成果的可用性。这套流程稍加修改，就可以应用于不同时间、不同区域的影像，实现大范围、长时间序列的水体动态监测。