GDAL提供了栅格矢量化等很给力的算法，但是好多算法都是通过Python脚本来提供的，对于没有安装Python环境的用户来说，这些非常有用的功能得到了很大程度的限制。GDAL工具中使用Python提供的就有栅格矢量化的功能，通过实验测试，将分类图进行矢量化后，能够很好的和原图进行匹配，而且也没有错误的多边形，下面就对GDAL中该功能做一个简单的说明。
GDAL栅格矢量化Python脚本分析，其位置在GDAL源代码目录下的/swig/python/scripts/gdal_polygonize.py，其代码如下：
1: 2: try: 3: from osgeo import gdal, ogr, osr 4: except ImportError: 5: import gdal, ogr, osr 6: 7: import sys 8: import os.path 9: 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: # Mainline
19: # =============================================================================
20: 21: format = 'GML'
22: options = [] 23: quiet_flag = 0 24: src_filename = None 25: src_band_n = 1 26: 27: dst_filename = None 28: dst_layername = None 29: dst_fieldname = None 30: dst_field = -1 31: 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 = 1 41: while i < len(argv):
42: arg = argv[i] 43: 44: if arg == '-f':
45: i = i + 1 46: format = argv[i] 47: 48: elif arg == '-q' or arg == '-quiet':
49: quiet_flag = 1 50: 51: elif arg == '-nomask':
52: mask = 'none'
53: 54: elif arg == '-mask':
55: i = i + 1 56: mask = argv[i] 57: 58: elif arg == '-b':
59: i = i + 1 60: 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 + 1 78: 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.Polygonize 90: 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 file
99: # =============================================================================
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 = None 113: 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 = None 127: 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 = None 144: 145: if dst_layer is None:
146: 147: srs = None 148: 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 = 0 160: 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 = None 172: else:
173: prog_func = gdal.TermProgress 174: 175: result = gdal.Polygonize( srcband, maskband, dst_layer, dst_field, options, 176: callback = prog_func ) 177: 178: srcband = None 179: src_ds = None 180: dst_ds = None 181: mask_ds = None同时该工具的说明文档见http://www.gdal.org/gdal_polygonize.html，英文的，不过都很容易看明白。查看Python代码发现，其中调用的最重要的函数是一个叫 gdal.Polygonize的函数，在GDAL算法说明文档中找到了这个函数说明，地址为：http://www.gdal.org/gdal__alg_8h.html#3f522a9035d3512b5d414fb4752671b1，如下：
1: CPLErr CPL_DLL CPL_STDCALL 2: GDALPolygonize( GDALRasterBandH hSrcBand, //输入栅格图像波段
3: GDALRasterBandH hMaskBand, //掩码图像波段，可以为NULL
4: 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 Shapefile
6: * @param pProgress 进度条指针
7: * @return 成功返回0
8: */
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 (poLayer == 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: else
11: printf("计算失败/n");
12: 13: delete pProgress; //释放资源
14: delete pTime;
15: system("pause");
16: }测试图像如下，
矢量化后的矢量（使用要素值渲染方式）
放大细节处对比（下左图为栅格图像，下右图为矢量化后的矢量）：
