【Python&GIS】基于Python实现栅格转面、面转栅格（栅格、矢量互转）

​

        各位好，我又来水文章了。最近因为同事在做生态服务相关的项目，需要对矢量数据进行操作，然后我就查了查相关资料，今天就和大家分享一下如何使用Python的GDAL库实现栅格转要素、要素转栅格（栅格、矢量互相转换）。其实我之前已经分享过栅格转面和计算要素面积的代码，大家感兴趣可以去看下：【Python&GIS】GDAL栅格转面&计算矢量面积

一、栅格转面

        这里栅格对应的值被写入矢量的属性字段Value中。注释很详细不多说。

def Raster_To_Vector(input_raster, output_shp):
    """
    :param input_raster: 输入需要转换的栅格数据
    :param output_shp: 输出矢量的路径
    :return：None
    """
    print("正在栅格转面。。。")
    # 栅格转面
    ds_raster = gdal.Open(input_raster)  # 读取路径中的栅格数据
    band_raster = ds_raster.GetRasterBand(1)  # 获取需要转为矢量的波段
    proj_raster = osr.SpatialReference()
    proj_raster.ImportFromWkt(ds_raster.GetProjection())
    # 将栅格数据的投影信息赋值给矢量
    driver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(output_shp):  # 若文件已经存在，则删除它继续重新做一遍
        driver.DeleteDataSource(output_shp)
    polygon = driver.CreateDataSource(output_shp)  # 创建数据资源
    layer_polygon = polygon.CreateLayer("Shp", srs=proj_raster, geom_type=ogr.wkbMultiPolygon)  # 创建图层，定义多面
    new_field = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
    layer_polygon.CreateField(new_field)
    gdal.FPolygonize(band_raster, None, layer_polygon, 0)  # 核心函数，执行的就是栅格转矢量操作
    polygon.SyncToDisk()
    polygon = None

二、面转栅格

        这里有个大问题！！！GDAL里面的转换函数里的参数是需要参考值的，所以我们需要输入一个栅格作为参考，才能将矢量转为栅格。这就导致这个函数的应用非常局限，而GIS中的面转栅格就不要参考栅格，甚至还可以选择栅格大小。

        无语......，但这个函数也不是一无是处，这个函数可以将我们的深度学习样本转为栅格，懂我意思？我们在画样本时都是矢量，但有些模型的输入是需要栅格数据的，这个时候这个函数还是有点用的。

def Vector_To_Raster(input_shp, refer_raster, attribute, output_raster):
    """
    :param input_shp: 输入需要转换的矢量数据
    :param refer_raster: 输入参考栅格数据
    :param attribute: 输入栅格值对应的矢量字段
    :param output_raster: 输出栅格数据路径
    :return: None
    """
    ds_raster = gdal.Open(refer_raster)
    ds_proj = ds_raster.GetProjection()  # 投影信息
    ds_trans = ds_raster.GetGeoTransform()  # 仿射地理变换参数
    ds_width = ds_raster.RasterXSize  # 获取宽度/列数
    ds_height = ds_raster.RasterYSize  # 获取高度/行数
    ds_raster = None  # 释放内存
    del ds_raster
    driver = ogr.GetDriverByName("ESRI Shapefile")
    ds = driver.Open(input_shp, 1)
    layer = ds.GetLayer()  # 获取图层文件对象
    result = gdal.GetDriverByName('GTiff').Create(output_raster, ds_width, ds_height, bands=1, eType=gdal.GDT_Byte)
    result.SetGeoTransform(ds_trans)  # 写入仿射地理变换参数
    result.SetProjection(ds_proj)  # 写入投影信息
    band = result.GetRasterBand(1)
    band.SetNoDataValue(0)  # 忽略背景值
    band.FlushCache()  # 清空数据缓存
    # options = ["ATTRIBUTE=attribute", "CHUNKY SIZE=0", "ALL_TOUCHED=False"]
    # 指定输入矢量数据的属性字段中的字段值作为栅格值写入栅格文件中，该值将输出到所有输出波段中。假如该值指定了，burn_Values参数的值将失效数可以设置为空。
    # 指定该运行操作的块的高度。该值越大所需的计算时间越小。如果该值没有设置或者设置为0则由GDAL的缓存大小根据公式：缓存所占的字节数/扫描函数的字节数得到。所以该值不会超出缓存的大小。
    # 设置为TRUE表示所有的像素接触到矢量中线或多边形，否则只是多边形中心或被Bresenham算法选中的部分。默认值为FALSE。简单来说，FALSE是八方向栅格化，TRUE是全路径栅格化。
    gdal.RasterizeLayer(result, [1], layer, burn_values=[1], options=[f"ATTRIBUTE={attribute}"])
    # 矢量转栅格函数。输出栅格，波段数，图层，栅格数值（可控失效），可选参数（栅格数值取字段）
    result = None
    del result, layer

三、完整代码

        没啥说的，注释很详细，用哪个函数就调用哪个函数即可！

# -*- coding: utf-8 -*-
"""
@Time ： 2023/11/3 14:37
@Auth ： RS迷途小书童
@File ：Convert Raster And Vector.py
@IDE ：PyCharm
@Purpose：栅格数据、矢量数据互相转换
"""
import os
from osgeo import gdal, osr, ogr


def Raster_To_Vector(input_raster, output_shp):
    """
    :param input_raster: 输入需要转换的栅格数据
    :param output_shp: 输出矢量的路径
    :return：None
    """
    print("正在栅格转面。。。")
    # 栅格转面
    ds_raster = gdal.Open(input_raster)  # 读取路径中的栅格数据
    band_raster = ds_raster.GetRasterBand(1)  # 获取需要转为矢量的波段
    proj_raster = osr.SpatialReference()
    proj_raster.ImportFromWkt(ds_raster.GetProjection())
    # 将栅格数据的投影信息赋值给矢量
    driver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(output_shp):  # 若文件已经存在，则删除它继续重新做一遍
        driver.DeleteDataSource(output_shp)
    polygon = driver.CreateDataSource(output_shp)  # 创建数据资源
    layer_polygon = polygon.CreateLayer("Shp", srs=proj_raster, geom_type=ogr.wkbMultiPolygon)  # 创建图层，定义多面
    new_field = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
    layer_polygon.CreateField(new_field)
    gdal.FPolygonize(band_raster, None, layer_polygon, 0)  # 核心函数，执行的就是栅格转矢量操作
    polygon.SyncToDisk()
    polygon = None


def Vector_To_Raster(input_shp, refer_raster, attribute, output_raster):
    """
    :param input_shp: 输入需要转换的矢量数据
    :param refer_raster: 输入参考栅格数据
    :param attribute: 输入栅格值对应的矢量字段
    :param output_raster: 输出栅格数据路径
    :return: None
    """
    ds_raster = gdal.Open(refer_raster)
    ds_proj = ds_raster.GetProjection()  # 投影信息
    ds_trans = ds_raster.GetGeoTransform()  # 仿射地理变换参数
    ds_width = ds_raster.RasterXSize  # 获取宽度/列数
    ds_height = ds_raster.RasterYSize  # 获取高度/行数
    ds_raster = None  # 释放内存
    del ds_raster
    driver = ogr.GetDriverByName("ESRI Shapefile")
    ds = driver.Open(input_shp, 1)
    layer = ds.GetLayer()  # 获取图层文件对象
    result = gdal.GetDriverByName('GTiff').Create(output_raster, ds_width, ds_height, bands=1, eType=gdal.GDT_Byte)
    result.SetGeoTransform(ds_trans)  # 写入仿射地理变换参数
    result.SetProjection(ds_proj)  # 写入投影信息
    band = result.GetRasterBand(1)
    band.SetNoDataValue(0)  # 忽略背景值
    band.FlushCache()  # 清空数据缓存
    # options = ["ATTRIBUTE=attribute", "CHUNKY SIZE=0", "ALL_TOUCHED=False"]
    # 指定输入矢量数据的属性字段中的字段值作为栅格值写入栅格文件中，该值将输出到所有输出波段中。假如该值指定了，burn_Values参数的值将失效数可以设置为空。
    # 指定该运行操作的块的高度。该值越大所需的计算时间越小。如果该值没有设置或者设置为0则由GDAL的缓存大小根据公式：缓存所占的字节数/扫描函数的字节数得到。所以该值不会超出缓存的大小。
    # 设置为TRUE表示所有的像素接触到矢量中线或多边形，否则只是多边形中心或被Bresenham算法选中的部分。默认值为FALSE。简单来说，FALSE是八方向栅格化，TRUE是全路径栅格化。
    gdal.RasterizeLayer(result, [1], layer, burn_values=[1], options=[f"ATTRIBUTE={attribute}"])
    # 矢量转栅格函数。输出栅格，波段数，图层，栅格数值（可控失效），可选参数（栅格数值取字段）
    result = None
    del result, layer


if __name__ == "__main__":
    shp = r"B:\Personal\彭俊喜\深度学习\CNN\0919小块实验/sample.shp"
    refer = r"B:\Personal\彭俊喜\深度学习\CNN\0919小块实验\example.tif"
    out = r"B:\Personal\彭俊喜\深度学习\CNN\0919小块实验\try.tif"
    Attribute = "AREA"
    Vector_To_Raster(shp, refer, Attribute, out)

        本文章主要是分享个人在学习Python过程中写过的一些代码。有些部分借鉴了前人以及官网的教程，如有侵权请联系作者删除，大家有问题可以随时留言交流，博主会及时回复。