Python读写修改Shapefile

小龙狗

已于 2022-03-20 20:26:42 修改

阅读量1w

点赞数 2

分类专栏： GIS Python编程文章标签： python gdal shapefile

于 2021-04-19 11:26:22 首次发布

本文链接：https://blog.youkuaiyun.com/ShyLoneGirl/article/details/115772352

版权

Python编程同时被 2 个专栏收录

42 篇文章

订阅专栏

GIS

6 篇文章

订阅专栏

说明

当我们需要借助编程手段对 Shapefile 进行操作时，会想到用 ArcGIS Engine、ArcGIS Runtime SDK或开源的 GDAL。在没有一定基础的情况下，这些弄起来属实挺麻烦(配环境、写程序等)。当我们做些简单的 Shapefile 数据处理，有没有好办法呢？这里非常推荐使用 Python 语言，介绍三个常用的开发库：

ArcPy
什么是快… ArcPy，下面是文档地址。
https://desktop.arcgis.com/zh-cn/arcmap/10.5/analyze/arcpy/what-is-arcpy-.htm
简言之就是 ArcGIS 的配套工具，好用，但是有限定条件：使用它的前提是安装好 ArcGIS 的环境，另外无论是运行还是打包后执行，都将永久依赖 ArcGIS 。
pyshp
由纯 Python 打造的读、写 Shapefile 的轻量级开源工具，Github 地址。
https://github.com/GeospatialPython/pyshp
如果你只做 Shapefile 读写的话，用它是非常合适的。
GDAL for Python
嗷它还是 GDAL，但是显然 Python 版本更适合于小型项目(可能叫工具更合适)，它比 C++版本的好配太多，也更省去源码编译的过程。功能强大且开源，选它选它。

1. GDAL安装

Python环境就不说了，先到这里下载 GDAL‑<gdal版本>‑cp<python版本>‑<平台架构>.whl 安装包，
https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
下载完成后安装

pip install xxx.whl

安装后试一下能否 import 进来包。

D:\test>python
Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32

>>> from osgeo import gdal
>>>

没有问题。

2. GDAL读取shp

#!/usr/bin/python
# -*- coding: UTF-8 -*-
from osgeo import gdal
from osgeo import ogr
from osgeo import osr

def ReadShp(_filename):
    # 支持中文路径
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    # 支持中文编码
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 注册所有的驱动
    ogr.RegisterAll()
    # 打开数据
    ds = ogr.Open(_filename, 0)
    if ds == None:
        return ("打开文件失败！")
    # 获取数据源中的图层个数，shp数据图层只有一个，gdb、dxf会有多个
    # iLayerCount = ds.GetLayerCount()
    # 获取第一个图层
    oLayer = ds.GetLayerByIndex(0)
    if oLayer == None:
        return ("获取图层失败！")
    # 对图层进行初始化
    oLayer.ResetReading()
    # 输出图层中的要素个数
    num = oLayer.GetFeatureCount(0)
    result_list = []
    # 获取要素
    for i in range(0, num):
        ofeature = oLayer.GetFeature(i)
        # shp 文件的 NAME 字段
        fieldName = ofeature.GetFieldAsString('NAME')
        # geom = str(ofeature.GetGeometryRef())
        result_list.append(fieldName)
    ds.Destroy()
    del ds
    return result_list
    
if __name__ == '__main__':
    
    shpFile = "TM_WORLD_BORDERS-0.3.shp"
    nameList = ReadShp(shpFile)
    for name in nameList:
        print (name)

3. GDAL创建shp

#!/usr/bin/python
# -*- coding: UTF-8 -*-

from osgeo import gdal
from osgeo import ogr
from osgeo import osr

def NewShp(_strDriverName):
    # 支持中文路径
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    # 属性表字段支持中文
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 注册驱动
    ogr.RegisterAll()
    # 创建shp数据
    strDriverName = "ESRI Shapefile"
    oDriver = ogr.GetDriverByName(strDriverName)
    if oDriver == None:
        return ("驱动不可用："+strDriverName)
    # 创建数据源
    oDS = oDriver.CreateDataSource(_strDriverName)
    if oDS == None:
        return ("创建文件失败"+_strDriverName)
    # 创建一个多边形图层，指定坐标系为WGS84
    papszLCO = []
    geosrs = osr.SpatialReference()
    geosrs.SetWellKnownGeogCS("WGS84")
    # 线：ogr_type = ogr.wkbLineString
    # 点：ogr_type = ogr.wkbPoint
    ogr_type = ogr.wkbPolygon
    # 面的类型为Polygon，线的类型为Polyline，点的类型为Point
    oLayer = oDS.CreateLayer("Polygon", geosrs, ogr_type, papszLCO)
    if oLayer == None:
        return ("图层创建失败！")
    # 创建属性表
    # 创建id字段
    oId = ogr.FieldDefn("id", ogr.OFTInteger)
    oLayer.CreateField(oId, 1)
    # 创建name字段
    oName = ogr.FieldDefn("name", ogr.OFTString)
    oLayer.CreateField(oName, 1)
    oDefn = oLayer.GetLayerDefn()
    # 创建要素
    # 数据集
    # wkt_geom id name
    features = ['test0;POLYGON((-1.58 0.53, -0.79 0.55, -0.79 -0.23, -1.57 -0.25, -1.58 0.53))', 'test1;POLYGON((-1.58 0.53, -0.79 0.55, -0.79 -0.23, -1.57 -0.25, -1.58 0.53))']
    for index, f in enumerate(features):
        oFeaturePolygon = ogr.Feature(oDefn)
        oFeaturePolygon.SetField("id",index)
        oFeaturePolygon.SetField("name",f.split(";")[0])
        geomPolygon = ogr.CreateGeometryFromWkt(f.split(";")[1])
        oFeaturePolygon.SetGeometry(geomPolygon)
        oLayer.CreateFeature(oFeaturePolygon)
    # 创建完成后，关闭进程
    oDS.Destroy()
    print ("数据集创建完成！")

    
if __name__ == '__main__':
    
    shpFile = "TestCreateShp.shp"
    NewShp(shpFile)

4. GDAL修改shp

#!/usr/bin/python
# -*- coding: UTF-8 -*-
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
def UpdateShp(_filename):
    # 支持中文路径
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    # 属性表字段支持中文
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 注册驱动
    ogr.RegisterAll()
    strDriverName = "ESRI Shapefile"
    oDriver = ogr.GetDriverByName(strDriverName)
    if oDriver == None:
        warninglist.append("驱动不可用："+strDriverName)
        return ("驱动不可用："+strDriverName)
    # 打开数据
    oDS = oDriver.Open(_filename, 1)
    oLayer = oDS.GetLayer()
    # 编辑属性字段
    for feature in oLayer: 
        feature.SetField('NAME', "EditName")
        oLayer.SetFeature(feature)
        feature.Destroy()
    oDS.Destroy()
    print ("编辑完成！")
    
if __name__ == '__main__':
    
    shpFile = "TM_WORLD_BORDERS-0.3.shp"
    UpdateShp(shpFile)