常州企业建站系统模板ps在线网站

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 osdef 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.profiletransform = src.transformcrs = src.crswith 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("未检测到淹没区域")returnmerged_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.areagdf_filtered = gdf[gdf['area'] > min_area]if gdf_filtered.empty:print("过滤后无有效淹没区域")returngdf_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.nanif 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 = Noneout_ds = Noneds1 = 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))