shapefile文件_入门-Python-读取Shapefile

前言

补一篇Shapefile的开源读取，Shapefile格式作为GIS的标准格式，几乎是各大GIS平台必须支持的一个交换格式，所以在各个平台上进行SHP的读写都是没问题的。但是如何在开发过程中读写SHP？是否需要在各大GIS平台中通过二次开发来读取？答案是不需要。本文介绍Python开源的读取方式，也顺带引出GDAL这个GIS最为知名的开源类库，为各位GIS开发初学者指一条阳关大道。

编程难就难在选择正确的大道，为什么？条条大路通罗马，编程去完成需求一样是八仙过海，但有些方法仅限于实现，有些方法没有扩展性，要知道实现与实用虽然是一字之差，但程序员在重构的道路上为此付出良多，如果一开始就能选择走正确的道路，那么以后就会走的轻松些。当然，如何选择，这就要靠经验的积累。

GDAL

1. 地理数据处理软件包GDAL简介​www.osgeo.cn

本段为摘抄：GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。该项目由Frank Warmerdam教授于1998年发起。 它利用抽象数据模型来表达所支持的各种文件格式。 它还有一系列命令行工具来进行数据转换和处理。 OGR(OpenGIS Simple Features Reference Implementation)是GDAL项目的一个子项目， 提供对矢量数据的支持。 一般把这两个库合称为GDAL/OGR，或者简称为GDAL。

简单理解：原GDAL是栅格数据的API，OGR是矢量数据的API，现GDAL=原GDAL+OGR。

GDAL类库是GIS专业做开发的同学，最应该学习和了解的类库，原因是该类库被各大GIS平台引用，所以无论将来开发在ArcGIS体系下做，还是在其它开源体系下，都会用到。第二个重要原因，开发无论是做的浅还是做的深，都可以学习，浅的可以只学习API进行调用，深的可以学习源代码了解shp的数据结构，了解GIS如何进行空间分析等等。
如果倒退十几年回到研究生学习GIS开发的岁月，肯定会优先学习它，这样的GIS开发者会有一条较为平滑的学习曲线，而不是一只懵懂的小鹿只会乱撞，我的开发学习路线是：MapGIS二次开发——ArcGIS全体系（AO+AE+ArcIMS+Server+JS API）——Geoserver——Leaflet，完全是凭兴趣爱好再硬撑，要不是ArcGIS体系永葆青春，估计早就倒在这崎岖的山坑之中。

请GIS专业的研究生，以后想从事开发工作的，学习该类库。学习语言有很多，可以是C++，Python甚至Javascript（NodeJS）。

GDAL的Python安装

该类库不能直接通过pip install。步骤最方便的可以按：

1.下载GDAL-Python-Windows安装包，链接地址，选择对应python版本和windows32或64版本，如python37和64位操作系统选择GDAL‑3.0.1‑cp37‑cp37m‑win_amd64.whl。

2.下载后，pip3 install GDAL‑3.0.1‑cp37‑cp37m‑win_amd64.whl，注意cmd执行路径，切换到GDAL‑3.0.1‑cp37‑cp37m‑win_amd64.whl所在目录。

接下来，可上一篇中配置的VS Code IDE中进行编码和类库调用

源码

import sys
from osgeo import gdal
from osgeo import ogr
#shapefile数据源
driver = ogr.GetDriverByName('ESRI Shapefile')
filename = 'D:CodePythonDatariver.shp'
dataSource = driver.Open(filename,0)
if dataSource is None:
    print('could not open')
    sys.exit(1)
#获取图层    
layer = dataSource.GetLayer(0)    
n = layer.GetFeatureCount()
print('feature count:', n)
#字段定义
layerDefinition = layer.GetLayerDefn()
#遍历要素
feat = layer.GetNextFeature()  
while feat:
    for i in range(layerDefinition.GetFieldCount()):
        fieldName =  layerDefinition.GetFieldDefn(i).GetName()
        value = feat.GetField(fieldName)
        print(fieldName + " : " + str(value))
    feat = layer.GetNextFeature()
#游标复位    
layer.ResetReading()  
print('done!')

源码很简单直观，按注释基本可了解整个读取过程，当然，如果要了解更多GDAL和OGR的API可以参考：

Welcome to the Python GDAL/OGR Cookbook!​pcjericks.github.io

重点来了，看了上面几个API，如果了解ArcGIS AO或AE开发的同学，是不是有似曾相识的感觉，没错！AO和AE的一些重要的数据类和结构的定义是参考GDAL得来的，IFeatureClass，IFeature，IFieldDef，IField等等这些ArcGIS接口定义几乎与GDAL中定义一样，非常容易理解。（说到这里，吐槽下Leaflet的类定义，Leaflet的Layer作为图层，它的含义竟然不是GIS的图层，而是类似于PhotoShop的图层，这让GIS科班出身的人有点难受。）

题外话

1.GDAL类库也有坐标转换的接口API，后期有空也可另文总结。

2.NodeJS对应node-gdal，直接npm install node-gdal。其它文章中，有关从零实现GIS Server，通过NodeJS来读取SHP的文章中已有提到，此处略过。

3.关于GIS开发的学习，一直在路上。问我，我也才刚起步。