目录

上一个主题

9. Python下面其他开源GIS库的使用

下一个主题

9.2. GeoJSON

关注公众号


常见问题

  1. Windows下的安装说明
  2. Jupyter免费在线实验环境
  3. 勘误与补充


>>> from helper import info; info()
待更新

9.1. 使用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

我真诚地希望这个简单地读取并写入数据,真诚地您日后可关注一些地理空间项目的具有挑战与有趣的部分。

9.1.1. 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 教程中模块的导入章节。

9.1.2. 读取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的格式。 库对文件扩展名并无严格限定。

从文件类对象读取shapefile

在任意Python文件对象中您还可以用关键字参数加载shapefile,可以指定三个文件中的任何一个。 此功能非常强大,您可以从url,zip文件,序列化对象或某些情况下的数据库中加载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文件访问形状记录一点更快,但会做没有它只是没有它。

“几何数据”的读取方法

“几何数据”通过Reader类的shapes( )和shape( )方法读取,二者的区别在于:shapes()方法不需要指定参数,其返回值是一个列表,包含该文件中所有的“几何数据”对象,而shape( )方法需要指定参数,返回的是所需要的“几何数据”对象。

shapefile的几何形状是表示物理位置的隐含弧形的点或形状的集合。 所有类型的shapefile只存储点。

你可以调用shapes()方法获得shapefile的几何列表。

>>> shapes = sf.shapes()

通过调用shapes()方法获得shapefile的几何的列表。Shapes是一个列表(相当于一维数组),存放的是该文件中所有的“几何数据”对象

>>> len(shapes)
1765

shapes方法返回的是描述每个形状记录的几何形状的Shape对象列表。

>>> len(list(sf.iterShapes()))
1765

您可以使用iterShapes()方法遍历shapefile的几何数据。

通过shapeType,bbox,points,parts方法返回每个“几何数据”对象的属性信息:

每个形状记录均包含以下属性:

>>> for name in dir(shapes[3]):
>>>     if not name.startswith('__'):
>>>         name
>>> shapes[3].shapeType
5

结果返回第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]
['-168.105', '7.206', '-55.628', '71.994']

以上结果返回第3个对象的数据范围(左下角的x,y坐标和右上角的x,y坐标)

bbox:如果形状类型包含多个点,则此元组描述左下(x,y)坐标和右上角(x,y)坐标。 如果shapeType是Null(shapeType == 0),结果就会出现AttributeError。

>>> shapes[3].parts
[0, 667, 671, 675, 679, 683, 688, 695, 715, 719, 724, 728, 732, 736, 740, 744, 748, 753, 758, 762, 766, 780, 784, 789, 793, 797, 802, 806, 810, 814, 818, 822, 826, 830, 835, 839, 843, 847, 852, 856, 861, 865, 870, 875, 879, 884, 888, 892, 896, 900, 904, 908, 912, 918, 923, 927, 931, 935, 939, 947, 954, 958, 962, 966, 970, 974, 978, 982, 988, 992, 996, 1001, 1005, 1014, 1022, 1026, 1030, 1034, 1038, 1042, 1046, 1051, 1055, 1059, 1063, 1067, 1071, 1075, 1079, 1083, 1087, 1091, 1095, 1099, 1104, 1108, 1112, 1127, 1132, 1136, 1140, 1144, 1148, 1152, 1156, 1160, 1164, 1168, 1172, 1176, 1180, 1184, 1189, 1194, 1198, 1203, 1207, 1211, 1216, 1220, 1229, 1233, 1282, 1286, 1293, 1298, 1302, 1306, 1310, 1315, 1319, 1323, 1327, 1336, 1340, 1344, 1348, 1352, 1356, 1360, 1364, 1368, 1372, 1383, 1387, 1391, 1395, 1399, 1403, 1408, 1412, 1416, 1420, 1424, 1429, 1442, 1446, 1450, 1454, 1458, 1463, 1467, 1471, 1475, 1479, 1483, 1487, 1491, 1495, 1499, 1503, 1507, 1513, 1517, 1521, 1525, 1529, 1533, 1537, 1541, 1548, 1552, 1556, 1560, 1564, 1569, 1573, 1577, 1581, 1585, 1589, 1599, 1603, 1607, 1611, 1615, 1619, 1625, 1629, 1635, 1643, 1647, 1651, 1655, 1659, 1667, 1671, 1679, 1683, 1687, 1691, 1695, 1699, 1703, 1707, 1711, 1715, 1719, 1726, 1730, 1734, 1740, 1744, 1748, 1752, 1756, 1760, 1764, 1768, 1773, 1777, 1781, 1785, 1789, 1793, 1797, 1801, 1805, 1809, 1813, 1817, 1821, 1825, 1829, 1833, 1837, 1843, 1847, 1851, 1855, 1859, 1863, 1868, 1872, 1876, 1880, 1886, 1892, 1896, 1900, 1904, 1908, 1912, 1917, 1921, 1925, 1929, 1933, 1937, 1941, 1945, 1949, 1953, 1957, 1961, 1965, 1969, 1973, 1977, 1982, 1986, 1990, 1994, 1998, 2002, 2006, 2010, 2014, 2018, 2037, 2041, 2045, 2050, 2054, 2059, 2063, 2067, 2071, 2075, 2079, 2083, 2089, 2093, 2097, 2101, 2105, 2110, 2114, 2118, 2122, 2126, 2130, 2139, 2144, 2148, 2152, 2156, 2160, 2164, 2168, 2172, 2187, 2191, 2195, 2199, 2203, 2207, 2211, 2215, 2219, 2223, 2227, 2231, 2235, 2239, 2244, 2248, 2252, 2257, 2262, 2266, 2270, 2274, 2278, 2282, 2286, 2290, 2295, 2299, 2303, 2307, 2311, 2315, 2319, 2323, 2327, 2331, 2335, 2339, 2343, 2347, 2351, 2355, 2359, 2363, 2367, 2371, 2376, 2380, 2384, 2388, 2392, 2400, 2404, 2408, 2412, 2417, 2421, 2425, 2430, 2435, 2439, 2443, 2447, 2451, 2455, 2459, 2463, 2467, 2471, 2475, 2479, 2483, 2487, 2491, 2496, 2500, 2506, 2510, 2514, 2519, 2523, 2528, 2533, 2537, 2541, 2546, 2550, 2554, 2559, 2563, 2567, 2571, 2575, 2580, 2584, 2588, 2594, 2612, 2616, 2621, 2625, 2629, 2633, 2637, 2641, 2645, 2649, 2653, 2657, 2661, 2665, 2669, 2673, 2677, 2682, 2686, 2691, 2696, 2700, 2704, 2708, 2712, 2716, 2720, 2724, 2728, 2733, 2737, 2741, 2747, 2751, 2755, 2759, 2763, 2767, 2771, 2775, 2779, 2783, 2787, 2791, 2796, 2800, 2804, 2808, 2812, 2816, 2820, 2824, 2828, 2832, 2836]

parts:parts将点集合组合为形状。 如果形状具有多个部分,则显示的是每个部分的第一点索引。 如果只有一个部分,则返回包含0的列表。

>>> len(shapes[3].points)
2840

此结果返回第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]
['-89.331', '69.245']
  • points:points属性包含一个元组列表,其中包含形状中每个点(x,y)坐标。

通过 shape() 读取

>>> 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]
['130.930', '-10.705', '150.877', '-0.339']

如果通过调用其索引来读取单个形状,请使用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
[('DeletionFlag', 'C', 1, 0),
 ['id', 'C', 80, 0],
 ['level', 'N', 10, 0],
 ['source', 'C', 80, 0],
 ['parent_id', 'N', 10, 0],
 ['sibling_id', 'N', 10, 0],
 ['area', 'F', 19, 11],
 ['Shape_Leng', 'F', 19, 11],
 ['Shape_Area', 'F', 19, 11]]

若查看上面的Reader对象(sf)字段,需要调用“fields”属性。

>>> records = sf.records()
>>>
>>> len(records)
1765

您可以通过调用records()方法来获取shapefile的记录列表。

>>> len(list(sf.iterRecords()))
1765

与geometry方法类似,可以使用iterRecords()方法对dbf记录进行迭代。

每个记录是与字段列表中的每个字段相对应的属性的列表。

例如,在blockgroups shapefile的记录中,第2和第3个字段是块组id和该旧金山块组的1990年人口计数:

>>> records[3][1:3]
[1, 'WVS']
>>> rec = sf.record(3)
>>>
>>> rec[1:3]
[1, 'WVS']

要读取单个记录,使用的是记录的索引调用record()方法:

“几何数据”和“属性数据”同时读取的方法

通过Reader类的shapeRecords( )和shapeRecord( )方法可以同时读取shapefile文件的“几何数据”和“属性数据”。

>>> shapeRecs = sf.shapeRecords()

调用shapeRecords()方法将返回所有形状的几何和属性值,以此作为ShapeRecord对象的列表。 每个ShapeRecord实例都有一个“shape”和“record”属性。 “shape”属性是一个ShapeRecord对象。 记录属性的是一个字段值列表:

>>> ShapeRecords = sf.shapeRecords( )
>>>
>>> ShapeRecords[0].shape.shapeType
5

返回第1个对象的“几何数据”的数据类型属性:

>>> shapeRecs[3].record[1:3]
[1, 'WVS']

读取同一条记录的前两点:

>>> ShapeRecords[0].record[1:3]
[1, 'WVS']

返回第1个对象的“属性数据”的第2和第3个属性值

>>> points = shapeRecs[3].shape.points[0:2]
>>>
>>> len(points)
2
>>> shapeRec = sf.shapeRecord(3)

shapeRecord()方法读取的是指定索引处的单个形状或记录。 要从blockgroups shapfile获取第四个形状记录,使用第三个索引.

>>> shapeRec.record[1:3]
[1, 'WVS']
>>> points = shapeRec.shape.points[0:2]
>>>
>>> len(points)
2

块组密钥和人口计数.

>>> shapeRecs = sf.iterShapeRecords()
>>> for shapeRec in shapeRecs:
>>>       # do something here
>>>          pass

最后,用iterShapeRecords()方法遍历所以文件.

9.1.3. 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"]
'Polygon'