[Python] GDAL/OGR操作矢量数据(shp、GeoJSON)

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

【相关链接】

1.GDAL/OGR例子合集
2.GDAL官网API
3.网上的学习笔记,教程来自于犹他州立大学
4.中国OSGeo教程

OGR中的数据结构

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常量
IntegerOFTInteger
List of integersOFTIntegerList
Floating point numberOFTReal
List of floating point numbersOFTRealList
StringOFTString
List of stringsOFTStringList
DateOFTDate
Time of dayOFTTime
Date and timeOFTDateTime
创建字段

在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)
【作用】用新字段的规则替换现有的字段
【参数说明】

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

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

需要更改的字段属性OGR常量
Field name onlyALTER_NAME_FLAG
Field type onlyALTER_TYPE_FLAG
Field width and/or precision onlyALTER_WIDTH_PRECISION_FLAG
All of the aboveALTER_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常量
PointwkbPoint
MultipointwkbMultiPoint
LinewkbLineString
MultilinewkbMultiLineString
PolygonwkbPolygon
MultipolygonwkbMultiPolygon
Unknown geometry typewkbUnknown(图层如果有多种几何类型就返回wkbUnknown)
No geometrywkbNone

查询

属性查询

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到磁盘
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页