统计矢量文件节点数的脚本

为什么80%的码农都做不了架构师？>>>   

用来做矢量化的工作量统计.昨天统计了一下,基本上使用绘图工具并需要做平移的建筑轮廓绘图,最慢的每分钟只有五个节点,这个有点太慢了,最快的是每分钟18个节点,这个速度其实也不是很快,因为毕竟是工具辅助.没有合理的工作量统计,每天就坐在那里不动,也看不出是不是在走神.这种状态可以结束了.

---

做了一点点更新，增加了字段填充率和图层范围的WKT的表示

#coding:utf-8
# 这个工具通过遍历指定路径下的矢量文件(shp)中的图元和图元要素数量来统计
# 20190216 完成第一版本
# 20190222 增加对字段内容填充数量和百分比的评估
from __future__ import division
import os
import sys
from osgeo import ogr
from osgeo import osr
os.environ["PYTHONIOENCODING"]="gbk"
# 从指定路径获取文件列表,文件夹递归,返回所有shape文件
def list_all_files(rootdir):
    import os
    _files = []
    list = os.listdir(rootdir) #列出文件夹下所有的目录与文件
    for i in range(0,len(list)):
           path = os.path.join(rootdir,list[i])
           if os.path.isdir(path):
              _files.extend(list_all_files(path))
           if os.path.isfile(path) and path.lower().endswith(".shp"):
              _files.append(path.decode('gbk').encode('utf-8'))
    return _files
# 递归计算图元的节点数,累加并返回
def getAllNodeNumber(geom):
	subgeomNum=geom.GetGeometryCount()
	totalNodeNum=0
	if subgeomNum==0:
		return geom.GetPointCount()
	elif subgeomNum>0:
		for i in range(subgeomNum):
			totalNodeNum+=getAllNodeNumber(geom.GetGeometryRef(i))
		return totalNodeNum

# 统计给定参数文件名列表中的每个图层的信息并返回统计结果
def countFile(vecFileName):
	driver = ogr.GetDriverByName('ESRI Shapefile')
	# driver = ogr.GetDriverByName('MapInfo File')
	dataSource = driver.Open(vecFileName) # 0 means read-only. 1 means writeable.
	geoTypeName=["Unknown","POINT","LINESTRING","POLYGON","MULTIPOINT"	]
	if dataSource is None:
		print 'Could not open %s' % (vecFileName)
	else:
		# print 'Opened %s' % (daShapefile)
		layer = dataSource.GetLayer()
		fieldsNumber=layer.GetLayerDefn().GetFieldCount()
		featureCount = layer.GetFeatureCount()
		feat = layer.GetNextFeature() 
		layerType = geoTypeName[layer.GetGeomType()]
		filledField=0
		nodeNum=0
		evlNodeNumber=0
		count=0
		e=layer.GetExtent()
		target = osr.SpatialReference()
		target.ImportFromEPSG(4326)
		source=layer.GetSpatialRef()
		if(source==None):
			source= osr.SpatialReference()
			source.ImportFromEPSG(4326)
		transform = osr.CoordinateTransformation(source, target)
		lyrExtent=ogr.CreateGeometryFromWkt( "POLYGON ((%f %f, %f %f,%f %f,%f %f))" % (e[0],e[2],e[0],e[3],e[1],e[3],e[1],e[2] ),layer.GetSpatialRef())
		lyrExtent.Transform(transform)
		lyrExtentArea=lyrExtent.Area()
		# lyrExtentArea=ogr.CreateGeometryFromWkt( "POLYGON (%s)" % str(lyrExtent),layer.GetSpatialRef())
		for fi in range(featureCount):
			feat=layer.GetFeature(fi)
			for fn in range(fieldsNumber):
				if feat.GetFieldAsString(fn)!="":
					 filledField+=1
			geom= feat.GetGeometryRef()	
			if geom == None:				
				continue
			nodeNum+=getAllNodeNumber(geom)
			count+=1
		if featureCount!=0 :
			evlNodeNumber=nodeNum/featureCount
		del driver
		del dataSource
		del geoTypeName
		filledRate=0.0
		filledRate=filledField*100/(fieldsNumber*featureCount) if fieldsNumber>0 and featureCount>0 else 0
		evlNodepreArea=nodeNum/lyrExtentArea  if lyrExtentArea>0  else 0
		print "\"%s\",\"%s\",%d,%.1f%%,%d,%d,%d,\"%s\",\"%s\"" % (os.path.basename(vecFileName.decode('utf-8').encode('gbk')),layerType,fieldsNumber,filledRate,featureCount,nodeNum,evlNodeNumber,os.path.abspath(vecFileName.decode('utf-8').encode('gbk')),lyrExtent)


# 主过程


daShapefiles =[]
# 判断参数是否存在
if len(sys.argv)!=2:
	print len(sys.argv)
	print "a path is nessary"
	exit()
# 从参数中提取路径或者文件名
path=sys.argv[1]
# 判断参数内容获取文件列表
if os.path.exists(path) and os.path.isfile(path) and path.lower().endswith(".shp"):
	daShapefiles =[path]
elif os.path.exists(path) and os.path.isdir(path):
	daShapefiles =list_all_files(path)
# 测试是否需要进行输出
if len(daShapefiles)>0:
	print u"\"文件名\",\"图元类型\",\"字段数\",\"字段填充率\",\"总图元数\",\"总节点数\",\"图元平均节点数\",\"路径\",\"图层范围\"" 
else:
	print "No shape file was found"
# 遍历文件列表并输出信息
for daShapefile in daShapefiles:
	countFile(daShapefile)

del daShapefiles

转载于:https://my.oschina.net/rodger/blog/3010857