基于Sentinel-1A GRD洪涝淹没范围提取（SDWI阈值法和OSTU自动阈值法）

0 前言

两幅灾前和灾后的遥感影像经过SARscape配准、滤波、辐射定标预处理之后，使用GDAL库分别使用SDWI阈值法和OSTU自动阈值法提取洪涝淹没范围

1 ENVI 5.6和SARscape5.6安装

通过网盘分享的文件：ENVI5(1).6
链接: https://pan.baidu.com/s/1mKcEkC3rJDxs4p_RuT64Gg?pwd=qwea 提取码: qwea （包含了ENVI5.6和SARscape 5.6）

ENVI 5.6软件安装教程 (安装ENVI5.6)

envi5.6+SARscape560安装(CSDN_20240623)_sarscape5.6安装教程-CSDN博客（安装SARscape）

安装很简单，对着教程一步步来，是不会出错的。

2 SARscape 对遥感影像进行预处理

基于雷达影像的洪水监测技术方法详解

我是对照着上面的文章进行预处理，但是在处理过程中发现了几个问题

2.1 没有下载精密轨道文件，处理后查看直方图发现值都是负数

下载精密轨道文件方法：3种方法下载Sentinel-1精密轨道数据_精密轨道数据下载-CSDN博客

假如下载的遥感影像日期为8月17日，哪么你要下载的轨道数据名称为

S1A_OPER_AUX_POEORB_OPOD_20230906T080631_V20230816T225942_20230818T005942.EOF，最后两个时间节点分别为成像前一天和成像后一天，最前面的哪个时间是生成轨道文件生成时间，这个不用管。

导入轨道文件

这个不能和原始SAR数据一起导入，SARscape会在导入SAR数据的时候在你设置的目录中去找轨道数据

下载Sentinel-1精密轨道数据(Poeorb)及如何使用-CSDN博客

2.2 导入DEM数据进行配准和辐射定标

[SARscape] 将DEM数据导入到SARscape软件_sarscape dem ftp地址-CSDN博客

一定要导入DEM数据

3 SWDI和OSTU

import rasterio
from rasterio import features
import numpy as np
from skimage.filters import threshold_otsu
from skimage.morphology import remove_small_objects, binary_closing, square
import geopandas as gpd
from shapely.geometry import shape
from shapely.ops import unary_union
import os

def extract_flood_area(pre_vv_path, post_vv_path, output_shp,
                       min_pixel_size=50, min_area=10000, closing_kernel=5):
    # ========== 1. 数据加载与差异计算 ==========
    with rasterio.open(pre_vv_path) as src:
        pre_vv = src.read(1)
        profile = src.profile
        transform = src.transform
        crs = src.crs

    with rasterio.open(post_vv_path) as src:
        post_vv = src.read(1)

    # 单位检查（可选）
    if np.nanmin(pre_vv) < 0 or np.nanmin(post_vv) < 0:
        print("检测到分贝单位（dB），确保两期数据单位一致。")

    # 计算差异（单位需一致）
    delta_vv = post_vv - pre_vv

    # ========== 2. 大津法阈值分割 ==========
    valid_pixels = delta_vv[~np.isnan(delta_vv)]
    if len(valid_pixels) == 0:
        raise ValueError("差异图无效，无有效数据用于阈值计算")
    thresh = threshold_otsu(valid_pixels)
    print(f"自动计算阈值: {thresh:.2f}")  # 修正单位标注

    # 生成二值掩膜
    mask = (delta_vv < thresh).astype(bool)

    # ========== 3. 后处理优化 ==========
    mask_clean = remove_small_objects(mask, min_size=min_pixel_size)
    mask_clean = binary_closing(mask_clean, square(closing_kernel))

    # ========== 4. 矢量转换与保存 ==========
    polygons = []
    for geom, val in features.shapes(mask_clean.astype(np.uint8), transform=transform, connectivity=8):
        if val == 1:
            polygons.append(shape(geom))

    if not polygons:
        print("未检测到淹没区域")
        return

    merged_geoms = unary_union([geom for geom in polygons if geom.area > 0])
    merged_polygons = list(merged_geoms.geoms) if merged_geoms.geom_type == 'MultiPolygon' else [merged_geoms]

    gdf = gpd.GeoDataFrame(geometry=merged_polygons, crs=crs)
    gdf['area'] = gdf.geometry.area
    gdf_filtered = gdf[gdf['area'] > min_area]

    if gdf_filtered.empty:
        print("过滤后无有效淹没区域")
        return

    gdf_filtered.geometry = gdf_filtered.geometry.simplify(5.0)
    gdf_filtered.to_file(output_shp)
    print(f"成功提取 {len(gdf_filtered)} 个淹没区域，保存至：{output_shp}")

# 输入参数（示例）

pre_vv = r""
post_vv = r""
output_shp = r""

# 执行提取
extract_flood_area(
    pre_vv_path=pre_vv,
    post_vv_path=post_vv,
    output_shp=output_shp,
    min_pixel_size=300,
    min_area=50000,
    closing_kernel=3
)

双极化水体指数：

import numpy as np
import os
from osgeo import gdal

# 设置参数
threshold =  # 根据实际情况调整NDPWI阈值
output_path = r''  # 输出结果路径


def extract_flood(ndpwi1_path, ndpwi2_path, output_path, threshold):
    # 确保输出目录存在
    os.makedirs(os.path.dirname(output_path), exist_ok=True)

    # 打开724期NDPWI影像
    ds1 = gdal.Open(ndpwi1_path, gdal.GA_ReadOnly)
    if ds1 is None:
        raise FileNotFoundError(f"无法打开文件：{ndpwi1_path}")
    band1 = ds1.GetRasterBand(1)
    ndpwi1 = band1.ReadAsArray().astype(np.float32)
    nodata1 = band1.GetNoDataValue()

    # 打开817期NDPWI影像
    ds2 = gdal.Open(ndpwi2_path, gdal.GA_ReadOnly)
    if ds2 is None:
        raise FileNotFoundError(f"无法打开文件：{ndpwi2_path}")
    band2 = ds2.GetRasterBand(1)
    ndpwi2 = band2.ReadAsArray().astype(np.float32)
    nodata2 = band2.GetNoDataValue()

    # 处理NoData值
    if nodata1 is not None:
        ndpwi1[ndpwi1 == nodata1] = np.nan
    if nodata2 is not None:
        ndpwi2[ndpwi2 == nodata2] = np.nan

    # 生成水体掩膜
    water_724 = (ndpwi1 > threshold).astype(np.uint8)
    water_817 = (ndpwi2 > threshold).astype(np.uint8)

    # 提取新增淹没区域（817期为水且724期非水）
    flood_mask = np.where((water_817 == 1) & (water_724 == 0), 1, 0)

    # 处理NaN区域（设置为0）
    flood_mask = np.where(np.isnan(ndpwi1) | np.isnan(ndpwi2), 0, flood_mask)

    # 创建输出文件（关键修正部分）
    driver = gdal.GetDriverByName('GTiff')
    if driver is None:
        raise RuntimeError("GDAL GeoTIFF驱动未加载，请检查环境配置")

    # 确保删除已存在的文件
    if os.path.exists(output_path):
        driver.Delete(output_path)

    # 创建新的输出数据集（添加OVERWRITE选项）
    out_ds = driver.Create(
        output_path,
        ds1.RasterXSize,
        ds1.RasterYSize,
        1,
        gdal.GDT_Byte,
        options=['COMPRESS=LZW', 'OVERWRITE=YES']  # 添加压缩和覆盖选项
    )

    if out_ds is None:
        raise RuntimeError(f"无法创建输出文件：{output_path}，请检查路径权限或磁盘空间")

    # 写入地理参考和投影
    out_ds.SetGeoTransform(ds1.GetGeoTransform())
    out_ds.SetProjection(ds1.GetProjection())

    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(flood_mask)
    out_band.SetNoDataValue(0)
    out_band.FlushCache()

    # 显式释放资源
    out_band = None
    out_ds = None
    ds1 = ds2 = None


# 执行提取
try:
    extract_flood(
        ndpwi1_path=r'',  # 注意检查路径拼写是否正确
        ndpwi2_path=r'',
        output_path=output_path,
        threshold=
    )
    print("洪涝淹没范围提取完成！输出文件：", output_path)
except Exception as e:
    print("处理失败:", str(e))