使用pyshp读写Shapefile
Python Shapefile库（pyshp）为Esri Shapefile格式提供读写支持。 Shapefile格式是Esri创建的较为流行的地理信息系统矢量数据格式。 更多信息，请阅读 http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf 的“ESRI Shapefile技术说明 - 1998年7月”。 Esri文档主要讲的是shp和shx文件格式。 但dbf格式其实很常用。 这种格式作为“XBase文件格式描述”常作为网络上记录，dbf格式是在20世纪60年代创建的，是较为简单的基于文件的数据库格式。
shapefile是GIS中非常重要的一种数据类型，在ArcGIS中被称为要素类(Feature Classes)，主要包括点(point)、线(polyline)和多边形(polygon)。Python脚本是ArcGIS官方推荐的脚本语言，Python脚本调用ArcGIS中各种工具和函数并且还能批量完成所需操作。本文介绍的Python Shapefile Library是一个Python库，在Python脚本中用于读取对ArcGIS中的Shapefile文件(.shp，.shx，.dbf等格式)。
Esri和XBase文件格式在设计和存储效率方面非常适宜，这是shapefile格式流行的原因之一，尽管当今存储和交换GIS数据的方法众多。
Pyshp与Python 2.4-3.x兼容。本文档提供的是pyshp读写shapefile的示例。 然而，更多的例子在GitHub pyshp wiki中，博客 http://GeospatialPython.com ，还可以通过 http://gis.stackexchange.com 搜索pyshp。
目前，示例中引用的样本人口普查块组shapefile可在GitHub项目网站 https://github.com/GeospatialPython/pyshp 上获得。 这些例子是开箱即用，你也可以轻松地运行它，只需要对你自己的shapefile小小的修改即可。
提示：如果您是GIS新手，您应该阅读有关地图投影的相关信息。 请访问：https：//github.com/GeospatialPython/pyshp/wiki/Map-Projections
我真诚地希望这个简单地读取并写入数据，真诚地您日后可关注一些地理空间项目的具有挑战与有趣的部分。
Python Shapefile Library的下载与安装
Python Shapefile Library下载地址：https://code.google.com/p/pyshp/ ，目前已托管到 GitHub 中。
Python Shapefile Library是一个纯 Python 库，使用时无需安装，只需在Python程序中导入该模块文件即可。
import shapefile
在做任何操作之前，您必须导入库。
以下示例将使用美国人口普查局Blockgroups数据集创建的shapefile，（美国人口普查局，位于加利福尼亚州旧金山附近），可在pyshp GitHub站点的github存储库中找到。具体导入方法请参考 Python 教程中模块的导入章节。
读取Shapefile
Python Shapefile Library提供了Reader类功能，可创建Reader类的对象，完成shapefile文件的读操作。
使用的是Python Shapefile Library读取shapefile文件的"几何数据"(Geometry)和"属性数据"(Attribute Record)
"几何数据"一般由多个几何对象组成，比如一个"点文件"，每个点就是一个对象；对于一个多边形文件，每个对象可能包含多个多边形，每个多边形又称为"块(parts)"，每个 "块"由多个点组成。
每个几何对象包含4个属性：数据类型(shapeType)，代表该"几何数据"对象的数据类型(点，shapeType=1，线，shapeType=3，多边形，shapeType=5)；数据范围(bbox)，只针对多点数据，代表该"几何数据"对象的边界范围；数据块(parts)，只针对线或者多边形，代表该"几何数据"对象各个块的第一个点的索引；点集(points)，代表该"几何数据"对象的所有点坐标。
"属性数据"即每个"几何数据"对象在属性表中的对应项。
若要读取shapefile，需创建一个新的“Reader”对象，命名为现有shapefile的名称。 shapefile格式可以写成以下三种文件格式：
import shapefile
sf = shapefile.Reader("/gdata/GSHHS_c")
您可以指定shapefile的基本文件名或任何shapefile组件的完整文件名。
或者
sf = shapefile.Reader("/gdata/GSHHS_c.shp")
或者
sf = shapefile.Reader("/gdata/GSHHS_c.dfb")
或任何一种均是shapefile的格式。 库对文件扩展名并无严格限定。
myshp = open("/gdata/GSHHS_c.shp", "rb")
mydbf = open("/gdata/GSHHS_c.dbf", "rb")
r = shapefile.Reader(shp=myshp, dbf=mydbf)
注意在述示例中，shx文件从不使用。 shx文件是shp文件中的可变长度记录的非常简单的固定记录索引。该文件的可读性是可选的。 如果可用pyshp将使用shx文件访问形状记录一点更快，但会做没有它只是没有它。
shapes = sf.shapes()
通过调用shapes（）方法获得shapefile的几何的列表。Shapes是一个列表(相当于一维数组)，存放的是该文件中所有的"几何数据"对象
len(shapes)
shapes方法返回的是描述每个形状记录的几何形状的Shape对象列表。
len(list(sf.iterShapes()))
您可以使用iterShapes（）方法遍历shapefile的几何数据。
通过shapeType，bbox，points，parts方法返回每个"几何数据"对象的属性信息：
每个形状记录均包含以下属性：
for name in dir(shapes[3]):
if not name.startswith('__'):
name
shapes[3].shapeType
结果返回第1个对象的数据类型属性(或者Shape.shapeType)。
shapeType：表示由shapefile规范定义的形状类型整数。
# Get the bounding box of the 4th shape.
# Round coordinates to 3 decimal places
bbox = shapes[3].bbox
['%.3f' % coord for coord in bbox]
以上结果返回第3个对象的数据范围(左下角的x，y坐标和右上角的x，y坐标)
bbox：如果形状类型包含多个点，则此元组描述左下（x，y）坐标和右上角（x，y）坐标。 如果shapeType是Null（shapeType == 0），结果就会出现AttributeError。
shapes[3].parts
parts：parts将点集合组合为形状。 如果形状具有多个部分，则显示的是每个部分的第一点索引。 如果只有一个部分，则返回包含0的列表。
len(shapes[3].points)
此结果返回第4个对象的所有点坐标
# Get the 8th point of the fourth shape
# Truncate coordinates to 3 decimal places
shape = shapes[3].points[7]
['%.3f' % coord for coord in shape]
- points：points属性包含一个元组列表，其中包含形状中每个点（x，y）坐标。
s = sf.shape(7)
Shape是第8个"几何数据"对象
# Read the bbox of the 8th shape to verify
# Round coordinates to 3 decimal places
['%.3f' % coord for coord in s.bbox]
如果通过调用其索引来读取单个形状，请使用shape（）方法。 索引形状的计数从0开始的。所以读取第8个形状记录时，你使用的是它的索引是7。
"属性数据"的读取方法
"属性数据"通过Reader类的records( )和record( )方法来读取，其区别与使用方法同shapes( )和shape( )方法。 "属性数据"的fields可通过Reader类的fields方法获取，其返回值为列表属性表中每个字段的名称、数据类型、数据长度等。
shapefile记录包含几何体集合体每个形状的属性，并存储在dbf文件中，隐含着shp几何文件和dbf属性文件的形状和对应记录顺序。
在读取shapefile时，可以使用shapefile的字段名称。你可以将shapefile的“字段”属性称为一个Python列表。每个字段的Python列表均包含以下信息：
- 字段名：描述此列索引处的数据名称。
- 字段类型：此列索引处的数据类型。类型可以是：字符，数字，长，日期或备忘录。备忘录在GIS中没有意义，而是作为xbase规范的一部分。
- 字段长度：此列索引处的数据长度。较旧的GIS软件可能会将此长度截短为“字符”字段的8或11个字符。
- 小数长度：“数字”字段中的小数位数。
fields = sf.fields
fields
若查看上面的Reader对象（sf）字段，需要调用“fields”属性。
records = sf.records()
len(records)
您可以通过调用records（）方法来获取shapefile的记录列表。
len(list(sf.iterRecords()))
与geometry方法类似，可以使用iterRecords（）方法对dbf记录进行迭代。
每个记录是与字段列表中的每个字段相对应的属性的列表。
例如，在blockgroups shapefile的记录中，第2和第3个字段是块组id和该旧金山块组的1990年人口计数：
records[3][1:3]
rec = sf.record(3)
rec[1:3]
要读取单个记录，使用的是记录的索引调用record（）方法：
shapeRecs = sf.shapeRecords()
调用shapeRecords（）方法将返回所有形状的几何和属性值，以此作为ShapeRecord对象的列表。 每个ShapeRecord实例都有一个“shape”和“record”属性。 “shape”属性是一个ShapeRecord对象。 记录属性的是一个字段值列表：
ShapeRecords = sf.shapeRecords( )
ShapeRecords[0].shape.shapeType
返回第1个对象的"几何数据"的数据类型属性：
shapeRecs[3].record[1:3]
读取同一条记录的前两点：
ShapeRecords[0].record[1:3]
返回第1个对象的"属性数据"的第2和第3个属性值
points = shapeRecs[3].shape.points[0:2]
len(points)
shapeRec = sf.shapeRecord(3)
shapeRecord（）方法读取的是指定索引处的单个形状或记录。 要从blockgroups shapfile获取第四个形状记录，使用第三个索引.
shapeRec.record[1:3]
points = shapeRec.shape.points[0:2]
len(points)
块组密钥和人口计数.
shapeRecs = sf.iterShapeRecords()
for shapeRec in shapeRecs:
# do something here
pass
最后，用iterShapeRecords（）方法遍历所以文件.
写shapefile
Python Shapefile Library提供了Writer类。
shapefile文件的创建
shapefile文件的创建分为2步：
1、创建"几何数据"，通过Writer类的point(x,y,z,m)方法创建点数据，通过poly(parts = [ [ [1,5], [5,5], [5,1], [3,3], [1,1] ] ])方法创建线(parts有2个点)或多边形数据(parts > 2个点)。
2、创建"属性数据"，首先通过field()方法创建属性表字段，然后通过record()方法为每个几何数据添加相应的属性值。
PyShp在编写shapefile时应尽可能灵活些，同时还需要有一定程度的自动验证功能，以确保不会写入无效文件。
PyShp可以只写一个组件文件，如shp或dbf文件，而不写其他。 因此，除了作为一个完整的shapefile库，它也可以用作基本的dbf（xbase）库。 Dbf文件是一种通用的数据库格式，通常用作独立的简单数据库格式。 甚至shp文件偶尔还会有用作一个独立的格式的情况。 一些基于Web的GIS系统使用用户上传的shp文件来指定感兴趣的区域。 许多精确的农业化学领域喷雾器也使用shp格式作为喷雾器系统的控制文件（通常结合自定义数据库文件格式）。
要创建shapefile，可以使用Writer类方法添加几何与属性，直到准备好保存文件。
import shapefile
w = shapefile.Writer()
以上方法是使用创建shapefile Writer类的实例：
设置形状类型
形状类型是shapefile中包含的几何类型。 所有形状必须与形状类型设置相匹配。
形状类型由shapefile规范定义的0和31之间的数字表示。 要注意，编号系统有几个尚未使用的保留编号，因此现有形状类型的编号并不是按顺序排列的。
有三种方法来设置形状类型：
- 在创建类实例时设置它。
- 通过为现有类实例分配值来设置它。
- 通过保存shapefile将其自动设置为第一个形状的类型来设置它。
通过创建Writer类的对象(如下面的sf)进行shapefile文件的write操作.
文件类型的确定
写shapefile文件时，首先要确定shapeType，可以通过以下两种方法确定:
1、在创建Writer类对象时直接确定shapeType，如上所示
2、通过为Writer类对象的shapeType属性赋值，如sf.shapeType = 1
w = shapefile.Writer(shapeType=1)
w.shapeType
要在创建Writer时手动设置Writer对象的形状类型：
w.shapeType = 3
w.shapeType
或者您可以在创建Writer后进行设置.
几何和记录平衡
因为每个形状都必须有其相应的记录，所以关键的是记录形状的数量以创建有效的shapefile。 为了防止意外出错，PSL具有“自动平衡”功能，以确保添加形状或记录时，方程式的两边对齐。 此功能在默认情况下不启用。
"几何数据"与"属性数据"的自动平衡
shapefile文件要求"几何数据"与"属性数据"一一对应，如果有"几何数据"而没有相应的属性值存在时，那么在打开ArcGIS软件时所创建的shapefile文件就会出错。为了避免这种情况的发生，可以设置 sf.autoBalance = 1，以确保每创建一个"几何数据"，该库会自动创建一个属性值(空的属性值)来对应。autoBalance默认为0。
w.autoBalance = 1
若要激活它，可以将属性autoBalance设置为1（True）.
您还可以在添加形状或记录时手动调用balance（），以确保另一方是最新的。当使用平衡功能时，在几何侧创建空形状，或者在属性数据处创建具有值“NULL”的每个字段的记录。因此，平衡选项保障构建shapefile的灵活性。
如果没有自动平衡，您可以随时添加几何或记录。 您可以创建所有形状，创建所有记录，反之亦然。 您可以在每次创建形状或记录后使用balance方法，并及时更新。 如果不使用平衡方法并忘记手动平衡几何和属性的情况下，shapefile将被大多数shapefile软件视为已损坏。
使用自动平衡，您可以添加形状或几何体，并根据需要在任一侧更新空白条目。 即使您忘记更新条目，shapefile也仍然有效，并且会被大多数shapefile软件视为正常，作正确处理。
w = shapefile.Writer()
w.point(122, 37) # No elevation or measure values
w.shapes()[0].points
w.point(118, 36, 4, 8)
w.shapes()[1].points
使用“点”方法添加点形状。 点由x，y和可选的z（高程）和m（测量）值指定。
w = shapefile.Writer()
w.poly(shapeType=3, parts=[[[122,37,4,9], [117,36,3,4]], [[115,32,8,8],[118,20,6,4], [113,24]]])
poly可以是多边形也可以是线。 Shapefile多边形至少得有4个点，最后一个点必须与第一个点相同。 而PyShp是自动强制关闭多边形。 线必须至少有两个点。 由于这两种形状类型之间有太多的相似性，因此它们统一使用“poly”这一单一方法创建。
w = shapefile.Writer()
w.null()
因为Null形状类型（形状类型为0）并无几何形状，所以在没有任何参数的情况下可以调用“null”方法。 这种类型的shapefile很少使用，但它是有效的。
assert w.shapes()[0].shapeType == shapefile.NULL
上述语句是writer对象的形状列表是一个null形状。
w = shapefile.Writer(shapefile.POINT)
w.point(1,1)
w.point(3,1)
w.point(4,3)
w.point(2,2)
w.field('FIRST_FLD')
w.field('SECOND_FLD','C','40')
w.record('First','Point')
w.record('Second','Point')
w.record('Third','Point')
w.record('Fourth','Point')
w.save('shapefiles/test/point')
w = shapefile.Writer(shapefile.POLYGON)
w.poly(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record('First','Polygon')
w.save('shapefiles/test/polygon')
w = shapefile.Writer(shapefile.POLYLINE)
w.line(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
w.poly(parts=[[[1,3],[5,3]]], shapeType=shapefile.POLYLINE)
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record('First','Line')
w.record('Second','Line')
w.save('shapefiles/test/line')
您还可以使用关键字参数添加属性，其中关键字是字段名称。
w = shapefile.Writer(shapefile.POLYLINE)
w.line(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record(FIRST_FLD='First', SECOND_FLD='Line')
w.save('shapefiles/test/line')
targetName = w.save()
assert("shapefile_" in targetName)
如果不指定任何文件名（即save（）），那么将生成一个唯一的文件名，其前缀为“shapefile_”，随后是用于所有三个文件的随机字符。 唯一的文件名将作为字符串返回。
try:
from StringIO import StringIO
except ImportError:
from io import BytesIO as StringIO
shp = StringIO()
shx = StringIO()
dbf = StringIO()
w.saveShp(shp)
w.saveShx(shx)
w.saveDbf(dbf)
# Normally you would call the "StringIO.getvalue()" method on these objects.
shp = shx = dbf = None
Python geo_interface
Python geo_interface提供了地理空间与Python库之间的数据交换接口。 接口以GeoJSON返回数据，这使得其他库与工具之间（包括Shapely，Fiona和PostGIS）具有良好的兼容性。 有关geo_interface协议的更多信息，请访问：https：//gist.giithub.com/sgillies/2217756。 有关GeoJSON的更多信息，请访问http://geojson.org。
s = sf.shape(0)
s.__geo_interface__["type"]
import helper; helper.info()