python与地理空间分析（3）shp文件操作

MeteoAI

发布于 2019-08-19 11:16:37

14.1K1

发布于 2019-08-19 11:16:37

文章被收录于专栏：MeteoAIMeteoAI

引言

在python与地理空间分析（1）与（2）中我们介绍了GIS中常用的数据类型、python在处理地理空间数据时用到的包以及给定经纬度计算空间距离的算法，本期我们主要介绍对地理空间分析中常用到的矢量数据shp文件的处理，在捍卫祖国领土从每一张地图开始我们也提供较为精准的包括南海九段线的中国地图，大家可以自行下载。

Shapefile文件介绍

地理空间分析中有多种矢量数据，大家最常见的是Shapefile（.shp）文件和GeoJSON文件（常用于前端网站发布使用的地理数据格式），GeoJSON文件几乎和python的字典和列表等数据类型一模一样，可以通过python自带的json库直接解析。Shapefile文件是地理信息软件公司Esri在1998年作为一种开放规范发布的矢量数据格式，并逐渐成为GIS数据的一种标准，目前几乎所有的地理空间分析软件都提供对Shapefile文件的支持。Shapefile文件的结构包括多个文件，最重要的文件包括.shp，.shx，.dbf以及.prj文件：

文件扩展	作用	备注
.shp	用于存储要素几何的主文件，其中包括几何图形	必要文件，有的软件只需要shp文件
.shx	形状要素索引文件，适当尺寸的几何元素索引信息可以加快访问速度	必要文件必须和shp文件在一起
.dbf	数据库文件，其中包含几何元素的属性信息	必要文件，可以通过excel打开，查看属性信息
.prj	地图投影信息	shp文件如果需要投影，必备
.sbn	空间bin文件，Shapefile的索引文件	包含一个特征的边框
.sbx	.sbn文件的索引记录文件	常用的空间索引的有序记录索引
.cpg	.dbf的代码文件	为.dbf文件提供国际化支持
.shp.xml	元数据	地理空间元数据.xml的容器

需要注意的是在拷贝shp文件时，需要至少把前4个文件一起拷贝，单独拷贝shp文件无法读取 python中提供了多种处理Shapefile文件的第三方包，例如PyShp，Shapely，Fiona，GeoPandas以及basemap和cartopy包对shp文件的可视化，但大多都是基于OGR库，因此本文主要介绍利用OGR库对shp文件的处理。

GDAL/OGR库

Geospatial Data Abstraction Library （GDAL）是使用C/C++语言编写的用于读写空间数据的一套跨平台开源库。现有的大部分GIS或者遥感平台，不论是商业软件ArcGIS，ENVI还是开源软件GRASS，QGIS，都使用了GDAL作为底层构建库。GDAL库由OGR和GDAL项目合并而来，GDAL主要用于空间栅格数据的读写，OGR主要用于空间要素矢量矢量数据的解析。此外，空间参考及其投影转换使用开源库 PROJ.4进行。OGR提供对矢量数据格式的读写支持，它所支持的文件格式包括：ESRI Shapefiles， S-57， SDTS， PostGIS，Oracle Spatial， Mapinfo mid/mif ， Mapinfo TAB。 GDAL的安装：

window:
step1: 在http://www.gisinternals.com/release.php根据编译器和操作系统位数，选择相应的gdal下载链接，下载GDAL Core和GDAL Bindings两个文件，
step2：系统变量-->Path变量，添加GDAL安装路径
linux（ubuntu）:
step1：sudo add-apt-repository ppa:ubuntugis && sudo apt-get update
setp2：sudo apt-get install gdal-bin
step3：sudo pip install gdal
MacOS:
step1：brew install gdal
step2：sudo pip install GDAL

OGR 处理Shapefile文件

OGR包括如下几部分：

•Geometry：类Geometry (包括OGRGeometry等类)封装了OpenGIS的矢量数据模型，并提供了一些几何操作，WKB(Well Knows Binary)和WKT(Well Known Text)格式之间的相互转换，以及空间参考系统(投影)。•Spatial Reference：类OGRSpatialReference封装了投影和基准面的定义。•Feature：类OGRFeature封装了一个完整feature的定义，一个完整的feature包括一个geometry和geometry的一系列属性。•Feature Definition：类OGRFeatureDefn里面封装了feature的属性，类型、名称及其默认的空间参考系统等。一个OGRFeatureDefn对象通常与一个层(layer)对应。•Layer：类OGRLayer是一个抽象基类，表示数据源类OGRDataSource里面的一层要素(feature)。•Data Source：类OGRDataSource是一个抽象基类，表示含有OGRLayer对象的一个文件或一个数据库。•Drivers：类OGRSFDriver对应于每一个所支持的矢量文件格式。类OGRSFDriver由类OGRSFDriverRegistrar来注册和管理。 OGR读取数据的流程（读取中国行政区划shp文件）：

#导入库
try: 
  from osgeo import ogr
except:
  import ogr
#加载相应数据类型的驱动,相当于初始化一个对象
driver = ogr.GetDriverByName('ESRI Shapefile')
#打开数据
fileName="./中国shp（包含九段线）/Province_9/Province_9.shp"
dataSource = driver.Open(fileName,0) #0是只读，1可写
if dataSource is None:
  print('could not open')
  sys.exit(1)
#获取图层
layer = dataSource.GetLayer(0)
#查看图层的信息
print("图层描述 :{0}".format(layer.GetDescription()))
print("图层范围 :{0}".format(layer.GetExtent()))
print("要素数量 :{0}".format(layer.GetFeatureCount()))
print("元数据描述 :{0}".format(layer.GetMetadata()))
print("空间参考 :{0}".format(layer.GetSpatialRef()))
'''
输出信息：
图层描述 :Province_9
图层范围 :(73.48734510162447, 134.77518963904342, 3.4026791910340695, 53.561602387407575)
要素数量 :4470
元数据描述 :{'DBF_DATE_LAST_UPDATE': '2019-08-05'}
空间参考 :GEOGCS["WGS 84",
  DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
          AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0,
      AUTHORITY["EPSG","8901"]],
  UNIT["degree",0.0174532925199433,
      AUTHORITY["EPSG","9122"]],
  AUTHORITY["EPSG","4326"]]
'''
#获取图层的字段信息
layerDefinition = layer.GetLayerDefn()
for i in range(layerDefinition.GetFieldCount()):
  fieldName =  layerDefinition.GetFieldDefn(i).GetName()
  fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
  fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
  fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
  GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()
  print("字段名：{0}  字段类型：{1}  字段长度：{2}  \
          字段精度：{3}".format(fieldName,fieldTypeCode,
                          fieldType,fieldWidth,GetPrecision))
'''
输出信息：
字段名：AREA  字段类型：2  字段长度：Real      字段精度：13
字段名：PERIMETER  字段类型：2  字段长度：Real 字段精度：13
字段名：SHENG_  字段类型：0  字段长度：Integer 字段精度：9
字段名：SHENG_ID  字段类型：0  字段长度：Integer  字段精度：9
字段名：SHENG  字段类型：0  字段长度：Integer  字段精度：4
字段名：name  字段类型：4  字段长度：String    字段精度：50
'''
#获取要素,并获取要素相应的字段
feat = layer.GetFeature(41)
s_id=feat.GetField('SHENG_ID') # s_id输出为22
geom = feat.GetGeometryRef()
'''
为多边形
POLYGON ((130.241590816154 42.8319851706773,130.241474452665 42.8319288004309,130.241344408453 42.8319760277996,130.241246218978 42.8321642867922,130.241215130657 42.8324318860684,130.241251657616 42.8325858252071,130.241319232586 42.8327018210134,130.24142516911 42.8327647300169,130.241541908539 42.8327337639682,130.241602683179 42.8326719024557,130.241718140689 42.8324253055427,130.241708835884 42.8321698350603,130.241590816154 42.8319851706773))
'''
dataSource.Destroy() #关闭数据源

读完数据，就是新建数据，以某一天的AQI观测数据为例（csv文件），新建shp点文件：

#导入库
import os
import pandas as pd
import osr
try: 
    from osgeo import ogr
except:
    import ogr
#读取AQI数据
AQI=pd.read_csv("AQI.csv")
#加载相应数据类型的驱动,相当于初始化一个对象
driver = ogr.GetDriverByName('ESRI Shapefile')
#创建shp文件
fileName='AQI_POINT.shp'
datasource = driver.CreateDataSource(fileName) #创建 Shape 文件
if os.access(fileName, os.F_OK ): #如文件已存在，则删除
    driver.DeleteDataSource(filename)
#定义空间参考
spatialref = osr.SpatialReference()
spatialref.ImportFromEPSG(4326)
#定义数据类型为POINT
geomtype = ogr.wkbPoint
#创建图层
layer = datasource.CreateLayer("AQI", srs=spatialref, geom_type=geomtype)
#将字段列表写入图层 [station,aqi,pm25,pm10,so2,no2,co,o3,latitude,longitude]
field_station = ogr.FieldDefn("station",ogr.OFTString)
field_station.SetWidth(24)  # 设置长度
layer.CreateField(field_station)
field_aqi = ogr.FieldDefn("aqi", ogr.OFTReal)
layer.CreateField(field_aqi)
layer.CreateField(ogr.FieldDefn("latitude", ogr.OFTReal)) 
layer.CreateField(ogr.FieldDefn("longitude", ogr.OFTReal))
#创建要素
for i in AQI.index:
    # 创建要素
    feature = ogr.Feature(layer.GetLayerDefn())
    # 和设置字段内容进行关联  ,从数据源中写入数据
    print(AQI.loc[i,'aqi'])
    feature.SetField("station", AQI.loc[i,'station'])
    feature.SetField("aqi", float(AQI.loc[i,'aqi']))
    feature.SetField("longitude", AQI.loc[i,'longitude'])
    feature.SetField("latitude", AQI.loc[i,'latitude'])
    # 创建WKT 文本点
    wkt = "POINT(%f %f)" % (
        float(AQI.loc[i,'longitude']), float(AQI.loc[i,'latitude']))
    # 生成实体点
    point = ogr.CreateGeometryFromWkt(wkt)
    # 使用点
    feature.SetGeometry(point)
    # 添加点
    layer.CreateFeature(feature)
    # 关闭 特征
    feature = None
datasource.Destroy() #关闭文件

AQI点shp文件展示

新建其他新的几何形状

#line 线
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
line.SetPoint(0,30,30)   #(10,10) -> (30,30)

#ring 环
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(0,0)
ring.AddPoint(100,0)
ring.AddPoint(100,100)
ring.AddPoint(0,100)
#结束的时候，用CloseRings关闭ring，或者将最后一个点的坐标设定为与第一个点相同
ring.CloseRings()
ring.AddPoint(0,0)

#polygon 多边形 由两个ring构成
outring = ogr.Geometry(ogr.wkbLinearRing)
outring.AddPoint(0,0)
outring.AddPoint(100,0)
outring.AddPoint(100,100)
outring.AddPoint(0,100)
outring.AddPoint(0,0)

inring = ogr.Geometry(ogr.wkbLinearRing)inring = ogr.Geometry(ogr.wkbLinearRing)
inring.AddPoint(25,25)
inring.AddPoint(75,25)
inring.AddPoint(75,75)
inring.AddPoint(25,75)
inring.CloseRings()

polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(outring)
polygon.AddGeometry(inring)

shp文件与其他数据类型的转换

#导入库
try: 
    from osgeo import ogr
except:
    import ogr
#加载相应数据类型的驱动,相当于初始化一个对象
driver = ogr.GetDriverByName('ESRI Shapefile')
#打开数据
fileName="./中国shp（包含九段线）/Province_9/Province_9.shp"
dataSource = driver.Open(fileName,0) #0是只读，1可写
#转换为geojson格式
json="Province_9.json"
ogr.GetDriverByName("GeoJSON").CopyDataSource(dataSource, json)
#转换为google earth的kml格式
kml="Province_9.kml"
ogr.GetDriverByName("KML").CopyDataSource(dataSource, kml)
#转换为csv格式，即把属性表导出和.dbf文件一致
csv="Province_9.csv"
ogr.GetDriverByName("CSV").CopyDataSource(dataSource, csv)