距离上一次博客更新，起码又是大半年，时光飞逝，我也已经老了。。。这一次，我解决了一个工程上的小问题，可能在行家看来简单，但是呢，它好像又没那么简单，就是我们通常用的栅格转矢量，
我们知道栅格转矢量，通常有以下方法：采用Arcgis进行栅格转矢量，然后工程化呢，就用arcpy实现，就可以了，或者用qgis，原理也差不多，编程的话，绕不过去的，当然是GDAL，这个大家都是直接
调用GDAL的栅格转矢量API就完事了，对于小数据，或者几百MB这种数据效率也不错，问题不大，但是当你的数据规模达到几个G，或者几百个G之后，在涉密的单机上，你怎么解决，只能重写，我这里
采用GDAL从源码层实现了超大影像的栅格转矢量，相比于直接调用GDAL的函数接口实现，效率要高太多了，我们看一下栅格转矢量的实现代码：
1 def RasterToPoly(rasterName, shpName): 2 """ 3 栅格转矢量 4 :param rasterName: 输入分类后的栅格名称 5 :param shpName: 输出矢量名称 6 :return: 7 """ 8 inraster = gdal.Open(rasterName) # 读取路径中的栅格数据 9 inband = inraster.GetRasterBand(1) # 这个波段就是最后想要转为矢量的波段，如果是单波段数据的话那就都是1 10 prj = osr.SpatialReference() 11 prj.ImportFromWkt(inraster.GetProjection()) # 读取栅格数据的投影信息，用来为后面生成的矢量做准备 12 13 outshp = shpName 14 drv = ogr.GetDriverByName("ESRI Shapefile") 15 if os.path.exists(outshp): # 若文件已经存在，则删除它继续重新做一遍 16 drv.DeleteDataSource(outshp) 17 Polygon = drv.CreateDataSource(outshp) # 创建一个目标文件 18 Poly_layer = Polygon.CreateLayer(shpName[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon) # 对shp文件创建一个图层，定义为多个面类 19 newField = ogr.FieldDefn('Value', ogr.OFTReal) # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value 20 Poly_layer.CreateField(newField) 21 22 gdal.FPolygonize(inband, None, Poly_layer, 0) # 核心函数，执行的就是栅格转矢量操作 23 Polygon.SyncToDisk() 24 Polygon = None 25 26 deleteBackground(shpName, 0) # 删除背景
然后转成矢量后，有一些Nodata图层需要删掉，调用如下代码解决：
def deleteBackground(shpName, backGroundValue): """ 删除背景,一般背景的像素值为0 """ driver = ogr.GetDriverByName('ESRI Shapefile') pFeatureDataset = driver.Open(shpName, 1) pFeaturelayer = pFeatureDataset.GetLayer(0) strValue = backGroundValue strFilter = "Value = '" + str(strValue) + "'" pFeaturelayer.SetAttributeFilter(strFilter) pFeatureDef = pFeaturelayer.GetLayerDefn() pLayerName = pFeaturelayer.GetName() pFieldName = "Value" pFieldIndex = pFeatureDef.GetFieldIndex(pFieldName) for pFeature in pFeaturelayer: pFeatureFID = pFeature.GetFID() pFeaturelayer.DeleteFeature(int(pFeatureFID)) strSQL = "REPACK " + str(pFeaturelayer.GetName()) pFeatureDataset.ExecuteSQL(strSQL, None, "") pFeatureLayer = None pFeatureDataset = None
我们来看一下时间对比情况（影像大小为9G）：
| 实现平台 | 时间 |
| GDAL 栅格转矢量API | 1小时43分钟 |
| 改进后API | 3分45秒 |
当然，我没有用高性能服务器跑，用的是家里的一台普通台式机，CPU为AMD 3600，（这里我要插一句，AMD YES！！！干死intel）。12个核心全部跑起来，执行效率真的是相当惊人。
好吧，没图我说什么(这是南京地区的高分辨率二值图影像，感谢南京朋友提供，我也忘记他的名字了。。)：
图 1 转换后的矢量
从这个来看，结果还是不错的。如果想技术交流或源码，QQ：1044625113，请备注超大影像处理。
最后，我想说几句心里话，我其实挺讨厌写那种狗屁论文的，啥用没有，就在哪里讲故事，吹牛逼，我之前也干过这种事，包括现在也在干这种事，唉，没办法，其实很多其他东西远远比
发论文重要的多，比如这种超大影像转矢量问题，其实还有很多其他工程问题，有待我们去解决的。比如沈阳的李国春老师做的事情，其实我很欣赏的，一句话来总结我的观点：“talk is cheap， show me your code.”