GDAL栅格矢量化

GDAL提供了栅格矢量化等很给力的算法，但是好多算法都是通过Python脚本来提供的，对于没有安装Python环境的用户来说，这些非常有用的功能得到了很大程度的限制。GDAL工具中使用Python提供的就有栅格矢量化的功能，通过实验测试，将分类图进行矢量化后，能够很好的和原图进行匹配，而且也没有错误的多边形，下面就对GDAL中该功能做一个简单的说明。

GDAL栅格矢量化Python脚本分析，其位置在GDAL源代码目录下的/swig/python/scripts/gdal_polygonize.py，其代码如下：

try:
    from osgeo import gdal, ogr, osr
except ImportError:
    import gdal, ogr, osr

import sys
import os.path

def Usage():
    print("""
gdal_polygonize [-o name=value] [-nomask] [-mask filename] raster_file [-b band]
                [-q] [-f ogr_format] out_file [layer] [fieldname]
""")
    sys.exit(1)

# =============================================================================
#     Mainline
# =============================================================================

format = 'GML'
options = []
quiet_flag = 0
src_filename = None
src_band_n = 1

dst_filename = None
dst_layername = None
dst_fieldname = None
dst_field = -1

mask = 'default'

gdal.AllRegister()
argv = gdal.GeneralCmdLineProcessor( sys.argv )
if argv is None:
    sys.exit( 0 )

# Parse command line arguments.
i = 1
while i < len(argv):
    arg = argv[i]

    if arg == '-f':
        i = i + 1
        format = argv[i]

    elif arg == '-q' or arg == '-quiet':
        quiet_flag = 1

    elif arg == '-nomask':
        mask = 'none'

    elif arg == '-mask':
        i = i + 1
        mask = argv[i]

    elif arg == '-b':
        i = i + 1
        src_band_n = int(argv[i])

    elif src_filename is None:
        src_filename = argv[i]

    elif dst_filename is None:
        dst_filename = argv[i]

    elif dst_layername is None:
        dst_layername = argv[i]

    elif dst_fieldname is None:
        dst_fieldname = argv[i]

    else:
        Usage()

    i = i + 1

if src_filename is None or dst_filename is None:
    Usage()

if dst_layername is None:
    dst_layername = 'out'

# =============================================================================
#     Verify we have next gen bindings with the polygonize method.
# =============================================================================
try:
    gdal.Polygonize
except:
    print('')
    print('gdal.Polygonize() not available.  You are likely using "old gen"')
    print('bindings or an older version of the next gen bindings.')
    print('')
    sys.exit(1)

# =============================================================================
#    Open source file
# =============================================================================

src_ds = gdal.Open( src_filename )

if src_ds is None:
    print('Unable to open %s' % src_filename)
    sys.exit(1)

srcband = src_ds.GetRasterBand(src_band_n)

if mask is 'default':
    maskband = srcband.GetMaskBand()
elif mask is 'none':
    maskband = None
else:
    mask_ds = gdal.Open( mask )
    maskband = mask_ds.GetRasterBand(1)

# =============================================================================
#       Try opening the destination file as an existing file.
# =============================================================================

try:
    gdal.PushErrorHandler( 'QuietErrorHandler' )
    dst_ds = ogr.Open( dst_filename, update=1 )
    gdal.PopErrorHandler()
except:
    dst_ds = None

# =============================================================================
#     Create output file.
# =============================================================================
if dst_ds is None:
    drv = ogr.GetDriverByName(format)
    if not quiet_flag:
        print('Creating output %s of format %s.' % (dst_filename, format))
    dst_ds = drv.CreateDataSource( dst_filename )

# =============================================================================
#       Find or create destination layer.
# =============================================================================
try:
    dst_layer = dst_ds.GetLayerByName(dst_layername)
except:
    dst_layer = None

if dst_layer is None:

    srs = None
    if src_ds.GetProjectionRef() != '':
        srs = osr.SpatialReference()
        srs.ImportFromWkt( src_ds.GetProjectionRef() )

    dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs )

    if dst_fieldname is None:
        dst_fieldname = 'DN'

    fd = ogr.FieldDefn( dst_fieldname, ogr.OFTInteger )
    dst_layer.CreateField( fd )
    dst_field = 0
else:
    if dst_fieldname is not None:
        dst_field = dst_layer.GetLayerDefn().GetFieldIndex(dst_fieldname)
        if dst_field < 0:
            print("Warning: cannot find field '%s' in layer '%s'" % (dst_fieldname, dst_layername))

# =============================================================================
#    Invoke algorithm.
# =============================================================================

if quiet_flag:
    prog_func = None
else:
    prog_func = gdal.TermProgress

result = gdal.Polygonize( srcband, maskband, dst_layer, dst_field, options,
                          callback = prog_func )

srcband = None
src_ds = None
dst_ds = None
mask_ds = None

同时该工具的说明文档见http://www.gdal.org/gdal_polygonize.html，英文的，不过都很容易看明白。查看Python代码发现，其中调用的最重要的函数是一个叫 gdal.Polygonize的函数，在GDAL算法说明文档中找到了这个函数说明，地址为：http://www.gdal.org/gdal__alg_8h.html#3f522a9035d3512b5d414fb4752671b1，如下：

CPLErr CPL_DLL CPL_STDCALL
GDALPolygonize( GDALRasterBandH hSrcBand,                 //输入栅格图像波段
                GDALRasterBandH hMaskBand,                //掩码图像波段，可以为NULL
                OGRLayerH hOutLayer,                      //矢量化后的矢量图层
                int iPixValField,                         //需要将像元DN值写入矢量属性字段的字段索引
                char **papszOptions,                      //算法选项，目前算法中没有用到，设置为NULL即可
                GDALProgressFunc pfnProgress,             //进度条回调函数
                void * pProgressArg );                    //进度条参数

通过对上面的函数参数进行分析，就可以自己直接调用该函数，使用C/C++语言来对其进行封装，达到我们的目的，下面是我对该函数的一个简单封装，指定一个输入的分类图像，指定一个输出的shp文件，就可以将该分类图像进行矢量化。

/**
* @brief 分类后处理之栅格矢量化
* @param pszSrcFile       输入的分类文件，如果输入文件是多波段文件，则只处理第一波段
* @param pszDstFile       输出结果文件，一个矢量文件
* @param pszFormat        输出文件格式，默认为ERSI Shapefile
* @param pProgress        进度条指针
* @return 成功返回0
*/
int ImagePolygonize(const char* pszSrcFile, const char* pszDstFile, const char* pszFormat, LT_Progress *pProgress)
{
    if (pProgress != NULL)
    {
        pProgress->ReSetting();
        pProgress->SetProgressCaption("提示");
        pProgress->SetProgressTip("开始计算栅格矢量化...");
    }

    GDALAllRegister();
    OGRRegisterAll();

    GDALDataset* poSrcDS = (GDALDataset*) GDALOpen(pszSrcFile, GA_ReadOnly);    //打开栅格图像
    if (poSrcDS == NULL)
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("不能打开指定文件，请检查文件是否存在！");

        return RE_NOFILE;
    }

    OGRSFDriver* poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszFormat);
    if (poDriver == NULL)
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("不能创建指定类型文件！");

        GDALClose((GDALDatasetH)poSrcDS);

        return RE_CREATEFILE;
    }

    OGRDataSource* poDstDS = poDriver->CreateDataSource(pszDstFile, NULL); //创建输出矢量文件
    if (poDstDS == NULL)
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("不能创建指定类型文件！");

        GDALClose((GDALDatasetH)poSrcDS);

        return RE_CREATEFILE;
    }

    OGRSpatialReference *poSpatialRef = new OGRSpatialReference(poSrcDS->GetProjectionRef());
    string strLayerName = GetLayerName(pszDstFile);

    OGRLayer* poLayer = poDstDS->CreateLayer(strLayerName.c_str(), poSpatialRef, wkbPolygon, NULL);
    if (poLayer == NULL)
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("创建矢量图层失败！");

        GDALClose((GDALDatasetH)poSrcDS);
        OGRDataSource::DestroyDataSource(poDstDS);

        delete poSpatialRef;
        poSpatialRef = NULL;

        return RE_CREATEFILE;
    }

    OGRFieldDefn ofieldDef("DN", OFTInteger);    //创建属性表，只有一个字段即“DN”，里面保存对应的栅格的像元值
    if (poLayer->CreateField(&ofieldDef) != OGRERR_NONE)
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("创建矢量图层属性表失败！");

        GDALClose((GDALDatasetH)poSrcDS);
        OGRDataSource::DestroyDataSource(poDstDS);

        delete poSpatialRef;
        poSpatialRef = NULL;

        return RE_CREATEFILE;
    }

    GDALRasterBandH hSrcBand = (GDALRasterBandH) poSrcDS->GetRasterBand(1);    //获取图像的第一个波段

    if (GDALPolygonize(hSrcBand, NULL, (OGRLayerH)poLayer, 0, NULL, GDALProgress, pProgress) != CE_None)//调用栅格矢量化
    {
        if (pProgress != NULL)
            pProgress->SetProgressTip("计算失败！");

        GDALClose((GDALDatasetH)poSrcDS);
        OGRDataSource::DestroyDataSource(poDstDS);

        delete poSpatialRef;
        poSpatialRef = NULL;

        return RE_FAILD;
    }

    GDALClose((GDALDatasetH)poSrcDS);    //关闭文件
    OGRDataSource::DestroyDataSource(poDstDS);

    delete poSpatialRef;
    poSpatialRef = NULL;

    if (pProgress != NULL)
        pProgress->SetProgressTip("计算成功！");

    return RE_SUCCESS;
}

在调用的时候，只需要指定两个文件路径即可，如下：

void main()
{
    LT_ConsoleProgress *pProgress = new LT_ConsoleProgress(); //进度条
    progress_timer *pTime = new progress_timer;    //计时

    int f = ImagePolygonize("C://Work//Data//Classify.tif","C://Work//Data//Classify.shp", "ESRI Shapefile", pProgress);

    if (f == RE_SUCCESS)
        printf("计算成功/n");
    else
        printf("计算失败/n");

    delete pProgress;    //释放资源
    delete pTime;
    system("pause");
}

测试图像如下，

矢量化后的矢量（使用要素值渲染方式）

放大细节处对比（下左图为栅格图像，下右图为矢量化后的矢量）：