python2.x之pyshp的使用
pyshp是python读写shape文件的一个很简单的库。下面记录其用法:
用法详见代码中:
1 #! /usr/bin/env python 2 # -*- coding:utf-8 -*- 3 4 import shapefile 5 6 sf = shapefile.Reader("shapefile/d_map_1000000.shp") 7 shapes = sf.shapes() # shapes方法返回描述每个形状记录的几何形状的Shape对象的列表。 8 9 print(len(shapes)) # 2130 10 print (shapes[5]) # <shapefile._Shape instance at 0x03785530> 11 12 for name in dir(shapes[5]): 13 if not name.startswith('__'): 14 print(name) 15 # 这个对象有4个属性 16 # bbox 17 # parts 18 # points 19 # shapeType 20 21 22 print("shapes[5].bbox:",shapes[5].bbox) # bbox:如果形状类型包含多个点,则此元组描述左下角(x,y)坐标和右上角坐标, 23 # 在点周围创建一个完整的框。 如果shapeType是Null(shapeType == 0),则会引发AttributeError。 24 # ('shapes[5].bbox:', [-60.0, -88.0, -36.0, -84.0]) 25 26 27 print("shapes[5].parts:",shapes[5].parts) # parts:只需将点集合分组为形状。 如果形状记录具有多个部分,则该属性包含每个部分的第一点的索引。 28 # 如果只有一个部分,则返回包含0的列表。 29 # ('shapes[5].parts:', [0]) 30 31 32 print("shapes[5].points:",shapes[5].points) # points属性包含一个元组列表,其中包含形状中每个点的(x,y)坐标 33 # ('shapes[5].points:', [(-60.0, -84.0), (-60.0, -88.0), (-36.0, -88.0), (-36.0, -84.0), (-60.0, -84.0)]) 34 35 print("shapes[5].points lenth:",len(shapes[5].points)) #('shapes[5].points lenth:', 5) 36 print("shapes[5].shapeType:",shapes[5].shapeType) # shapeType:表示由shapefile规范定义的形状类型的整数。 37 # ('shapes[5].shapeType:', 5) 38 39 fields = sf.fields # 字段列表 40 print("多少字段:%s" % len(fields)) # 多少字段:9 41 print(fields) 42 # [('DeletionFlag', 'C', 1, 0), ['name', 'C', 32, 0], ['column', 'C', 32, 0], ['row', 'C', 32, 0], ['NS', 'C', 32, 0], 43 # ['type', 'C', 32, 0], ['latitude', 'C', 32, 0], ['lat_diff', 'C', 32, 0], ['lon_diff', 'C', 32, 0]] 44 # 字段名:描述此列索引处的数据的名称。 45 # 字段类型:此列索引处的数据类型。类型可以是:字符,数字,长,日期或备忘。 “备忘”类型在GIS中没有意义,而是xbase规范的一部分。 46 # 字段长度:在此列索引处找到的数据的长度。较旧的GIS软件可能会将此长度截短为“字符”字段的8或11个字符。 47 # 小数长度:在“数字”字段中找到的小数位数。 48 49 50 records1 = sf.records() # 通过调用records()方法获取shapefiles记录的列表。 51 print(len(records1)) # 2130 52 print(records1[5]) # ['DS2206', '22', '06', 'S', 'D', '-88--76', '4', '24'] 53 54 55 shapeRecs = sf.shapeRecords() 56 # 同时读取几何和记录 57 # 您希望同时检查记录的几何和属性。 shapeRecord()和shapeRecords()方法让你做到这一点。 58 # 59 # 调用shapeRecords()方法将返回所有形状的几何和属性作为ShapeRecord对象的列表。 60 # 每个ShapeRecord实例都有一个“shape”和“record”属性。 形状属性是一个ShapeRecord对象,在第一部分“阅读几何”中被分割。 61 # 记录属性是如“读取记录”部分中所示的字段值列表。 62 shapeRecs[3].record[1:3] 63 points = shapeRecs[3].shape.points[0:2]
上述主要将的是shape文件的读,下面我举个例子来说明怎么写shape文件:
1:1000000地形图分幅:
1 #! /usr/bin/env python 2 # -*- coding:utf-8 -*- 3 4 from public import create_shapefile 5 import settings 6 7 # 1:1 000 000 地形图分幅 8 class CreateShapefile(object): 9 10 def __init__(self,filename): 11 self.filename = filename 12 self.map_type = settings.MAP_TYPE 13 self.fields = settings.MAP_FIELDS 14 self.w = create_shapefile(self.fields) 15 for item in settings.LAT_LON_RANGE: 16 self.create_poly(item) 17 self.w.save(self.filename) 18 19 def create_poly(self, args): 20 lat_start,lat_final,lat_diff = args[0],args[1],args[2] 21 lon_start,lon_final,lon_diff = args[3],args[4],args[5] 22 latitude = "{0}/{1}".format(str(lat_start),str(lat_final)) # 纬度范围 23 for lat in range(lat_start,lat_final,lat_diff): 24 for lon in range(lon_start,lon_final,lon_diff): 25 26 a = [lon, lat + lat_diff] # 左上角的点坐标 27 b = [lon, lat] # 左下角的点坐标 28 c = [lon + lon_diff, lat] # 右上角的点坐标 29 d = [lon + lon_diff, lat + lat_diff] # 右下角的点坐标 30 31 row, NS = self.row_1_1000000(lat,lat_diff) # 行号和南北半球 32 column = self.column_1_1000000(lat,lon,lon_diff) # 列号 33 name = "D{0}{1}{2}".format(NS, row, column) # 分幅编码 34 35 self.w.poly(parts=[[a, b, c, d, a]]) # 为多边形添加点 36 self.w.record(name,name, row, column, NS, self.map_type, latitude, lat_diff, lon_diff) # 添加字段的值 37 38 39 def row_1_1000000(self,lat,lat_diff): 40 """ 41 功能:计算1:1000000行号和南北半球。 42 输入:lat 纬度 43 lat_diff 纬差 44 返回:row 行号 45 NS 南北半球 46 """ 47 if lat < 0: 48 NS = "S" 49 row = abs(lat) / lat_diff 50 else: 51 NS = "N" 52 row = (lat + lat_diff) / lat_diff 53 if row < 10: 54 row = "0{0}".format(row) 55 return row,NS 56 57 def column_1_1000000(self,lat,lon,lon_diff): 58 """ 59 功能:计算1:1000000列号。 60 输入:lat 纬度 61 lon 经度 62 lon_diff 经差 63 返回:column列号 64 """ 65 if -60 <= lat < 60:x = 31 66 elif 60 <= lat < 76 or -76 <= lat < -60:x = 16 67 else:x = 8.5 68 column = int(float(lon)/lon_diff + x) 69 if column < 10: column = "0{0}".format(column) 70 return column 71 72 if __name__ == "__main__": 73 shape_file = CreateShapefile("shapefile/d_map_1_1000000")
作者:ZingpLiu
出处:http://www.cnblogs.com/zingp/
本文版权归作者和博客园共有,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接。