gdal_polygonize 示例

gdal_polygonize 是一个用于将栅格数据转换为多边形要素的GDAL工具函数，可以将栅格图像转换为矢量数据，可以用于许多不同类型的GIS分析。

示例代码

import gdal
import ogr
import osr

# 定义输入数据文件路径
input_raster = "input.tif"

# 打开输入数据文件
ds = gdal.Open(input_raster)

# 获取栅格图像的元数据
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
driver = ds.GetDriver().ShortName
projection = ds.GetProjection()
geotransform = ds.GetGeoTransform()

# 获取栅格数据
band = ds.GetRasterBand(1)
data = band.ReadAsArray()

# 定义输出矢量数据文件名
output_vector = "output.shp"

# 创建输出矢量数据的数据源
driver = ogr.GetDriverByName("ESRI Shapefile")
ds_out = driver.CreateDataSource(output_vector)

# 定义输出矢量数据的投影和坐标系
srs = osr.SpatialReference()
srs.ImportFromWkt(projection)

# 创建输出矢量数据的图层
layer_name = os.path.splitext(os.path.basename(output_vector))[0]
layer = ds_out.CreateLayer(layer_name, srs=srs)

# 添加要素的字段
field_id = ogr.FieldDefn("id", ogr.OFTInteger)
field_id.SetWidth(10)
layer.CreateField(field_id)

# 使用 gdal_polygonize 将栅格数据转换成多边形要素
gdal.Polygonize(band, None, layer, 0, [], callback=None)

# 保存矢量数据
ds_out.Destroy()

详细解释

打开输入数据文件

input_raster = "input.tif"
ds = gdal.Open(input_raster)

首先需要指定要处理的栅格图像的文件名，然后使用 gdal.Open 打开该文件，并将返回的数据集赋给变量 ds。

获取栅格图像的元数据

cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
driver = ds.GetDriver().ShortName
projection = ds.GetProjection()
geotransform = ds.GetGeoTransform()

可以使用 GDAL 提供的各种函数来获取栅格图像的元数据，例如图像的行数、列数、波段数、驱动程序名称、投影信息和地理变换信息等。

获取栅格数据

band = ds.GetRasterBand(1)
data = band.ReadAsArray()

通过 ds.GetRasterBand(1) 来获取第一个波段（如果有多个波段），然后使用 band.ReadAsArray() 将其读取到内存中，结果存储在名为 data 的 NumPy 数组中。

创建输出矢量数据的数据源

output_vector = "output.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
ds_out = driver.CreateDataSource(output_vector)

需要指定要生成的矢量数据文件的名称和类型，然后使用 ogr.GetDriverByName 获取合适的驱动程序，并将其用于创建数据源。

定义输出矢量数据的投影和坐标系

srs = osr.SpatialReference()
srs.ImportFromWkt(projection)

在创建矢量数据之前，需要明确的定义输出矢量数据的投影和坐标系信息。可以使用 osr.SpatialReference() 创建一个空间参考对象，然后使用 ImportFromWkt 函数将投影信息导入该对象。

创建输出矢量数据的图层

layer_name = os.path.splitext(os.path.basename(output_vector))[0]
layer = ds_out.CreateLayer(layer_name, srs=srs)

使用 os.path.splitext 和 os.path.basename 来分别获取输出矢量数据的名称和扩展名，然后使用 ds_out.CreateLayer 函数创建一个新的图层对象，该函数需要指定一个名称和一个空间参考对象。

添加要素的字段

field_id = ogr.FieldDefn("id", ogr.OFTInteger)
field_id.SetWidth(10)
layer.CreateField(field_id)

为了将新的要素添加到矢量数据中，可能需要在图层中定义一些字段。可以使用 ogr.FieldDefn 创建一个名为 id 的整数类型字段，并通过 layer.CreateField 函数将其添加到图层中。

转换栅格数据到多边形要素

gdal.Polygonize(band, None, layer, 0, [], callback=None)

使用 gdal.Polygonize 函数来将栅格数据转换成多边形要素，并将其添加到指定的图层中。可以通过设置不同的参数来控制结果的输出方式、精度和格式。

保存矢量数据

ds_out.Destroy()

最后，需要在程序末尾调用 ds_out.Destroy() 来关闭输出矢量数据源，并将它们保存到磁盘上。