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

GDAL栅格矢量化Python脚本分析，其位置在GDAL源代码目录下的

/swig/python/scripts/gdal_polygonize.py，其代码如下：1:2: try:3: from osgeo import gdal, ogr, osr4: except ImportError:5: import gdal, ogr, osr6:7: import sys8: import os.path9:10: def Usage():11: print("""12: gdal_polygonize [-o name=value] [-nomask] [-mask filename] raster_file [-b band]13: [-q] [-f ogr_format] out_file [layer] [fieldname]14: """)15: sys.exit(1)16:17: # =============================================================================18: # Mainline19: # =============================================================================20:21: format = 'GML'22: options = []23: quiet_flag = 024: src_filename = None25: src_band_n = 126:27: dst_filename = None28: dst_layername = None29: dst_fieldname = None30: dst_field = -131:32: mask = 'default'33:34: gdal.AllRegister()35: argv = gdal.GeneralCmdLineProcessor( sys.argv )36: if argv is None:37: sys.exit( 0 )38:39: # Parse command line arguments.40: i = 141: while i < len(argv):42: arg = argv[i]43:44: if arg == '-f':45: i = i + 146: format = argv[i]47:48: elif arg == '-q' or arg == '-quiet':49: quiet_flag = 150:51: elif arg == '-nomask':52: mask = 'none'53:54: elif arg == '-mask':55: i = i + 156: mask = argv[i]57:58: elif arg == '-b':59: i = i + 160: src_band_n = int(argv[i])61:62: elif src_filename is None:63: src_filename = argv[i]64:65: elif dst_filename is None:66: dst_filename = argv[i]67:68: elif dst_layername is None:69: dst_layername = argv[i]70:71: elif dst_fieldname is None:72: dst_fieldname = argv[i]73:74: else:75: Usage()76:77: i = i + 178:79: if src_filename is None or dst_filename is None:80: Usage()81:82: if dst_layername is None:83: dst_layername = 'out'84:85: # =============================================================================86: # Verify we have next gen bindings with the polygonize method.87: # =============================================================================88: try:89: gdal.Polygonize90: except:91: print('')92: print('gdal.Polygonize() not available. You are likely using "old gen"')93: print('bindings or an older version of the next gen bindings.')94: print('')95: sys.exit(1)96:97: # =============================================================================98: # Open source file99: # =============================================================================100:101: src_ds = gdal.Open( src_filename )102:103: if src_ds is None:104: print('Unable to open %s' % src_filename)105: sys.exit(1)106:107: srcband = src_ds.GetRasterBand(src_band_n)108:109: if mask is 'default':110: maskband = srcband.GetMaskBand()111: elif mask is 'none':112: maskband = None113: else:114: mask_ds = gdal.Open( mask )115: maskband = mask_ds.GetRasterBand(1)116:117: # =============================================================================118: # Try opening the destination file as an existing file.119: # =============================================================================120:121: try:122: gdal.PushErrorHandler( 'QuietErrorHandler' )123: dst_ds = ogr.Open( dst_filename, update=1 )124: gdal.PopErrorHandler()125: except:126: dst_ds = None127:128: # =============================================================================129: # Create output file.130: # =============================================================================131: if dst_ds is None:132: drv = ogr.GetDriverByName(format)133: if not quiet_flag:134: print('Creating output %s of format %s.' % (dst_filename, format))135: dst_ds = drv.CreateDataSource( dst_filename )136:137: # =============================================================================138: # Find or create destination layer.139: # =============================================================================140: try:141: dst_layer = dst_ds.GetLayerByName(dst_layername)142: except:143: dst_layer = None144:145: if dst_layer is None:146:147: srs = None148: if src_ds.GetProjectionRef() != '':149: srs = osr.SpatialReference()150: srs.ImportFromWkt( src_ds.GetProjectionRef() )151:152: dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs )153:154: if dst_fieldname is None:155: dst_fieldname = 'DN'156:157: fd = ogr.FieldDefn( dst_fieldname, ogr.OFTInteger )158: dst_layer.CreateField( fd )159: dst_field = 0160: else:161: if dst_fieldname is not None:162: dst_field = dst_layer.GetLayerDefn().GetFieldIndex(dst_fieldname)163: if dst_field < 0:164: print("Warning: cannot find field '%s' in layer '%s'" % (dst_fieldname, dst_layername))165:166: # =============================================================================167: # Invoke algorithm.168: # =============================================================================169:170: if quiet_flag:171: prog_func = None172: else:173: prog_func = gdal.TermProgress174:175: result = gdal.Polygonize( srcband, maskband, dst_layer, dst_field, options,176: callback = prog_func )177:178: srcband = None179: src_ds = None180: dst_ds = None181: mask_ds = None
1: CPLErr CPL_DLL CPL_STDCALL2: GDALPolygonize( GDALRasterBandH hSrcBand, //输入栅格图像波段3: GDALRasterBandH hMaskBand, //掩码图像波段，可以为NULL4: OGRLayerH hOutLayer, //矢量化后的矢量图层5: int iPixValField, //需要将像元DN值写入矢量属性字段的字段索引6: char **papszOptions, //算法选项，目前算法中没有用到，设置为NULL即可7: GDALProgressFunc pfnProgress, //进度条回调函数8: void * pProgressArg ); //进度条参数
通过对上面的函数参数进行分析，就可以自己直接调用该函数，使用C/C++语言来对其进行封装，达到我们的目的，下面是我对该函数的一个简单封装，指定一个输入的分类图像，指定一个输出的shp文件，就可以将该分类图像进行矢量化。1: /**2: * @brief 分类后处理之栅格矢量化3: * @param pszSrcFile 输入的分类文件，如果输入文件是多波段文件，则只处理第一波段4: * @param pszDstFile 输出结果文件，一个矢量文件5: * @param pszFormat 输出文件格式，默认为ERSI Shapefile6: * @param pProgress 进度条指针7: * @return 成功返回08: */9: int ImagePolygonize(const char* pszSrcFile, const char* pszDstFile, const char* pszFormat, LT_Progress *pProgress)10: {11: if (pProgress != NULL)12: {13: pProgress->ReSetting();14: pProgress->SetProgressCaption("提示");15: pProgress->SetProgressTip("开始计算栅格矢量化...");16: }17:18: GDALAllRegister();19: OGRRegisterAll();20:21: GDALDataset* poSrcDS = (GDALDataset*) GDALOpen(pszSrcFile, GA_ReadOnly); //打开栅格图像22: if (poSrcDS == NULL)23: {24: if (pProgress != NULL)25: pProgress->SetProgressTip("不能打开指定文件，请检查文件是否存在！");26:27: return RE_NOFILE;28: }29:30: OGRSFDriver* poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszFormat);31: if (poDriver == NULL)32: {33: if (pProgress != NULL)34: pProgress->SetProgressTip("不能创建指定类型文件！");35:36: GDALClose((GDALDatasetH)poSrcDS);37:38: return RE_CREATEFILE;39: }40:41: OGRDataSource* poDstDS = poDriver->CreateDataSource(pszDstFile, NULL); //创建输出矢量文件42: if (poDstDS == NULL)43: {44: if (pProgress != NULL)45: pProgress->SetProgressTip("不能创建指定类型文件！");46:47: GDALClose((GDALDatasetH)poSrcDS);48:49: return RE_CREATEFILE;50: }51:52: OGRSpatialReference *poSpatialRef = new OGRSpatialReference(poSrcDS->GetProjectionRef());53: string strLayerName = GetLayerName(pszDstFile);54:55: OGRLayer* poLayer = poDstDS->CreateLayer(strLayerName.c_str(), poSpatialRef, wkbPolygon, NULL);56: if (poDstDS == NULL)57: {58: if (pProgress != NULL)59: pProgress->SetProgressTip("创建矢量图层失败！");60:61: GDALClose((GDALDatasetH)poSrcDS);62: OGRDataSource::DestroyDataSource(poDstDS);63:64: delete poSpatialRef;65: poSpatialRef = NULL;66:67: return RE_CREATEFILE;68: }69:70: OGRFieldDefn ofieldDef("DN", OFTInteger); //创建属性表，只有一个字段即“DN”，里面保存对应的栅格的像元值71: if (poLayer->CreateField(&ofieldDef) != OGRERR_NONE)72: {73: if (pProgress != NULL)74: pProgress->SetProgressTip("创建矢量图层属性表失败！");75:76: GDALClose((GDALDatasetH)poSrcDS);77: OGRDataSource::DestroyDataSource(poDstDS);78:79: delete poSpatialRef;80: poSpatialRef = NULL;81:82: return RE_CREATEFILE;83: }84:85: GDALRasterBandH hSrcBand = (GDALRasterBandH) poSrcDS->GetRasterBand(1); //获取图像的第一个波段86:87: if (GDALPolygonize(hSrcBand, NULL, (OGRLayerH)poLayer, 0, NULL, GDALProgress, pProgress) != CE_None)//调用栅格矢量化88: {89: if (pProgress != NULL)90: pProgress->SetProgressTip("计算失败！");91:92: GDALClose((GDALDatasetH)poSrcDS);93: OGRDataSource::DestroyDataSource(poDstDS);94:95: delete poSpatialRef;96: poSpatialRef = NULL;97:98: return RE_FAILD;99: }100:101: GDALClose((GDALDatasetH)poSrcDS); //关闭文件102: OGRDataSource::DestroyDataSource(poDstDS);103:104: delete poSpatialRef;105: poSpatialRef = NULL;106:107: if (pProgress != NULL)108: pProgress->SetProgressTip("计算成功！");109:110: return RE_SUCCESS;111: }
在调用的时候，只需要指定两个文件路径即可，如下：1: void main()2: {3: LT_ConsoleProgress *pProgress = new LT_ConsoleProgress(); //进度条4: progress_timer *pTime = new progress_timer; //计时5:6: int f = ImagePolygonize("C://Work//Data//Classify.tif","C://Work//Data//Classify.shp", "ESRI Shapefile", pProgress);7:8: if (f == RE_SUCCESS)9: printf("计算成功/n");10: else11: printf("计算失败/n");12:13: delete pProgress; //释放资源14: delete pTime;15: system("pause");16: }


测试图像如下，


矢量化后的矢量(使用要素值渲染方式)


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