近期由于工作原因需要读写shape文件，主要涉及几何数据和属性数据，在网上学习了一番，有些问题网上也没解决方案，不过最终自己还是解决了，今天有空在此作一番整理，希望与各位进行分享
一、前期准备
python读写shp文件所需的库为pyshp，对应的安装命令为 pip install pyshp
注：arcgis 10.2 及以下版本的shp文件是gbk编码，10.2.1及以上版本是utf8编码，pyshp保存的文件是utf8编码，中文乱码由编码不一致导致
二、写出shape文件
本来想先写读取shape文件的，但是手头没有测试用的shp数据，就决定先写写出shp数据吧，这样就可以造一个读取用的测试数据了
直接上代码吧，通过代码来介绍会比较方便
1、需要导入库，import shapefile
2、定义一个写出shp文件的函数 write_shp，包含众多参数，下面一个个介绍
3、target和shape_type是初始化shape.Writer的两个必要参数
target是生成的shape文件的路径，具体到文件名即可，不用带后缀，程序会自动生产shp、dbf、shx三个文件
shape_type是指定几何数据类型，是一个INT类型的数据，取值范围见代码，我们主要用1、3、5，分别对应点、线、面
autoBalance设置为1可以避免几何数据行数与属性行数不匹配的问题，此时系统会自动填充属性，保持数量对应，默认是False，需要设置
4、geo_data为几何数据，点线面的结构不同，在本函数中：
points为点集（不是几何类型的点集，是我们包含很多行记录，多行记录组成的点集，下同），用列表表示，列表中的每一项为一个二维的点坐标，同样为列表类型，包含两个元素，分别为xy坐标，一个坐标对应一行记录；
polylines为线集，用列表表示，列表中的每一项为一条线，一条线对应一行记录；一条线也用列表表示，其中的每一项为线的每一段（一条线可能分几段，一般就一段）；一段也用列表表示，其中的每一项为一个二维点，同point；
polygons为面集，用列表表示，列表中的每一项为一个面，实际上也是一条线，只是首尾相连；同线有好几段，面也可能有几条线组成，其他类似，所以线和面的结构相同；
pointzs为三维点集，相比于points多了Z坐标，polylinezs和polygonzs也相同。
5、后面五个参数是属性相关参数：
field_names为属性名数组，个数与属性个数一致
field_types为属性的数据类型数组，个数与属性个数一致，具体见代码，主要用“C”字符类型，“N”数字类型，“F”浮点数类型
field_sizes为属性的数据长度，个数与属性个数一致，这个参数一般用于“C”字符类型的设置，默认为50
field_decimals为属性值的保留小数位数，个数与属性个数一致，这个参数一般用于“F”浮点型的设置，默认为0
field_values为属性值数组，个数与geo_data(行数)一致；每一项为一条记录，数值顺序为属性的顺序，个数与geo_data(行数)一致，属性值不存在的用None表示
6、最后说一下文件的保存，之前看网上有用save()函数，但是我调用之后发现无此函数，同样查看源码后发现save()函数被注释了，写出文件已经集成到了record()函数，所以最后不用调save()
#coding=utf-8
import shapefile
def write_shp(target, shape_type, geo_data, field_names, field_types, field_sizes, field_decimals, field_values):
"""
生成的文件名，自动生成.shp .dbf .shx
target = "D:\\test"
'''
几何类型
NULL = 0
POINT = 1
POLYLINE = 3
POLYGON = 5
MULTIPOINT = 8
POINTZ = 11
POLYLINEZ = 13
POLYGONZ = 15
MULTIPOINTZ = 18
POINTM = 21
POLYLINEM = 23
POLYGONM = 25
MULTIPOINTM = 28
MULTIPATCH = 31
'''
shapeType = 1
几何数据
points = [[1,1],[2,2]]
polylines = [[[[1,1],[2,2]],[[2,2],[3,3]]],[[[1,2],[2,4]]]]
polygons = [[[[1,1],[2,2],[3,1]],[[3,1],[2,2],[3,2]]],[[[1,2],[2,3],[1,3]]]]
pointzs = [[1,1,1],[2,2,2]]
polylinezs = [[[[1,1,1],[2,2,2]],[[2,2,3],[3,3,4]]],[[[1,2,5],[2,4,6]]]]
polygonzs = [[[[1,1,1],[2,2,2],[3,1,3]],[[3,1,4],[2,2,5],[3,2,6]]],[[[1,2,7],[2,3,8],[1,3,9]]]]
属性名称
field_names = ['test1','test2','test3']
'''
属性类型
“C”：字符，文字。
“N”：数字，带或不带小数。
“F”：浮动（与“N”相同）。
“L”：逻辑，表示布尔值True / False值。
“D”：日期。
“M”：备忘录，在GIS中没有意义，而是xbase规范的一部分。
'''
field_types = ['C','N','F']
字段长度，默认50
field_sizes = [50,50,50]
字段精度，默认为0
field_decimals = [0,0,2]
属性值
field_values = [['test',1,1.111],['测试',2,2.226]]
"""
try:
w=shapefile.Writer(target=target, shapeType=shape_type, autoBalance=1)
# 添加属性
for idx_field in range(len(field_names)):
w.field(field_names[idx_field], field_types[idx_field], field_sizes[idx_field], field_decimals[idx_field])
for idx_record in range(len(geo_data)):
# 写出几何数据
if shape_type == 1:
w.point(*geo_data[idx_record])
elif shape_type == 3:
w.line(geo_data[idx_record])
elif shape_type == 5:
w.poly(geo_data[idx_record])
elif shape_type == 11:
w.pointz(*geo_data[idx_record])
elif shape_type == 13:
w.linez(geo_data[idx_record])
elif shape_type == 15:
w.polyz(geo_data[idx_record])
else:
continue
# 写出属性数据
record = []
for idx_field in range(len(field_names)):
record.append(field_values[idx_record][idx_field])
w.record(*record)
return "success"
except shapefile.ShapefileException as e:
return repr(e)
field_names = ['test1','test2','test3']
field_types = ['C','N','F']
field_sizes = [50,50,50]
field_decimals = [0,0,2]
field_values = [['test',1,1.111],['测试',2,2.226]]
target = "D:\\data\\shp\\point"
shape_type = 1
points = [[1,1],[2,2]]
print(write_shp(target, shape_type, points, field_names, field_types, field_sizes, field_decimals, field_values))
target = "D:\\data\\shp\\polyline"
shape_type = 3
polylines = [[[[1,1],[2,2]],[[2,2],[3,3]]],[[[1,2],[2,4]]]]
print(write_shp(target, shape_type, polylines, field_names, field_types, field_sizes, field_decimals, field_values))
target = "D:\\data\\shp\\polygon"
shape_type = 5
polygons = [[[[1,1],[2,2],[3,1]],[[3,1],[2,2],[3,2]]],[[[1,2],[2,3],[1,3]]]]
print(write_shp(target, shape_type, polygons, field_names, field_types, field_sizes, field_decimals, field_values))
target = "D:\\data\\shp\\pointz"
shape_type = 11
pointzs = [[1,1,1],[2,2,2]]
print(write_shp(target, shape_type, pointzs, field_names, field_types, field_sizes, field_decimals, field_values))
target = "D:\\data\\shp\\polylinez"
shape_type = 13
polylinezs = [[[[1,1,1],[2,2,2]],[[2,2,3],[3,3,4]]],[[[1,2,5],[2,4,6]]]]
print(write_shp(target, shape_type, polylinezs, field_names, field_types, field_sizes, field_decimals, field_values))
target = "D:\\data\\shp\\polygonz"
shape_type = 15
polygonzs = [[[[1,1,1],[2,2,2],[3,1,3]],[[3,1,4],[2,2,5],[3,2,6]]],[[[1,2,7],[2,3,8],[1,3,9]]]]
print(write_shp(target, shape_type, polygonzs, field_names, field_types, field_sizes, field_decimals, field_values))
函数和测试代码都如上所示，可直接运行，得到对应的文件，如下图所示
使用ArcMap打开二维shp，如下图所示
几何数据显示正常，接口调用的参数是polygons，读者也可以改成points或polylines
属性也显示正常，因为field_decimals设置的浮点数保留2位小数，所以输入的3位小数四舍五入成2位
使用ArcScene打开三维shp，如下图所示
三、读取shp文件
通过上一节，我们得到了带有几何数据和属性数据的shp文件，这一节我们来读取它们
3.1 读取属性列表
首先来读取属性列表，直接上代码
#coding=utf-8
import shapefile
def get_shp_field_list(path):
try:
try:
file = shapefile.Reader(path)
except UnicodeDecodeError:
file = shapefile.Reader(path, encoding="gbk")
fields = file.fields
res = []
for field in fields:
res.append(field[0])
return res
except shapefile.ShapefileException as e:
return repr(e)
print(get_shp_field_list("D:\\data\\shp\\polygon"))
path是shp文件的路径，读取使用的是shapefile.Reader，第一个参数是shp路径，只需写到文件名即可，不用带后缀（带也无所谓），encoding需要设置为‘gbk’，用来解析中文
调用此函数，可得到
3.2 读取几何和属性值
获取到属性列表后，可继续获取每一行的几何数据和你想要的属性的属性值，代码如下。
path是shp文件路径，同上
fields是需要获取的属性名列表
#coding=utf-8
import shapefile
def get_shp_shape_records(path, fields):
try:
try:
file = shapefile.Reader(path)
shape_records = file.shapeRecords()
except UnicodeDecodeError:
file = shapefile.Reader(path, encoding="gbk")
shape_records = file.shapeRecords()
shp_types = []
shp_datas = []
shp_fields = []
for i in range(len(fields)):
field_values = []
shp_fields.append(field_values)
for shape_record in shape_records:
# type
shp_type = shape_record.shape.shapeType
if shp_type == 1:
shp_types.append('point')
points = shape_record.shape.points
points_order = []
points_order.append(len(points))
for point in points:
points_order.append(point[0])
points_order.append(point[1])
shp_datas.append(points_order)
elif shp_type == 3:
shp_types.append('polyline')
points = shape_record.shape.points
parts = shape_record.shape.parts
polyline_record = []
for idx in range(0,len(parts)-1):
polyline_record_part = points[parts[idx]:parts[idx+1]]
polyline_record.append(polyline_record_part)
polyline_record_part = points[parts[-1]:]
polyline_record.append(polyline_record_part)
shp_datas.append(polyline_record)
elif shp_type == 5:
shp_types.append('polygon')
points = shape_record.shape.points
parts = shape_record.shape.parts
polygon_record = []
for idx in range(0,len(parts)-1):
polygon_record_part = points[parts[idx]:parts[idx+1]]
polygon_record.append(polygon_record_part)
polygon_record_part = points[parts[-1]:]
polygon_record.append(polygon_record_part)
shp_datas.append(polygon_record)
elif shp_type == 11:
shp_types.append('pointz')
points = shape_record.shape.points
zs = shape_record.shape.z
points_order = []
points_order.append(len(points))
for idx in range(len(points)):
points_order.append(points[idx][0])
points_order.append(points[idx][1])
points_order.append(zs[idx])
shp_datas.append(points_order)
elif shp_type == 13:
shp_types.append('polylinez')
points = shape_record.shape.points
zs = shape_record.shape.z
parts = shape_record.shape.parts
polylinez_record = []
for idx in range(0,len(parts)-1):
polylinez_record_part = []
for idx_pt in range(parts[idx],parts[idx+1]):
pt = list(points[idx_pt])
pt.append(zs[idx_pt])
polylinez_record_part.append(pt)
polylinez_record.append(polylinez_record_part)
polylinez_record_part = []
for idx_pt in range(parts[-1],len(points)):
pt = list(points[idx_pt])
pt.append(zs[idx_pt])
polylinez_record_part.append(pt)
polylinez_record.append(polylinez_record_part)
shp_datas.append(polylinez_record)
elif shp_type == 15:
shp_types.append('polygonz')
points = shape_record.shape.points
zs = shape_record.shape.z
parts = shape_record.shape.parts
polygonz_record = []
for idx in range(0,len(parts)-1):
polygonz_record_part = []
for idx_pt in range(parts[idx],parts[idx+1]):
pt = list(points[idx_pt])
pt.append(zs[idx_pt])
polygonz_record_part.append(pt)
polygonz_record.append(polygonz_record_part)
polygonz_record_part = []
for idx_pt in range(parts[-1],len(points)):
pt = list(points[idx_pt])
pt.append(zs[idx_pt])
polygonz_record_part.append(pt)
polygonz_record.append(polygonz_record_part)
shp_datas.append(polygonz_record)
else:
return "undefined type"
# field
for idx, field in enumerate(fields):
shp_fields[idx].append(shape_record.record[field])
return shp_types, shp_datas, shp_fields
except shapefile.ShapefileException as e:
return repr(e)
print(get_shp_shape_records("D:\\data\\shp\\point", ['test1','test3']))
print(get_shp_shape_records("D:\\data\\shp\\polyline", ['test1','test3']))
print(get_shp_shape_records("D:\\data\\shp\\polygon", ['test1','test3']))
print(get_shp_shape_records("D:\\data\\shp\\pointz", ['test1','test3']))
print(get_shp_shape_records("D:\\data\\shp\\polylinez", ['test1','test3']))
print(get_shp_shape_records("D:\\data\\shp\\polygonz", ['test1','test3']))
调用此函数，可得到一个元组，共三项
第一项是每一行的几何类型
第二项是每一行的几何数据，结构同写出shp时的介绍
第三项是每一行的属性数据，列表中的元素数量与输入的属性名个数相同，每个属性的列表中的元素个数与几何数据的个数（行数）相同
四、shp文件追加记录
有时候我们需要为已经存在的shp文件增加记录，我翻看了shapefile库，未找到可以追加的方法，只能采用大家都能想得到的方法了，方法流程和代码如下
先用read的方法获取原shp文件数据
然后用write的方法创建新的shp文件对象，把read出来的数据写入该对象
最后将追加的数据写入该对象
#coding=utf-8
import shapefile
def add_shp_records(path, geo_data, field_values):
try:
# 先读取原文件的几何数据和属性数据
try:
file = shapefile.Reader(path)
shape_records = file.shapeRecords()
shape_fields = file.fields
shape_type = file.shapeType
except UnicodeDecodeError:
file = shapefile.Reader(path, encoding="gbk")
shape_records = file.shapeRecords()
shape_fields = file.fields
shape_type = file.shapeType
# 创建与原文件同名的导出文件
w=shapefile.Writer(target=path, shapeType=shape_type, autoBalance=1)
# 创建与原文件相同的属性
for shape_field in shape_fields:
w.field(*shape_field)
# 将原文件的几何数据和属性数据添加到新文件
for shape_record in shape_records:
w.shape(shape_record.shape)
w.record(*shape_record.record)
# 在新文件追加待添加的几何数据和属性数据
for idx_record in range(len(geo_data)):
# 追加几何数据
if shape_type == 1:
w.point(*geo_data[idx_record])
elif shape_type == 3:
w.line(geo_data[idx_record])
elif shape_type == 5:
w.poly(geo_data[idx_record])
elif shape_type == 11:
w.pointz(*geo_data[idx_record])
elif shape_type == 13:
w.linez(geo_data[idx_record])
elif shape_type == 15:
w.polyz(geo_data[idx_record])
else:
continue
# 追加属性数据
w.record(*field_values[idx_record])
return "success"
except shapefile.ShapefileException as e:
return repr(e)
path = 'D:\data\shp\point'
points = [[3,3],[4,4]]
field_values = [['add',3,None],['增加',4,None]]
print(add_shp_records(path, points, field_values))
原文件arcgis显示
追加记录后arcgis显示
扫一扫
3854
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
