Python-GDAL教程：二值图轮廓转为面矢量数据，并计算每个面矢量的面积

一 介绍

如下图，如何把该二值图的轮廓提取出来，转化为矢量文件，并计算每个面矢量的面积？

二 算法实现

1 轮廓转矢量

代码和注释如下：

# _*_ coding: utf-8 _*_

import sys
import os
import io
import cv2
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf8')  # 打印出中文字符
try:
    from osgeo import gdal
    from osgeo import ogr
    from osgeo import osr
except ImportError:
    import gdal
    import ogr
    import osr

img_path = 'D:\\CodePython\\rasters2vector\\test4.png'
img0 = cv2.imread(img_path)
img1 = cv2.cvtColor(img0, cv2.COLOR_BGR2GRAY)
ret, img2 = cv2.threshold(img1, 200, 255, cv2.THRESH_BINARY)
contours, hierarchy = cv2.findContours(
    img2, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)  # 提取二值图的所有轮廓边界
# print(contours)
# img = cv2.drawContours(img0, contours, -1, (0, 255, 0), 3)
# cv2.imshow("rotation", img)
# cv2.waitKey()
# cv2.destroyAllWindows()

gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")  # 为了支持中文路径
gdal.SetConfigOption("SHAPE_ENCODING", "CP936")  # 为了使属性表字段支持中文
strVectorFile = "D:\\CodePython\\rasters2vector\\test\\test4.shp"  # 定义写入路径及文件名
ogr.RegisterAll()  # 注册所有的驱动
strDriverName = "ESRI Shapefile"  # 创建数据，这里创建ESRI的shp文件
oDriver = ogr.GetDriverByName(strDriverName)
if oDriver == None:
    print("%s 驱动不可用！\n", strDriverName)

oDS = oDriver.CreateDataSource(strVectorFile)  # 创建数据源
if oDS == None:
    print("创建文件【%s】失败！", strVectorFile)

srs = osr.SpatialReference()  # 创建空间参考
srs.ImportFromEPSG(4326)  # 定义地理坐标系WGS1984
papszLCO = []
# 创建图层，创建一个多边形图层,"TestPolygon"->属性表名
oLayer = oDS.CreateLayer("TestPolygon", srs, ogr.wkbPolygon, papszLCO)
if oLayer == None:
    print("图层创建失败！\n")

'''下面添加矢量数据，属性表数据、矢量数据坐标'''
oFieldID = ogr.FieldDefn("FieldID", ogr.OFTInteger)  # 创建一个叫FieldID的整型属性
oLayer.CreateField(oFieldID, 1)

oDefn = oLayer.GetLayerDefn()  # 定义要素
gardens = ogr.Geometry(ogr.wkbMultiPolygon)  # 定义总的多边形集
i = 0
for contour in contours:
    area = cv2.contourArea(contour)
    if area > 100:  # 面积大于n才保存
        # print(area)
        box1 = ogr.Geometry(ogr.wkbLinearRing)
        i += 1
        for point in contour:
            x_col = float(point[0, 1])
            y_row = float(point[0, 0])
            box1.AddPoint(y_row, x_col)
        oFeatureTriangle = ogr.Feature(oDefn)
        oFeatureTriangle.SetField(0, i)
        garden1 = ogr.Geometry(ogr.wkbPolygon)  # 每次重新定义单多变形
        garden1.AddGeometry(box1)  # 将轮廓坐标放在单多边形中
        gardens.AddGeometry(garden1)  # 依次将单多边形放入总的多边形集中
gardens.CloseRings()  # 封闭多边形集中的每个单多边形，后面格式需要

geomTriangle = ogr.CreateGeometryFromWkt(str(gardens))  # 将封闭后的多边形集添加到属性表
oFeatureTriangle.SetGeometry(geomTriangle)
oLayer.CreateFeature(oFeatureTriangle)
oDS.Destroy()
print("数据矢量创建完成！\n")

上面代码只是将所有矢量数据添加到一个要素中，生成矢量图片如下图：

2 整个要素转为多个面矢量

下面代码将会把整个要素变为我们理解的一个个面矢量

# '''矢量炸开'''
# ogr2ogr参数说明 ogr2ogr -f "ESRI Shapefile" -nlt LINESTRING -explodecollections  output.shp input.shp
os.system("ogr2ogr -f \"ESRI Shapefile\" -explodecollections  \"D:\\CodePython\\rasters2vector\\test\\test5.shp\" \"D:\\CodePython\\rasters2vector\\test\\test4.shp\"")

3 计算每个面矢量的面积

下面将计算每个面矢量的面积，并将面积值添加到属性表中，代码和注释如下：

'''计算面积'''
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(
    "D:\\CodePython\\rasters2vector\\test\\test5.shp", 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
new_field.SetPrecision(2)  # 设置面积精度
layer.CreateField(new_field)
for feature in layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea() #计算面积
    # m_area = (area/(0.0089**2))*1e+6  # 单位由十进制度转为米
    # print(m_area)
    feature.SetField("Area", area) #将面积添加到属性表中
    layer.SetFeature(feature)
dataSource = None

在arcgis中打开.shp文件效果如下图：