遥感图像处理:栅格转矢量面、重新赋予栅格投影坐标系

---

前言

仅仅是分享我在做大棚识别觉得有意义的几个方法，希望别过期

---

一、栅格转矢量面

import os

from osgeo import gdal, ogr, osr
import cv2
def raster_to_shapefile(raster_path, shapefile_path, threshold_value=128):
    # Open the raster file
    raster_ds = gdal.Open(raster_path)
    band = raster_ds.GetRasterBand(1)
    array = band.ReadAsArray()

    # Apply threshold to create binary image
    binary_array = (array > threshold_value).astype(int)

    # Create a temporary raster file with the binary image
    temp_raster_path = 'temp_binary.tif'
    driver = gdal.GetDriverByName('GTiff')
    temp_raster_ds = driver.Create(temp_raster_path, raster_ds.RasterXSize, raster_ds.RasterYSize, 1, gdal.GDT_Byte)
    temp_raster_ds.SetGeoTransform(raster_ds.GetGeoTransform())
    temp_raster_ds.SetProjection(raster_ds.GetProjection())
    temp_raster_ds.GetRasterBand(1).WriteArray(binary_array)
    temp_raster_ds.FlushCache()

    # Create shapefile
    shp_driver = ogr.GetDriverByName('ESRI Shapefile')
    if shp_driver is None:
        raise Exception("Shapefile driver not available.")

    # Create new shapefile
    if os.path.exists(shapefile_path):
        shp_driver.DeleteDataSource(shapefile_path)
    shp_ds = shp_driver.CreateDataSource(shapefile_path)

    # Define the spatial reference system
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromWkt(raster_ds.GetProjection())

    # Create a new layer
    layer = shp_ds.CreateLayer('raster_layer', srs=spatial_ref, geom_type=ogr.wkbPolygon)

    # Create a field for IDs
    field_def = ogr.FieldDefn('ID', ogr.OFTInteger)
    layer.CreateField(field_def)

    # Raster to polygon conversion
    gdal.Polygonize(temp_raster_ds.GetRasterBand(1), None, layer, 0, [], callback=None)

    # Clean up
    temp_raster_ds = None
    raster_ds = None
    shp_ds = None
    os.remove(temp_raster_path)

以及我项目中存在一个最大矢量面的问题

def remove_largest_polygon_and_convert_to_shapefile(raster_path, shapefile_path, threshold_value=128):
    # 打开栅格文件
    raster_ds = gdal.Open(raster_path)
    if raster_ds is None:
        raise Exception(f"无法打开栅格文件: {raster_path}")

    # 获取栅格波段并读取数组数据
    band = raster_ds.GetRasterBand(1)
    array = band.ReadAsArray()

    # 应用阈值进行二值化处理
    binary_array = (array > threshold_value).astype(int)

    # 创建临时栅格文件
    temp_raster_path = 'temp_binary.tif'
    driver = gdal.GetDriverByName('GTiff')
    temp_raster_ds = driver.Create(temp_raster_path, raster_ds.RasterXSize, raster_ds.RasterYSize, 1, gdal.GDT_Byte)

    # 设置临时栅格的地理变换和投影
    temp_raster_ds.SetGeoTransform(raster_ds.GetGeoTransform())
    temp_raster_ds.SetProjection(raster_ds.GetProjection())
    temp_raster_ds.GetRasterBand(1).WriteArray(binary_array)
    temp_raster_ds.FlushCache()

    # 创建 shapefile
    shp_driver = ogr.GetDriverByName('ESRI Shapefile')
    if shp_driver is None:
        raise Exception("Shapefile 驱动不可用")

    # 如果文件已存在，则删除
    if os.path.exists(shapefile_path):
        shp_driver.DeleteDataSource(shapefile_path)
    shp_ds = shp_driver.CreateDataSource(shapefile_path)

    # 从栅格中提取空间参考系统
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromWkt(raster_ds.GetProjection())

    # 打印空间参考系统信息
    spatial_ref_wkt = spatial_ref.ExportToWkt()
    print("空间参考系统 WKT 信息:")
    print(spatial_ref_wkt)

    # 创建图层并应用空间参考系统
    layer = shp_ds.CreateLayer('raster_layer', srs=spatial_ref, geom_type=ogr.wkbPolygon)

    # 为图层创建一个 ID 字段
    field_def = ogr.FieldDefn('ID', ogr.OFTInteger)
    layer.CreateField(field_def)

    # 将栅格转换为多边形
    gdal.Polygonize(temp_raster_ds.GetRasterBand(1), None, layer, 0, [], callback=None)

    # 删除最大多边形
    layer_defn = layer.GetLayerDefn()
    max_area = 0
    max_feature = None
    for feature in layer:
        geom = feature.GetGeometryRef()
        area = geom.GetArea()
        if area > max_area:
            max_area = area
            max_feature = feature

    if max_feature:
        layer.DeleteFeature(max_feature.GetFID())

    # 清理临时数据
    temp_raster_ds = None
    raster_ds = None
    shp_ds = None
    os.remove(temp_raster_path)
    print(f"已创建 shapefile 并删除最大多边形: {shapefile_path}")

二、赋予矢量栅格投影坐标系

import os
from osgeo import gdal, ogr, osr
import cv2
def copy_geoCoordSys(tif_proj, tif_geotrans,Image_Line):
    def def_geoCoordSys(read_path, img_transf, img_proj):
        array_dataset = gdal.Open(read_path)
        img_array = array_dataset.ReadAsArray(0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)
        if 'int8' in img_array.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in img_array.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        if len(img_array.shape) == 3:
            img_bands, im_height, im_width = img_array.shape
        else:
            img_bands, (im_height, im_width) = 1, img_array.shape
        filename = read_path[:-4] + '_proj' + read_path[-4:]
        driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动
        dataset = driver.Create(filename, im_width, im_height, img_bands, datatype)
        dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数
        dataset.SetProjection(img_proj)  # 写入投影

        # 写入影像数据
        if img_bands == 1:
            dataset.GetRasterBand(1).WriteArray(img_array)
        else:
            for i in range(img_bands):
                dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
        print(read_path, 'geoCoordSys get!')

    cv2.imwrite('pre.tif',img=Image_Line)
    def_geoCoordSys('pre.tif', tif_geotrans, tif_proj)

使用案例

来源于大棚识别的项目的部分代码

# 线化处理：使用霍夫变换直线检测来提取直线
	dataloader, original_width, original_height, tif_proj, tif_geotrans = get_dataset(law_tif_path, stride=128, size=(256, 256))
    median_img = cv2.medianBlur(head_img-edge_img, 5)
    # 二值化
    ret, BW_Original = cv2.threshold(median_img, 128, 255, cv2.THRESH_BINARY)
    # 生成二值化的边缘分割图以及地块图
    cv2.imwrite(r"pre2.tif", BW_Original)
    cv2.imwrite(r"edge.tif", edge_img)
    # 将投影坐标系注入进生成的pre_proj.tif文件中，并生成output.shp的矢量数据
    copy_geoCoordSys(tif_proj, tif_geotrans, Image_Line=median_img)
    remove_largest_polygon_and_convert_to_shapefile("pre_proj.tif", "output.shp")

---

总结

最终，它们的分割图会被赋予原来的坐标系 ，也刚好在底图上显示，分割效果不怎么好其实。