土地报备坐标txt文件转shp遇到的坑以及该功能的 Python(Arcpy) 实现 / 开普饭

前言：在各种项目中，将国土报备坐标 txt 文本转换为 shp 是很常见的，一般都是通过 ArcGIS 添加XY数据的方式来实现。但是效率较低，同时存在着导出的 shp 线段出错、到处飞线的情况。所以这里结合 Python 不仅能快速将 txt 文件转换为 shp 文件，同时还能解决遇到的各种奇怪问题...一. 使用 Python(ArcPy) 绘制shp

1.什么是ArcPyArcPy 是一个安装 ArcGIS 会附带的站点包，通过 Python 实现。简言之，能通过 Python 直接调用 arcpy 执行地理数据分析、数据转换、数据管理和地图自动化以及多种客制化的需求。ArcGIS 帮助链接：https://desktop.arcgis.com/zh-cn/arcmap/10.3/analyze/arcpy/what-is-arcpy-.htm需要注意的是 ArcPy 是通过 Python2 实现的，不兼容 Python3！需要注意的是 ArcPy 是通过 Python2 实现的，不兼容 Python3！需要注意的是 ArcPy 是通过 Python2 实现的，不兼容 Python3！除非使用 ArcGIS pro，要花钱哦！2.如何构造shp(面)主要涉及类Array 类ArcPy 中自带的数组对象，可包含点和数组，用于构造几何对象。官方文档：https://desktop.arcgis.com/zh-cn/arcmap/10.3/analyze/arcpy-classes/array.htmPolygon 类可以使用 Polygon 类从头开始创建几何对象。接收 Array 类中数据然后使用 Polygon 类完成要素类的创建。官方文档：https://desktop.arcgis.com/zh-cn/arcmap/10.3/analyze/arcpy-classes/polygon.htm构造面点连成线，线组成了一个面。只需在 arcpy.Polygon() 类中传入一个点集就能形成一个面，点集（数组）通过 arcpy.Array() 构造传入（最后一个点不必与首个点相同，就算没有也会自动闭合）。[[20.0, 20.0],[30.0, 20.0],[30.0, 10.0],[20.0, 10.0]]

多个面就是多个点集。两个点集就是两个面、两个图形。[[[5.0, 3.0],[3.0, 3.0],[3.0, 5.0],[5.0, 3.0]],[[7.0, 5.0],[5.0, 5.0],[5.0, 7.0],[7.0, 5.0]]]

那么如果是环岛、孔洞等图形呢？不过是多个点集，当某个点集在另一个点集范围内时，arcpy.Polygon() 类可以自动把重叠部分完全扣除掉，实现环岛、孔洞等。所以使用该方法绘制图形是非常方便的。[[[40.0, 40.0], [50.0, 40.0], [50.0, 30.0],[40.0, 30.0]],[[45.0, 35.0], [48.0, 35.0], [48.0, 36.0],[45.0, 36.0]]]

以下是使用 arcpy.Array() 和 arcpy.Polygon() 实现简单的点集转 shp 代码示例：# -*- coding:utf-8 -*-# ---------------------------------------------------------------------------# Author: Lcc hygnic# Created on: 2020/10/16 15:55# Reference:"""Description:Python2.7Usage:"""# ---------------------------------------------------------------------------import arcpyimport osdef draw_poly(coord_list, sr, y, x):"""创建多边形coord_list(List)：多个点组成的坐标sr: 投影系y(Int): y坐标列x(Int): x坐标列"""parts = arcpy.Array()yuans = arcpy.Array()yuan = arcpy.Array()for part in coord_list:for pnt in part:if pnt:yuan.add(arcpy.Point(pnt[y], pnt[x]))else:# null point - we are at the start of a new ringyuans.add(yuan)yuan.removeAll()# we have our last ring, add ityuans.add(yuan)yuan.removeAll()# if we only have one ring: remove nestingif len(yuans) == 1:yuans = yuans.getObject(0)parts.add(yuans)yuans.removeAll()# 只有一个，单个图形if len(parts) == 1:parts = parts.getObject(0)return arcpy.Polygon(parts, sr)if __name__ == '__main__':# 点集，将之前的示例的三个图形放到一起poly = [[[20.0, 20.0], [30.0, 20.0],[30.0, 10.0], [20.0, 10.0]],[[5.0, 3.0], [3.0, 3.0],[3.0, 5.0], [5.0, 3.0]],[[7.0, 5.0], [5.0, 5.0],[5.0, 7.0], [7.0, 5.0]],[[40.0, 40.0], [50.0, 40.0],[50.0, 30.0], [40.0, 30.0]],[[45.0, 35.0], [48.0, 35.0],[48.0, 36.0], [45.0, 36.0]]]arcpy.env.overwriteOutput = Trueoutput_folder = os.getcwd()# 创建空白 shpname = "NewShp"# ▶注释1◀blank_shp = arcpy.CreateFeatureclass_management(output_folder, name, "Polygon", spatial_reference=None)# 创建、写入面Rows = arcpy.da.InsertCursor(blank_shp, "SHAPE@")f = polyp = draw_poly(f, sr=None, y=0, x=1)Rows.insertRow([p])write_sth = "Output path: " + os.path.join(output_folder, name)print write_sthdel Rows运行该程序会创建一个名为 NewShp.shp 的文件。▶注释1◀：arcpy.CreateFeatureclass_management() 创建一个空要素类，该方法可以传入参考系参数 spatial_reference。在这里没有设置参考系，如果之后需要的话，需要使用定义投影工具定义投影系或者提前就设置好投影系。二. 土地报备坐标txt文件（坐标交换数据）

数据主要依据标准《勘测定界界址点坐标交换格式》📚我们拿到的 TXT 文件中的数据基本按照以下标准：附录一：《勘测定界界址点坐标交换格式》...坐标交换格式具有两种格式，分别如下：文本格式[属性描述]格式版本号=数据产生单位=数据产生日期=坐标系=几度分带=投影类型=计量单位=带号=精度=转换参数=X平移,Y平移,Z平移,X旋转,Y旋转,Z旋转,尺度参数[地块坐标]界址点数,地块面积,地块编号,地块名称,记录图形属性(点、线、面),图幅号,地块用途,地类编码,@{点号,地块圈号,X坐标,Y坐标......点号,地块圈号,X坐标,Y坐标}界址点数,地块面积,地块编号,地块名称,记录图形属性(点、线、面),图幅号,地块用途,地类编码,@{点号,地块圈号,X坐标,Y坐标......点号,地块圈号,X坐标,Y坐标}注意：所有的逗号分隔符都必须是英文输入法状态下的逗号；地块圈号不能小于零；数据产生日期的格式为：2000-12-12；坐标系为54北京坐标系或80国家大地坐标系；投影类型为高斯克吕格或等角多圆锥；几度分带为3或6；带号、精度、转换参数、界址点数、地块面积、地块圈号,X坐标,Y坐标必须为数字型；且不能用该（9999,000,000）方式表示；地块编号、地块名称、记录图形属性(点、线、面)、图幅号、地块用途、地类编码、点号的每项里不能含有“,” 、“@”符号。点号中不能有字母，如J01等，建议为1。（预审报盘坐标格式要求）...例子：（该例子中数据经过处理，如果使用该数据生成shp,结果可能不符合预期）[属性描述]格式版本号=数据产生单位=国土资源部数据产生日期=2020-9-10坐标系=2000国家大地坐标系几度分带=6投影类型=高斯克吕格计量单位=米带号=35精度=0.001转换参数=0,0,0,0,0,0,0[地块坐标]8137,559.4952,1,火星乡地球村改造项目,面,,,,@J1,1,3433398.731421,35572460.011417J2,1,3433413.686436,35572361.307118J3,1,3433422.659303,35572214.746216J4,1,3433425.650143,35572095.104671J5,1,3433431.631890,35571897.696121J6,1,3433446.586915,35571804.9739J7,1,3433437.613643,35571673.368316J8,1,3433377.792817,35571493.906440J1,1,3433398.731421,35572460.01141718400,998.7312,1,火星乡地球村改造项目,面,,,,@J1,1,3437490.7824,35577811.1468J2,1,3437488.2225,35577838.46J3,1,3437488.2225,35577864.9191J4,1,3437503.5857,35577884.5502J5,1,3437514.6817,35577899.9133J6,1,3437514.6818,35577928.0804J7,1,3437506.9998,35577940.0295J8,1,3437506.5118,35577944.1285J9,1,3437503.9028,35577960.8496J1,1,3437490.7824,35577811.1468数据结构说明：

三.转换过程中的坑

那么把土地报备坐标 txt 文件转换成 shp，主要做的就是把TXT文件中的点正确的分成一个个点集（数组），然后传入 arcpy.Polygon() 类就行了。然而在实际操作中，由于各种原因，有很多坑，故可能出现下面这些情况：

这种情况出现的原因通常只有一个，就是没有正确的把一个个点集分开，导致本应该闭合的一个图形没有结束，然后最后一个点又连接到了下一个图形，出现这种线“乱飞乱连”的错误情况。而为什么没有正确把点集分开，是因为对数据的理解不全面（或者说数据有几种样式，txt 文本自己本人从相关官网下载，不存在其他人私自修改的可能），上文所说的标准是我在网上找到的，感觉说的不甚清楚且不全。实际在实践中，出现了以下几种样式（情况）：情况一无地块描述，如 18400,998.7312,1,火星乡地球村改造项目,面,,,,@ 。这种用 @ 符号结尾的段落，是每一个单独地块都有的描述符，本可以作为地块之间区分的标志物。然而实际上部分 txt 文件中直接就没有地块描述符 @ ，就如下图两个单独地块之间没有用地块描述符 @ 分开。

情况二序号混乱。下图的这种情况可以理解为多部件（一个地块实际是由两个单独地块组成的），就像这样

；如果第二个部件在第一个部件范围内，就会变成孔洞、环岛，就像这样

。从 txt 文本可以看到部件从1变成2，序号是连续的，从J1019 到 J1020（中间的 J1 是部件1的起始点，可有可无，在这里不影响）。

然后奇怪的事情来了，或者说，“另一种情况吧”。下图也是多部件，部件从1变成2，但是在这个 txt 文本中，序号没有紧接上面从 J6077 到 J6078，而是变成了 J1，然后 J2，重新开始了。

情况三更多的奇怪情况。从上面来看，由于序号规则不详、不完善或者说数据错误，非常容易混淆多部件和地块之间的编号。更糟糕的是，我还遇到过以下几种情况，正如下面三个不同 txt 文本的截图所示:从左到右，第一个是直接开头序号都不为1了，直接从 J9677 开始；第二个是部件开头部不为1，为0；第三个是部件计数跳跃，直接从1变成174，完全不知道为什么会这样！

四.解决思路 🤔

那么如何处理这些情况呢？咋一看上去非常复杂。实际上非常简单，那就是不去纠结部件和地块序号，直接全部当成不同的地块处理（相当于有多部件的地块被打散了，这完全不影响出图，反正在一张shp文件上）。到了这里，我已经逐渐理解一切，关键在于序号列和部件列，一旦序号列从1重新开始或者部件发生改变，数据就进行分割。分好的数据（一个地块）是一个列表，然后一个 txt 文件中的所有数据组成一个大列表进行处理。五.实现代码

完整实现代码（脚本形式）: ✌️好耶# -*- coding:utf-8 -*-# ---------------------------------------------------------------------------# arcgis 10.3# Author: hygnic lcc# Created on: 2020/4/8 16:33# Reference:"""Description: 将国土土地报备坐标txt文本处理生成shp文件原始数据如下： [属性描述] 格式版本号= 数据产生单位= 数据产生日期=2077 坐标系=2000ff 几度分带= 投影类型=高斯克吕格计量单位=米带号=ffff 精度=ffff 转换参数=fffffff [地块坐标] 41,0,1,0,面,0,0,,@ - J1,1,3.0, 8.0 J2,1,1.0, 8.0 J3,1,2.0, 10.0 J4,1,3.0, 8.0 555,0,1,0,面,0,0,,@ J1,1,9.0, 11.0 J2,1,9.0, 8.0 J3,1,6.0, 8.0 J4,1,6.0, 11.0 J1,2,7.0, 10.0 - J2,2,7.0, 9.0 - J3,2,8.0, 9.0 - J4,2,8.0, 10.0处理为如下列表： [ [ [3.0, 8.0], [1.0, 8.0], [2.0, 10.0], [3.0, 8.0] ] , [ [9.0, 11.0], [9.0, 8.0], [6.0, 8.0], [6.0, 11.0], [9.0, 11.0] ] , [ [7.0, 10.0], [7.0, 9.0], [8.0, 9.0], [8.0, 10.0], [7.0, 10.0] ] ]每一个会闭合的线都组成一个单独的列表，防止因为首尾不接导致闭合面错误Usage:"""# ---------------------------------------------------------------------------import reimport arcpyimport osdef points_genarator(txt_file): """ points_genarator(txt_file) return list txt_file：文本文件地址将txt转换为可以使用的点集列表 """ fffs = open(txt_file, "r") input_data = fffs.readlines() # input_data.remove("\n") # 去除换行符 # input_data = [x.strip() for x in input_data] # input_data = input_data[2:] # 去除前13行 #▶注释1◀ 去除带@的行 input_data = [x.strip() for x in input_data if "@" not in x] # 去除非J行 while True: # 去除前13行 # 如果首行字符不在下列元组中 first = ("J", "1", "0", "2", "3", "4", "5", "6", "7", "8", "9") if input_data[0][0] not in first: del input_data[0] # elif : # del input_data[0] else: break # input_data = input_data[12:] # 去除前13行 fffs.close() # 每一根闭合线单独组成一个列表 line_closed = [] # 多根线条组成面 polygon_list = [] #▶注释2◀ while input_data: row = input_data.pop(0) # ['J1', '1', '3405133.8969', '35353662.0113'] row1 = row.split(",") # 去除字母和数字组合中的字母如 "J1" 只保留 "1" row_num = re.findall(r'[0-9]+|[a-z]+', row1[0]) # ['1'] row_num2 = row1[1] # '1' 文本最后有空行会报错 if not line_closed: # 为空时 line_closed.append(row1) flag1 = int(row_num[0]) # 1 第一列数字 # 第二列数字 flag2 = row1[1] # '1' elif int(row_num[0]) > flag1: if row_num2 != flag2: # 第二列出现不一样的，新一轮开始 flag2 = row_num2 polygon_list.append(line_closed) line_closed = [] # 未开始新一轮就接着上一个，如果开始新一轮那么这里就是第一个起始点 line_closed.append(row1) elif int(row_num[0]) == flag1: # 新一轮开始 polygon_list.append(line_closed) line_closed = [] line_closed.append(row1) if line_closed: polygon_list.append(line_closed) return polygon_listdef draw_poly(coord_list, sr, y, x): """ 创建多边形 coord_list(List)：多个点组成的坐标 sr: 投影系 y(Int): y坐标列 x(Int): x坐标列 """ parts = arcpy.Array() yuans = arcpy.Array() yuan = arcpy.Array() for part in coord_list: for pnt in part: if pnt: yuan.add(arcpy.Point(pnt[y], pnt[x])) else: # null point - we are at the start of a new ring yuans.add(yuan) yuan.removeAll() # we have our last ring, add it yuans.add(yuan) yuan.removeAll() # if we only have one ring: remove nesting if len(yuans) == 1: yuans = yuans.getObject(0) parts.add(yuans) yuans.removeAll() # 只有一个，单个图形 if len(parts) == 1: parts = parts.getObject(0) return arcpy.Polygon(parts, sr)def main(info, txt_folder, output_folder): """ 功能实现的主函数 info(List): 存放信息的列表，作为独立脚本使用时没啥用 txt_folder(Unicode/String): 包含txt文件的文件夹 output_folder(Unicode/String): 导出文件夹 """ arcpy.env.overwriteOutput = True txts = os.listdir(txt_folder) for txt in txts: if txt[-3:].lower() == "txt": txt_p = os.path.join(txt_folder, txt) f = points_genarator(txt_p) name = os.path.splitext(os.path.basename(txt_p))[0] # 创建空白shp blank_shp = arcpy.CreateFeatureclass_management( output_folder,name, "Polygon", spatial_reference=None) # create the polygons and write them Rows = arcpy.da.InsertCursor(blank_shp, "SHAPE@") #▶注释3◀ p = draw_poly(f, sr=None, y=3, x=2) Rows.insertRow([p]) del Rows write_sth = "Export succeed: " + os.path.join( output_folder, name) print write_sth info.append(write_sth) # 给新建的shp文件添加 MC 字段和值 shp_name = name + ".shp" newfield_name = "MC" fresh_layer = arcpy.mapping.Layer( os.path.join(output_folder, shp_name)) arcpy.AddField_management( fresh_layer, newfield_name, "TEXT", field_length=100) cursor2 = arcpy.da.UpdateCursor(fresh_layer, newfield_name) for roww in cursor2: roww[0] = name cursor2.updateRow(roww) del cursor2if __name__ == '__main__': """脚本单独使用""" """----------------------------------------------""" """---------------------PARA---------------------""" arcpy.env.overwriteOutput = True #▶注释4◀ dir_p = ur"G:\高标准分布图\xx县\测试" # txt文件夹 output_dir = ur"G:\高标准分布图\xx县\文件夹"# 输出文件夹 """----------------------------------------------""" """----------------------------------------------""" infos = [] main(infos, dir_p, output_dir)points_genarator 方法将输入的 txt 文件正确分割成多个点集（数组）；然后使用 draw_poly 方法将该点集绘制成图形。▶注释1◀：由于我们不再详细区分多部件等，所以不再需要描述符 @，况且有的 txt 文件直接都没有该描述符。所以直接将存在 @ 的行删除。▶注释2◀：这个循环中的代码看上去比较复杂，其作用是判断每一行的部件或者序号的顺序是否发生改变，从来将不同点集直接分开。▶注释3◀：文本中的第三行和第四行对应的不是 x y 坐标，而是 y x 坐标。所以 y=3, x=2。▶注释4◀：输入 txt 文件所在地址和输出的 shp 文件的地址，然后运行即可。Note: 之前有的朋友问我为什么运行报错，我检查了很久，最后发现 txt 文件最后有空行。所以如果有报错的话可以先检查一下是否有空行，一般从自然资源局等官网下载下来的 txt 文件是没有的空行的。正确分割数据得到的正确 shp 文件与错误 shp 文件的对比：

结束语

参考文章关于构造面：https://gis.stackexchange.com/questions/14952/arcgis-10-0-arcpy-how-to-create-a-polygon-geometry-from-inner-and-outer-ring-po使用版本ArcGIS10.3 Python2.7源代码及其PDF文档下载地址：链接：https://pan.baidu.com/s/1ihBKcufLuRkAmXlJFFsDDA 提取码：3t9e如果下载地址失效的话，可以回复 d1 获取最新下载地址。分享GIS，不止于Python。很难的教程、有趣的文章、好看的地图，这里都有！荟GIS精粹，关注公众号：GIS荟，带你飞!很难的文章系列：《Python和ArcGIS自动化制图完全指南》——共六章，附带教程源码、数据和文档手册下载。《深入制图表达》——深入制图表达的机制和运用实现。《ArcGIS制图高级技巧》——ArcGIS制图高级技巧，发现 ArcGIS 不为人知的各种制图美化技巧。（持续更新中！）《基于Python的ArcGIS(ArcPy)多进程自动出图》——使用多进程快速大批量出图。《用Python创建你第一个GIS程序》——简单易懂，手把手教你用 Python 搭建你的第一的 GIS 程序。《GIS 进阶成神》 ——进阶成神之路！...不难的有趣文章系列：《从地图发现世界》——从地图，发现奇特的、美丽的、我的世界。（持续更新中！）《制图艺术》——制作不出优美的地图？进来看看是不是缺点东西（持续更新中！）...更多文章可以使用搜索哦欢迎交流

土地报备坐标txt文件转shp遇到的坑以及该功能的 Python(Arcpy) 实现

相关推荐