[Python] GDAL/OGR操作矢量数据（shp、GeoJSON）

最新推荐文章于 2024-07-30 10:51:07 发布

geodoer

最新推荐文章于 2024-07-30 10:51:07 发布

阅读量1.6w

点赞数 29

分类专栏： GIS实践文章标签： OGR GDAL python

本文链接：https://blog.csdn.net/summer_dew/article/details/87930241

版权

GIS实践专栏收录该内容

86 篇文章 10 订阅

订阅专栏

GDAL项目旨于地理数据抽象模型对地理数据文件进行读写管理；而其项目下有两大类模块：GDAL和OGR
OGR提供操作矢量数据的API，GDAL模块提供栅格数据的API

OGR中的数据结构

driver解析数据的驱动

# 【方法一】从已有数据源中获取驱动变量
ds = ogr.open(r'D:\....\...shp')
driver = ds.GetDriver()
# 【方法二】通过名称获取
json_driver = ogr.GetDrverByName('GeoJSON')
    # 获得驱动名的两种方法
    # 1. OGR网站上有介绍，通过GDAL/ORG自带的ogrinfo
    # 2. 代码中提供的print_drivers函数来获取驱动程序的名字

datasource管理数据源

（1）创建数据源
有了驱动程序之后，提供数据源名称来使用它创建一个空的数据源。新创建的数据源会自动打开等待写入

# 创建一个功能齐全的SpatialLite数据源，而不是使用SQLite
fn = "filename"

driver = ogr.GetDriverByName("SQLite")
# 创建新的数据源时，不能覆盖现有的数据源。如果你的代码可能会覆盖现有数据，那么在创建新数据之前需要删除旧数据
if os.path.exists(fn):
	driver.DeleteDataSource(fn)
ds= driver.CreateDataSource(fn, 
	options=[ #创建选项：传递一个字符串列表（参数在OGR网站上有文档介绍）
		'SPATIALITE=yes'
	]
)
if ds is None:
	#【注意】如果数据源没有创建成功，那么CreateDataSource函数会返回为空。如果为空对象，之后使用会报AttributeError错误
	sys.exit('Could not create {0}'.formaat(json_fn ) )

（2）删除数据源

datasource.Destory() #关闭
del datasource # 删除datasource
	""" 删除ds变量强制关闭文件，并将所有的编辑成果写入磁盘中
	【注意】删除图层变量并不会触发这个操作，必须关闭数据源才行
	【备注】如果想保持数据源打开的话，可以通过图层对象或者数据源对象调用ds.SyncToDisk()
	【警告】为了使你的编辑写入到磁盘中，必须关闭文件或者调用SyncToDisk函数。如果没有这么做，并且在交互环境中还打开数据源，那么你会很失望地发现创建了一个空的数据集
	"""

layer管理图层

layer = datasource.GetLayer(0) #根据下标来获取图层（对于shp而言只有一个图层）
layer.GetGeomType() #图层几何类型
n = layer.GetFeatureCount() #要素数量
extent = layer.GetExtent() #上下左右边界
readedNum = GetFeaturesRead() #已经读取多少条Feature
layer.ResetReading() #重置内部feature指针，指向FID=0的要素
layer.GetSpatialRef() #使用wkt文本表示空间参考系统的实例
layer.schema #获取FieldDefn对象列表（属性名）

feature管理要素类

从Layer中获得feature

# 1. 【获得图层中的要素】
feature = layer.GetFeature(0) 
# 2. 【获取图层中所有要素】
feat = layer.GetNextFeature()
while feat:
	feat = layer.GetNextFeature()
layer.ResetReading()

更新要素

要素更新操作也大同小异，只不过你操作的是图层中的现有要素，而不是一个空白的对象

# 【例子】为图层的每一个要素添加一个唯一的ID
layer.CreateField( # 为图层添加一个ID字段
	ogr.FieldDefn('ID', ogr.OFTInteger)
) 
n = 1
for feat in layer:
	feat.SetField('ID', n) #设置字段的值
	layer.SetFeature(feat) #将数值传递给SetFeature函数，来更新图层中的要素信息（这里不是CreateFeature函数）
	n += 1

删除要素

删除要素更容易，需要知道的就是待删除要素的FID

#【例子】删除字段City_Name为"Seattle"的要素
# 找到目标的FID进行删除
for feat in layer:
	if feat.GetField('City_Name')=='Seattle':
		layer.DeleteFeature(feat.GetFID() ) #删除操作

"""
【但是】
1. 某些特定的数据格式使用此方式并不能完全删除要素。
2. 你可能看不出，有时要素只是被标记为删除了，而不是被完全抛出，它仍然还潜伏着。
3. 正因如此，你不会看到其他要素被分配到刚才删除的FID
4. 这意味着，如果已经删除了很多要素，在文件中可能存在大量不必要的已用空间
【所以】删除这些要素可以回收对应的空间
【如何回收对应空间】如果你有一些关系数据库的经验，应该熟悉这一点。它类似于在微软的Access数据库上运行压缩和修复或在PostgreSQL数据库中使用重组功能（VACUUM）
【具体做法】打开数据源，然后执行一条SQL语句来压缩数据库（不同的数据格式回收的方法不同）
"""
ds.ExecuteSQL('REPACK' + layer.GetName() ) #shapefile格式回收方式
ds.ExecuteSQL('RECOMPUTE EXTENT ON' + layer.GetName() ) #确保空间范围更新
	# 当现有要素发生改变或被删除后，它不会更新元数据中的空间范围

关闭要素

feature.Destory() #关闭

创建要素

往现有图层中添加新要素和往全新的图层中添加要素的操作一样，具体请看示例

【步骤】
1.基于图层字段创建一个空要素，填充它
2.然后把它插入到该图层

geometry管理要素的几何属性

Geometry的几何属性

【org.Geometry中的几何类型】ogr.wkbPoint、org.wkbLineString、org.wkbPolygon等，请看wkt

# 获得Geometry的几何属性
layer.GetGeomType() == ogr.wkbPolygon
layer.GetGeometryName()

从Feature获取Geometry

# 两种方法
geometry = feature.geometry()
geometry = feat.GetGeometryRef()

常用方法

geom.GetGeometryName() #要素类型
geom.GetPointCount() #要素点个数
geom.GetX(0);geom.GetY(0);geom.GetZ(0) #获得第一个坐标点的X、Y、Z
print(geom) #打印出所有点
geom.ExportToWkt() #导出WKT

创建Geometry

请看示例

field管理要素的属性字段

【字段类型常量】

有一些不同的属性字段类型，但并不是每一种字段类型都支持所有的数据格式。这时候用于描述各种数据格式的在线文档就派上用场了。

字段数据类型	OGR常量
Integer	OFTInteger
List of integers	OFTIntegerList
Floating point number	OFTReal
List of floating point numbers	OFTRealList
String	OFTString
List of strings	OFTStringList
Date	OFTDate
Time of day	OFTTime
Date and time	OFTDateTime

创建字段

在Layer中新建属性字段

# 创建并添加第一个属性字段
coord_fld = ogr.FieldDefn('X', ogr.OFTReal)
# 设置限制
coord_fld.SetWidth(8) 
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
# 重用FieldDefn对象来创建第二个字段
coor_fld.SetName('Y')
out_lyr.CreateField(coord_fld)

更改属性定义

【函数】AlterFieldDefn(iField, field_def, nFlags)
【作用】用新字段的规则替换现有的字段
【参数说明】

iField是想改变的属性字段对应的索引值
field_def是新属性字段的定义对象
nFlags是一个整数，是下表中的一个或多个常数的总和

【用于指明字段定义的哪个属性可以更改的标记，如果想使用多个，一起添加即可】

需要更改的字段属性	OGR常量
Field name only	ALTER_NAME_FLAG
Field type only	ALTER_TYPE_FLAG
Field width and/or precision only	ALTER_WIDTH_PRECISION_FLAG
All of the above	ALTER_ALL_FLAG

# 【例子一】 更改原字段的名字
index = layer.GetLayerDefn().GetFieldIndex('NAME') #获取字段的索引值
fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString) #创建新属性的字段定义
layer.AlterFieldDefn(index, fld_defn, ogr.ALTER_NAME_FLAG) 

# 【例子二】 更改字段的多个属性，例如名称和精度
lyr_defn = layer.GetLayerDefn()
index = lyr_defn.GetFieldIndex('X')
width = lyr_defn.GetFieldDefn(index).GetWidth()
fld_defn = ogr.FieldDefn('X_coord', ogr.OFTReal)
fld_defn.SetWidth(width)
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
layer.AlterFieldDefn(index, fld_defn, flag)
	"""【注意】在此例子中创建新的字段定义时使用的是原始字段宽度
	如果没有为原始数据设置足够大的字段宽度，结果将不正确
	为了解决这个问题，需要使用与原始数据一样的字段宽度
	为了使更改的精度生效，所有记录必须重写
	为数据设置一个比它自身更高的精度并不会提高精度，因为数据不能凭空产生，但如果设置的精度不够高，精度就会降低
	"""

其他

feature.GetFieldCount() # 属性总数
feature.keys() # 属性名
feat.GetField('AREA') # 获取字段
feature.GetFieldAsString("FID") # 获取字段成string

经验

1.没能为GeoJSON文件设置一个精度
2.想在shapefile中设置字段精度，必须设置字段宽度
3.设置字段的默认值：FieldDefn.SetDefault("我是默认值")

wkt

WKB（二进制Well-KnownBinary），用于不同软件程序间进行几何要素类型转换的一种二进制表示标准。
因为它是二进制格式，所以人们无法直接阅读获取其表示的内容，但是熟知文本格式（Well-Known Text，WKT）可以阅读。

几何要素类型	对应OGR常量
Point	wkbPoint
Multipoint	wkbMultiPoint
Line	wkbLineString
Multiline	wkbMultiLineString
Polygon	wkbPolygon
Multipolygon	wkbMultiPolygon
Unknown geometry type	wkbUnknown（图层如果有多种几何类型就返回wkbUnknown）
No geometry	wkbNone

查询

属性查询

from osgeo import ogr
import os
shpfile = r'C:\tmp\data.shp'
ds = ogr.Open(shpfile)
layer = ds.GetLayer(0) #得到图层
lyr_count = layer.GetFeatureCount()
print(lyr_count) #原先的要素总数
layer.SetAttributeFilter("AREA > 800000") #属性查询条件
lyr_count = layer.GetFeatureCount() #查询过后，layer被更新
print(lyr_count) #满足条件的layer
# 将该layer保存成shp
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = r'C:\tmp\data.shp'
if os.access( extfile, os.F_OK ):
	driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
# 遍历，复制
feat = layer.GetNextFeature()
while feat is not None:
	layernew.CreateFeature(feat)
	feat = layer.GetNextFeature()
newds.Destroy()

空间查询

SQL查询

from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = r'C:\tmp\test.shp'
world_ds = ogr.Open(world_shp)
world_layer = world_ds.GetLayer()
world_layer_name = world_layer.GetName()
result = world_ds.ExecuteSQL("select * from %s where prov_id = '22' order by BNDRY_ID desc" % (world_layer_name)) # ) # ExecuteSQL是基于数据集进行的，而不是图层
resultFeat = result.GetNextFeature ()
out_shp = r'C:\tmp\test\test_sql_result.shp'
create_shp_by_layer(out_shp, result) #保存结果
# 对查询结果进行遍历
while resultFeat :
    print resultFeat.GetField('BNDRY_ID')
    resultFeat = result.GetNextFeature ()
#执行下一条SQL语句之前一定要先释放
world_ds.ReleaseResultSet(result)

其他

设置编码

# 全局设置
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") #支持文件名称及路径内的中文
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")        #支持属性字段中的中文
# 为图层指定编码
out_layer = ds.CreateLayer( "CommercialHousing",
    srs = in_layer.GetSpatialRef(),
    geom_type= ogr.wkbPoint,
    options  = [ 
        "ENCODING=UTF-8"
    ]
)

示例

打开shp

from osgeo import ogr
inshp_path = 'D:\mycode\GISandPython\0data\park_point_shp\xiamen_20181128_park.shp'
driver = ogr.GetDriverByName('ESRI Shapefile') #查找一个特定的驱动程序
datasource = driver.Open(inshp_path, 0) #0只读，1可写
dir( datasource) #使用Python的内省函数dir()查看所有方法

删除数据

# 删除矢量数据文件
if os.path.exists(out_shp)
	driver.DeleteDataSource(out_shp)
ds2 = driver.CreateDataSource(out_shp)

遍历所有属性值

#【遍历所有属性值】
for i in range(feature.GetFieldCount() ):
	print( feature.GetField(i) )
### 【查看表的结构，各个字段的名称等信息】在layer附加信息中看
layerdef = layer.GetLayerDefn()
for i in range(layerdef.GetFieldCount() ):
	defn = layerdef.GetFieldDefn(i)
	print(defn.GetName(), defn.GetWidth(), defn.GetType(), defn.GetPrecision() )

创建

创建点状数据集

from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'point_demo.shp'
point_coors = [[300,450], [750, 700], [1200, 450], [750, 200], [750, 450]]
print point_coors
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds  = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('point',None,ogr.wkbPoint)
for aa in point_coors:
    wkt = 'POINT (' + str(aa[0]) + ' ' + str(aa[1]) + ')'
    geom = ogr.CreateGeometryFromWkt(wkt)
    feat = ogr.Feature(layernew.GetLayerDefn())
    feat.SetGeometry(geom)
    layernew.CreateFeature(feat)
newds.Destroy()

创建线状数据集

from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'line_demo.shp'
point_coors = [300,450, 750, 700, 1200, 450, 750, 200]
print point_coors
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
newds  = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('point',None,ogr.wkbLineString)
wkt = 'LINESTRING (%f %f,%f %f,%f %f,%f %f,%f %f)' % (point_coors[0],point_coors[1],
    point_coors[2],point_coors[3], point_coors[4],point_coors[5],
    point_coors[6],point_coors[7], point_coors[0],point_coors[1])
print(wkt)
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
newds.Destroy()

创建多边形数据集

from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'rect_demo.shp'
if os.access( extfile, os.F_OK ):
 driver.DeleteDataSource( extfile )
extent = [400, 1100, 300, 600]
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
width = math.fabs(extent[1]-extent[0])
height = math.fabs(extent[3]-extent[2])
tw = width/2
th = width/2
extnew = extent[0]+tw
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
    extent[1],extent[3], extent[1],extent[2],
    extent[0],extent[2], extent[0],extent[3])
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
newds.Destroy()

创建属性

from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile") #shp驱动器
extfile = 'rect_field_demo.shp' #输出的shp路径
if os.access( extfile, os.F_OK ): #文件是否已存在
    driver.DeleteDataSource( extfile )

extent = [400, 1100, 300, 600] #范围
newds = driver.CreateDataSource(extfile) #按照驱动器创建数据源
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon) #创建图层
# 创建字段
fieldcnstr = ogr.FieldDefn("fd",ogr.OFTString) #创建字段(字段名，类型)
fieldcnstr.SetWidth(32) #设置宽度
layernew.CreateField(fieldcnstr) #将字段设置到layer
fieldf = ogr.FieldDefn("f",ogr.OFTReal) 
layernew.CreateField(fieldf)
# wkt格式的polygon
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
    extent[1],extent[3], extent[1],extent[2],
    extent[0],extent[2], extent[0],extent[3])
geom = ogr.CreateGeometryFromWkt(wkt) #根据wkt创建geometry
feat = ogr.Feature( layernew.GetLayerDefn() ) #从layer中得到feature的要求，创建一个Feature对象
feat.SetField('fd',"这里是字段的值") #设置字段值
feat.SetGeometry(geom) #设置要素的geometry属性
layernew.CreateFeature(feat) #在图层中创建该feature
newds.Destroy() #删除（包括保存）

创建点Geometry

from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20) #向Point添加多个点，不会报错，但最终只会用最后添加的点

创建线Geometry

from osgeo import ogr
line = ogr.Geometry(ogr.wkbLineString) # 创建
line.AddPoint(10, 20) # 添加点
line.SetPoint(1, 50, 50) # 修改点坐标
line.GetPointCount() # 得到所有点数目
line.GetX(0); line.GetY(0) # 读取0号点的坐标

创建多边形Geometry

from osgeo ipmort ogr
# 创建环
ring = ogr.Geometry(ogr.wkbLinearRing) #创建
ring.AddPoint(10,10) #添加点
ring.CloseRings() #闭合 或 添加的最后一个点和第一个相同
# 将环添加到多边形中
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(outring)
polygon.AddGeometry(inring)
ring = polygon.GetGeometryRef(0) #得到多边形内的几何对象

拷贝

datasource层次的拷贝

from osgeo import ogr
import os,math
inshp = "C:\tmp\test.shp"
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = 'C:\tmp\test_copy.shp'
if os.access( outputfile, os.F_OK) : #已经有该文件
	driver.DeleteDataSource(outputfile)  #删除
pt_cp = driver.CopyDataSource(ds, outputfile) #根据数据源ds创建数据，返回指针
pt_cp.Release() #释放指针，将数据写入到磁盘

layer层次的拷贝

from osgeo import ogr
import os,math
inshp = "C:\tmp\test.shp"
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = "C:\tmp\test_copy2.shp"
if os.access( outputfile, os.F_OK):
	driver.DeleteDataSourse( outputfile )
layer = ds.GetLayer()
newds = driver.CreateDataSource(outputfile) #首先得有数据源
pt_layer  = newds.CopyLayer(layer, 'abcd') #复制图层，返回指针
newds.Destroy() #对newds进行Destroy()操作，才能将数据写入磁盘

feature层次的拷贝

from osgeo import ogr
import os,math
inshp = 'C:\tmp\test.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = 'c:\tmp\test.shp'
if os.access( outputfile, os.F_OK ):
    driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile) #按outputfile的数据源创建一个newds
layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString) #创建一个图层
layer = ds.GetLayer() #得到原数据源的图层
extent = layer.GetExtent() #图层范围
# 遍历旧图层的所有feature，进行复制
feature = layer.GetNextFeature()
while feature is not None: 
    layernew.CreateFeature(feature)
    feature = layer.GetNextFeature()
newds.Destroy() # 包括了将数据flush到磁盘