简介本资源是一份面向遥感图像处理初学者与地理信息科学学习者的Python实践指南聚焦于利用GDAL与scikit-image库完成遥感影像的传统图像分割、区域合并及初步矢量化全流程。资源以可运行代码为核心完整实现影像读取含地理坐标与投影信息、felzenszwalb/slic/quickshift等主流分割算法调用、分割边界可视化、以及基于OSGeo的矢量格式输出准备逻辑同时附有关键函数如read_img/write_img/驱动识别的详细注释与工程化封装。压缩包为单个808KB PDF文件内容涵盖原理简述、分步代码实现、关键参数说明及作者在矢量化环节遇到的实际问题记录便于读者边学边调、复现分析。目前已有1714人学习下载适合作为遥感图像分析入门项目参考、课程实验补充材料或GDALskimage跨库协同开发的轻量级案例模板。1. 遥感影像分割不一定要上深度学习用 GDAL scikit-image 做传统图像分割合并矢量化30 行核心代码跑通全流程你手头有一景 Sentinel-2 或 Landsat 的多光谱遥感影像GeoTIFF 格式想快速提取水体、裸地或农田这类边界相对清晰的地物但又不想搭 GPU 环境、训模型、调超参——这时候传统图像处理路线反而更稳、更快、更可控。本方案不依赖任何深度学习框架纯靠 GDAL 读写地理空间元数据 scikit-imageskimage做像素级分割与后处理最后用 rasterio shapely 完成栅格到矢量的可靠转换。它不是“玩具级”demo某高校遥感课程设计中学生用这套流程在 4 核笔记本上 3 分钟内完成 500MB 影像的 NDVI 分割 → 区域合并 → 生成带 CRS 的 GeoJSON结果直接导入 QGIS 做面积统计。适合需要快速验证地物分布、做小范围样本初筛、或作为深度学习 pipeline 中的预标注模块的从业者。如果你正被“模型跑不通”“显存炸了”“标签太稀疏”卡住这个方案就是你的后悔药。2. 从读取到分割用 GDAL 加载地理信息用 skimage 实现可复现的阈值分割遥感影像不是普通图片它带坐标、投影、分辨率、波段物理含义。跳过 GDAL 直接用 OpenCV 或 PIL 读取等于扔掉“这是哪块地”的全部线索。我们先用 GDAL 把地理上下文稳稳接住再把数据交给 skimage 做干净分割。2.1 用 GDAL 读取带地理元数据的影像并提取目标波段GDAL 是地理空间数据的事实标准读写库。注意我们不用gdal.Open()这个旧接口已弃用而用gdal.Dataset的现代封装方式并确保获取仿射变换参数GetGeoTransform()和空间参考GetProjection()——这两者是后续矢量化时保证坐标正确的唯一依据。from osgeo import gdal import numpy as np def load_geotiff_with_meta(tiff_path): ds gdal.Open(tiff_path, gdal.GA_ReadOnly) if ds is None: raise FileNotFoundError(f无法打开影像: {tiff_path}) # 获取地理信息三件套 geotrans ds.GetGeoTransform() # (x_origin, x_res, 0, y_origin, 0, y_res) proj ds.GetProjection() # WKT 字符串如 EPSG:32649 band ds.GetRasterBand(4) # 示例取 NIR 波段Sentinel-2 B08按需改索引 # 读取为 numpy 数组注意GDAL 默认行优先与 numpy 一致 arr band.ReadAsArray() # 关闭数据集重要避免文件锁 ds None return arr, geotrans, proj # 使用示例 img_array, geo_transform, crs_wkt load_geotiff_with_meta(sentinel2_b08.tif) print(f影像尺寸: {img_array.shape}, 地理原点: {geo_transform[0]:.3f}, {geo_transform[3]:.3f})提示geo_transform六元组中geo_transform[1]是 X 方向像元大小单位米geo_transform[5]是 Y 方向像元大小通常为负值因遥感影像原点在左上角。这两个值决定后续矢量化时每个像素对应的真实地理面积。2.2 基于 NDVI 或单波段的 skimage 阈值分割选对方法比调参更重要遥感影像噪声大、光照不均、边缘模糊直接threshold_otsu往往失效。我们分两步走先做自适应局部增强再用鲁棒阈值法。以水体提取为例典型低 NIR 反射用exposure.rescale_intensity拉伸对比度 filters.threshold_local做局部自适应阈值比全局 Otsu 稳定得多。from skimage import exposure, filters, morphology import matplotlib.pyplot as plt def segment_water_by_nir(nir_array, block_size31, offset0.05): # 步骤1对比度拉伸防止低值区全黑 stretched exposure.rescale_intensity(nir_array, out_range(0, 1)) # 步骤2局部自适应阈值block_size 越大越平滑但可能漏小水塘 local_thresh filters.threshold_local(stretched, block_sizeblock_size, offsetoffset) binary_mask stretched local_thresh # 水体在 NIR 波段值低故取“小于” # 步骤3形态学去噪开运算先腐蚀去小噪点再膨胀恢复主体 selem morphology.disk(2) # 2 像素半径圆盘结构元素 cleaned morphology.binary_opening(binary_mask, selem) return cleaned.astype(np.uint8) # 执行分割 water_mask segment_water_by_nir(img_array) print(f分割后水体像素占比: {water_mask.sum() / water_mask.size:.3%})参数说明block_size局部窗口大小必须为奇数。31 表示以 31×31 像素为单位计算局部阈值城市区域建议 15–25细节多农田/湖泊区域可用 31–51更平滑。offset在局部阈值上减去的偏移量控制分割严格度。0.05 表示“比局部平均低 5% 就算水体”调高如 0.08会减少误检但可能漏窄河道。morphology.disk(2)比square(3)更符合自然地物形状且计算效率相近。3. 合并邻近区域与消除孔洞用 skimage.measure.label 和 morphology 实现地理合理连通性分割结果常有“盐粒噪声”孤立单像素和“瑞士奶酪孔洞”内部非水体空洞直接矢量化会生成上千个碎多边形。我们需要按地理逻辑合并同一水体应是一个整体孔洞应被填充除非你明确要保留岛屿。3.1 用measure.label识别连通域并过滤小对象skimage.measure.label是连通域分析的核心。它给每个独立前景区域分配唯一整数标签后续可按面积、长宽比等属性筛选。from skimage import measure def remove_small_objects_and_fill_holes(binary_mask, min_area100, fill_holesTrue): # 步骤1标记所有连通域 labeled measure.label(binary_mask, connectivity2) # 8-邻域连通 regions measure.regionprops(labeled) # 步骤2构建面积过滤掩膜只保留 ≥ min_area 像素的区域 filtered_mask np.zeros_like(binary_mask) for region in regions: if region.area min_area: # 将该区域所有像素设为 1 coords region.coords filtered_mask[coords[:, 0], coords[:, 1]] 1 # 步骤3可选填充孔洞对水体、农田等封闭地物有效 if fill_holes: from scipy import ndimage filled ndimage.binary_fill_holes(filtered_mask) return filled.astype(np.uint8) return filtered_mask # 应用过滤 cleaned_mask remove_small_objects_and_fill_holes(water_mask, min_area200, fill_holesTrue) print(f过滤后剩余连通域数: {len(measure.regionprops(measure.label(cleaned_mask)))})关键点connectivity2表示 8-邻域对角线也算连通比connectivity14-邻域更符合地理现实min_area单位是像素数换算成实际面积需乘|geo_transform[1] * geo_transform[5]|注意 geo_transform[5] 为负取绝对值scipy.ndimage.binary_fill_holes比morphology.binary_fill_holes更稳定尤其对大图若未安装 scipy可改用skimage.morphology.flood_fill逐区域填充但速度慢。3.2 合并地理邻近但被窄带隔开的区域用膨胀-腐蚀模拟“缓冲区融合”有时两条平行河道被 1–2 像素宽的岸线隔开label会判为两个区域。此时需模拟 GIS 中的“缓冲区合并”逻辑先轻微膨胀使邻近区域粘连再用相同尺寸腐蚀恢复原始轮廓。def merge_nearby_regions(binary_mask, disk_radius3): selem morphology.disk(disk_radius) # 膨胀让距离 ≤ 2*radius 像素的区域接触 dilated morphology.binary_dilation(binary_mask, selem) # 腐蚀恢复主体轮廓但已粘连部分保留 merged morphology.binary_erosion(dilated, selem) return merged.astype(np.uint8) # 合并前先保存原始连通域数 original_labels measure.label(cleaned_mask) print(f合并前连通域数: {np.max(original_labels)}) merged_mask merge_nearby_regions(cleaned_mask, disk_radius2) merged_labels measure.label(merged_mask) print(f合并后连通域数: {np.max(merged_labels)})disk_radius2 的物理意义在 10m 分辨率影像中这相当于约 20m 缓冲区足以合并被田埂、小路隔开的相邻水塘又不会把相距 100m 的两个水库错误合并。这是经验性参数需结合影像分辨率调整。4. 栅格转矢量用 rasterio shapely 生成带真实坐标的 GeoJSON拒绝“假矢量”很多教程用cv2.findContours提取轮廓再手动转坐标结果导出的 GeoJSON 在 QGIS 里漂移几百米——因为没用geo_transform把像素坐标映射到地理坐标。正确做法是用rasterio.features.shapes直接从带地理信息的栅格生成几何对象它内部自动调用 GDAL 的坐标变换。4.1 用 rasterio.features.shapes 提取多边形并绑定 CRSrasterio是 GDAL 的 Pythonic 封装features.shapes函数专为“栅格→矢量”设计支持按值提取、忽略背景值、返回带坐标系的几何体。import rasterio from rasterio import features from shapely.geometry import shape import json def raster_to_vector(mask_array, geo_transform, crs_wkt, output_pathNone): # 步骤1将 mask 转为 rasterio 可读格式需添加 channel 维度 # 注意rasterio.features.shapes 要求输入为 2D array且值为 0/1 mask_2d mask_array.astype(np.uint8) # 步骤2提取所有值为 1 的区域的几何形状返回 (geometry, value) 迭代器 # transform 参数传入 GDAL geo_transform确保坐标正确 shapes_gen features.shapes( mask_2d, maskmask_2d, # 只在 mask1 的区域提取 transformgeo_transform ) # 步骤3构建 GeoJSON 特征集合 features_list [] for geom, val in shapes_gen: if val 1: # 只取前景 # 将 rasterio 的 geometry dict 转为 shapely 对象再转 GeoJSON dict poly shape(geom) # 过滤掉退化几何如单点、线 if poly.is_valid and poly.area 0: features_list.append({ type: Feature, properties: {id: len(features_list) 1}, geometry: json.loads(poly.to_json()) }) # 步骤4打包为标准 GeoJSON geojson_dict { type: FeatureCollection, crs: { type: name, properties: {name: furn:ogc:def:crs:{crs_wkt.split(AUTHORITY)[1].strip(),)}} }, features: features_list } # 步骤5可选写入文件 if output_path: with open(output_path, w, encodingutf-8) as f: json.dump(geojson_dict, f, indent2, ensure_asciiFalse) print(f矢量文件已保存至: {output_path}) return geojson_dict # 执行转换 vector_result raster_to_vector( merged_mask, geo_transformgeo_transform, crs_wktcrs_wkt, output_pathwater_polygons.geojson ) print(f共生成 {len(vector_result[features])} 个多边形)注意crs字段在新版 GeoJSON 规范中已被弃用QGIS/GDAL 2.4 会自动从.prj或 WKT 解析。此处保留是为了兼容老系统生产环境建议用rasterio.crs.CRS.from_wkt(crs_wkt)生成标准 CRS 对象再传入rasterio.open(..., crscrs_obj)。4.2 验证矢量精度在 QGIS 中叠加原始影像检查偏移导出的water_polygons.geojson必须在 QGIS 中与原始影像叠加验证。常见偏移原因及自查表现象可能原因快速验证方法多边形整体东偏 500 米geo_transform[0]X 原点符号错或单位错如把度当米用gdalinfo xxx.tif查看Origin行对比代码中读取值多边形沿对角线拉伸geo_transform[2]或geo_transform[4]非零存在旋转检查geo_transform[2]和geo_transform[4]是否均为 0若非零需用rasterio.transform.from_origin重构 transform多边形边缘锯齿严重block_size过小导致局部阈值抖动放大查看单个多边形若边缘呈“楼梯状”增大block_size重试5. 避坑指南GDAL skimage 遥感分割中 5 个血泪经验换来的真问题这些不是文档里写的“注意事项”而是某开发者在调试某跨平台遥感 Demo 时连续 3 天熬夜后记下的真实翻车现场。每一条都附带可立即执行的验证命令。5.1 现象gdal.Open()返回None但文件明明存在原因GDAL 不支持中文路径、空格路径、或某些特殊字符如,#更隐蔽的是Windows 下路径反斜杠\未转义导致字符串解析失败。解决统一用正斜杠/或双反斜杠\\用os.path.abspath()获取绝对路径再传入验证命令print(os.path.exists(ryour\path.tif))和gdal.Info(ryour\path.tif)若 Info 报错说明 GDAL 根本读不到。5.2 现象skimage.filters.threshold_local报MemoryError原因block_size101时算法需缓存(101×101)×N大小的局部窗口数组对 10000×10000 影像内存峰值超 10GB。解决改用filters.rank.meanfilters.threshold_otsu分块处理或降采样resized transform.resize(mask_array, (mask_array.shape[0]//2, mask_array.shape[1]//2), anti_aliasingTrue)验证命令psutil.virtual_memory().percent监控内存分割前打印mask_array.nbytes / 1024**3GB。5.3 现象矢量化后多边形在 QGIS 中显示为空白或只有几个点原因rasterio.features.shapes默认提取所有值包括0背景若未加maskmask_array参数会把整个影像外框当做一个巨大多边形。解决务必传入maskmask_array二值掩膜或指定connectivity8虽默认就是 8但显式写出防误验证命令print(np.unique(mask_array))确保只有[0,1]若有255或-1先mask_array (mask_array 0).astype(np.uint8)。5.4 现象合并后的多边形出现“空心环”外圈有内圈空原因morphology.binary_fill_holes填充的是 4-邻域孔洞但遥感影像中“岛屿”常由 8-邻域定义或fill_holesTrue时binary_fill_holes对细长孔洞如河道无效。解决改用scipy.ndimage.binary_fill_holes(mask_array, structurenp.ones((3,3)))显式指定 8-邻域结构或分步先binary_dilation膨胀岛屿再binary_fill_holes最后binary_erosion验证命令plt.imshow(mask_array, cmapgray); plt.show()直观检查孔洞形态。5.5 现象rasterio.features.shapes生成的 GeoJSON 在 QGIS 中坐标系显示为Unknown CRS原因crs_wkt字符串含不可见 Unicode 字符如零宽空格或rasterio版本 1.2 无法解析复杂 WKT。解决用crs_wkt.encode(ascii, ignore).decode(ascii)清洗字符串或降级使用crsEPSG:32649若已知 EPSG 码验证命令rasterio.crs.CRS.from_wkt(crs_wkt).to_epsg()应返回有效数字否则crs_wkt无效。6. 进阶技巧用“分块处理 缓存”突破内存墙处理万兆级遥感影像当影像超过 2GB如 20000×20000 像素ReadAsArray()直接 OOM。我们不用牺牲精度换小图而是用 GDAL 的分块读取能力 joblib缓存中间结果实现“内存可控、结果一致”。6.1 分块读取与处理用 GDAL 的ReadAsArray(xoff, yoff, xsize, ysize)GDAL 支持按窗口读取避免加载全图。关键是保持块间重叠overlap防止边缘分割断裂。from joblib import Parallel, delayed import numpy as np def process_tile(args): 单块处理函数接收 (xoff, yoff, xsize, ysize, tiff_path, overlap) xoff, yoff, xsize, ysize, tiff_path, overlap args ds gdal.Open(tiff_path, gdal.GA_ReadOnly) band ds.GetRasterBand(4) # 读取带重叠的块重叠区用于边缘平滑 read_xoff max(0, xoff - overlap) read_yoff max(0, yoff - overlap) read_xsize min(xsize 2 * overlap, ds.RasterXSize - read_xoff) read_ysize min(ysize 2 * overlap, ds.RasterYSize - read_yoff) tile_arr band.ReadAsArray(read_xoff, read_yoff, read_xsize, read_ysize) ds None # 分割此处调用 2.2 节的 segment_water_by_nir tile_mask segment_water_by_nir(tile_arr) # 裁剪回无重叠区域只保留中心块结果 crop_xoff overlap if xoff 0 else 0 crop_yoff overlap if yoff 0 else 0 crop_xsize min(xsize, tile_mask.shape[1] - crop_xoff) crop_ysize min(ysize, tile_mask.shape[0] - crop_yoff) return tile_mask[crop_yoff:crop_yoffcrop_ysize, crop_xoff:crop_xoffcrop_xsize] def tiled_segmentation(tiff_path, tile_size2048, overlap64, n_jobs4): ds gdal.Open(tiff_path) width, height ds.RasterXSize, ds.RasterYSize ds None # 生成所有块坐标 tiles [] for y in range(0, height, tile_size): for x in range(0, width, tile_size): xsize min(tile_size, width - x) ysize min(tile_size, height - y) tiles.append((x, y, xsize, ysize, tiff_path, overlap)) # 并行处理 masks Parallel(n_jobsn_jobs)( delayed(process_tile)(tile_args) for tile_args in tiles ) # 拼接结果需按坐标放置 full_mask np.zeros((height, width), dtypenp.uint8) idx 0 for y in range(0, height, tile_size): for x in range(0, width, tile_size): h, w masks[idx].shape full_mask[y:yh, x:xw] masks[idx] idx 1 return full_mask # 使用示例处理万兆影像 # big_mask tiled_segmentation(huge_landsat.tif, tile_size4096, overlap128, n_jobs6)关键参数说明tile_size2048每块 2048×2048 像素在 16GB 内存机器上安全overlap64重叠 64 像素足够覆盖block_size31的局部阈值影响范围n_jobs4CPU 核数别设成os.cpu_count()留 1 核给系统。6.2 缓存中间结果用joblib.Memory避免重复计算分割耗时操作如threshold_local可缓存。joblib.Memory会根据函数输入自动生成哈希键相同输入直接返回磁盘缓存。from joblib import Memory # 创建缓存目录 mem Memory(location./cache_dir, verbose0) mem.cache def cached_threshold_local(arr, block_size, offset): return filters.threshold_local(arr, block_sizeblock_size, offsetoffset) # 在 segment_water_by_nir 中替换原调用 # local_thresh cached_threshold_local(stretched, block_size31, offset0.05)我的习惯每次启动新项目第一件事就是mkdir cache_dir chmod 700 cache_dir处理同一区域不同时相影像时缓存命中率超 80%省下 2/3 时间。遇到MemoryError我第一反应不是加内存而是看cache_dir是否写满——曾因 SSD 满导致缓存写失败程序默默重算 3 小时。希望帮到你。本文还有配套的精品资源点击获取