用地CAD转GIS一直都是老大难的问题,主要办法是通过FME等工具。
GIS中读取的CAD是分为点、线、面几个图层,与GSI的数据集分类是一致的,这个里面并没有填充面。
基于ArcGIS的转换有两个思路,一是读取dxf文件中的hatch信息,然后在GIS中创建面。
二是通过GIS打开DWG,读取其中面相关的信息,创建面。
读取dxf文件
DXF是AutoCAD 绘图交换文件。DXF 是Autodesk(欧特克)公司开发的用于AutoCAD与其它软件之间进行CAD数据交换的CAD数据文件格式。
DXF文件由标题段、表段、块段、实体段和文件结束段5部分组成,好处是直接用python可以读取其中的信息,只是格式比较复杂有很多搞不懂的数据。
看官方的dxf帮助文件 ,这个是组码代码,搞懂了组码基本上可以获取信息了。
主要需要获取的几个组码如下
8是图层名称,需要获取hatch的图层名称得到用地类型
10、20就是每个折点的XY坐标,也可以是圆心之类的坐标
62是颜色,是CAD的色号,
CAD最坑爹的就是会把RGB颜色转换成近似色号储存在62里,直接导入GIS的图层字段color就是这个颜色号
真正的RGB色彩在这里420里面,将RRRGGGBBB的数字转成16进制的RRGGBB,然后再整体转成十进制的数字,藏的好深。
接下来看代码
读取dxf文件
dxf_file = 'dxf文件路径'
with open(dxf_file, mode="r") as f: #,encoding='utf-8'
cad = f.read()
data= cad.split('\n')#按行拆开
如果在python3里运行,建议加上encoding=‘utf-8’,这样中文可以正常显示,在GIS的python2.7里,会有时会出现编码错误。
然后读取hatch,找到hatch就找到了段落,看完会发现,这里没有颜色的62或者420,对的,如果颜色是bylayer那图元里就没有颜色标记了,得去layer段里找颜色。
图层信息是这样的
先获取layer的所有颜色
def read_layer(data):
layers = {}
for i in range(len(data)):
if data[i] == 'LAYER' and data[i-1]==' 0': #检查前一个0,保证上一个结束
#print(data[i],i)
for j in range(i,len(data)): #检查LAYER段
if data[j] == ' 2': #图层
layer = data[j+1]
elif data[j]==' 62': #色号
color = str(int(data[j+1]))
elif data[j]=='420': #有RGB用RGB,420都是后出现的
color = to_RGB(data[j+1])
elif data[j]==' 0':
layers[layer] = color
break
return layers
这里需要颜色转换,10进制整体转16进制,再两两转10进制就是RGB
def to_RGB(color):
c_16 = hex(int(color, 10))[2:]#整体转16进制,去掉开头的0x
return '%s,%s,%s'%(int(c_16[:2],16),int(c_16[2:4],16),int(c_16[4:6],16))
#两两16转10进制,成为RGB色号
最后是把hatch的数据读出来
#读取dxf中的hatch属性
def read_hatch(data):
hatchs = []
for i in range(len(data)):
if data[i] == 'HATCH': #找到hatch
#print(data[i],i)
rings=[]
layer,color = '',''
for j in range(i,len(data)):
if data[j] == ' 8':#图层
layer = data[j+1]
elif data[j]== ' 93' : #坐标开始
for a in range(j,len(data)):
if data[a]==' 10' and data[a+2] ==' 20' and data[a+1]!='0.0' and data[a+3]!='0.0':
rings.append([data[a+1],data[a+3]])
elif data[a]==' 97' and data[a+1]==' 0':
break
elif data[j]=='420':
color = to_RGB(data[j+1])
elif data[j]==' 62':
color = data[j+1]
elif data[j]==' 0':
if color == '':
color = layers[layer]
hatchs.append([layer,color,rings])
break
return hatchs
但测试之后发现还有一些特殊情况,线段边的是72 1 开头的,而曲线,由72 2 开头的,圆心位置的XY坐标显著偏差。这个不能用正常的折线点坐标的方式来记录。
生成边界看一下确实是有一段很短的弧线!!
目前在GIS中没有找到arcpy生成弧线的办法,所以,这个方式中断了
在GIS读取含弧线面的坐标,结果出来了一串坐标点,也就是说GIS中弧线是以多线段表达的。
下篇写另一个方式