GDAL栅格矢量化

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: }

测试图像如下，

image

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

image

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

imageimage

posted on 2011-02-23 21:12  王大王  阅读(593)  评论(0)    收藏  举报

导航