超大影像栅格转矢量快速实现
距离上一次博客更新,起码又是大半年,时光飞逝,我也已经老了。。。这一次,我解决了一个工程上的小问题,可能在行家看来简单,但是呢,它好像又没那么简单,就是我们通常用的栅格转矢量,
我们知道栅格转矢量,通常有以下方法:采用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 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) 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.”