osgeo：OGR

OGR数据结构层次

from osgeo import ogr
path = './osgeo_data/Test.shp'

一、驱动(driver)$\Rightarrow$数据源(datasource)$\Rightarrow$图层(layer)

• driver：[ogr.GetDriverByName(‘ESRI Shapefile’)]
  • Open(self, path, update=0)
    • self $\Rightarrow$ [ogr] or [driver]
    • path $\Rightarrow$ 绝对路径
    • update: 0 $\Rightarrow$ 只读, 1 $\Rightarrow$ 可写
• datasource：[ogr.Open()] or [driver.Open()]
  • GetDriver：获取驱动
  • GetLayer：获取图层
  • GetLayerCount：图层个数
  • GetLayerByName | GetLayerByIndex：获取指定名称/索引图层

datasource = ogr.Open(path)
driver = datasource.GetDriver()
# driver = ogr.GetDriverByName('ESRI Shapefile')
# datasource = driver.Open(path, 0)
driver.name, datasource.GetRefCount()

(‘ESRI Shapefile’, 1)

二、图层$\Rightarrow$图层字段结构(layerdefn)|要素(feature)

• layer
  • GetExtent：图层范围
  • GetMetadata：图层信息
  • GetSpatialRef：图层空间参考

  ---

  • GetLayerDefn：获取图层结构

  ---

  • GetFeatureCount：图层要素个数
  • GetFeature：获取图层指定要素
  • GetNextFeature：获取图层下一个可用的要素

  ---

• layerdefn:[Layer.GetLayerDefn()]
  • GetName：所有字段名称
  • GetWidth：所有字段宽度
  • GetType：所有字段类型
  • GetPrecision：所有字段的格式精度
  • GetFieldDefn：所有字段个数
  • AddFieldDefn：添加FieldDefn

layer = datasource.GetLayer(0)
layer.GetLayerDefn(), layer.GetFeature(0)
layer.GetExtent(), layer.GetMetadata(), layer.GetSpatialRef()

((376523.5365316905, 406542.0224471837, 4453450.847271131, 4479300.099031694),
{‘DBF_DATE_LAST_UPDATE’: ‘2023-03-07’},
<osgeo.osr.SpatialReference; proxy of <Swig Object of type ‘OSRSpatialReferenceShadow *’ at 0x000001BB901566F0> >)

layerdefn = layer.GetLayerDefn()
for i in range(layerdefn.GetFieldCount()):
    defn = layerdefn.GetFieldDefn(i)
    print(defn.GetName(), defn.GetWidth(), defn.GetType(), defn.GetPrecision())

Id 6 0 0

三、要素$\Rightarrow$字段(field)|形状(geometry)

• feature

  • GetField([‘key’/index])：要素字段值

  • GetFieldCount：要素字段个数

  • keys：所有字段名称

  • GetGeometryRef：要素几何形状

• field

• geometry

  • GetGeometryName：几何要素名称
  • GetGeometryCount：几何要素个数
  • GetPointCount：点的个数
  • ExportToWkt：将几何图形转换为文本格式
  • GetSpatialReference：空间参考
  • GetEnvelope：几何要素外接矩形地理范围(不是整个图层)

feature = layer.GetFeature(0)
feature.keys(), feature.GetField('Id'), feature.GetFieldCount()

([‘Id’], 0, 1)

geom = feature.GetGeometryRef()
geom.GetGeometryName(), geom.GetGeometryCount(), geom.GetPointCount()

(‘POLYGON’, 1, 0)

geom.GetGeometryType(), geom.ExportToWkt()[:60]

(3, ‘POLYGON ((377649.229753523 4460455.16065141,378316.307218312’)

geom.GetEnvelope()

(376523.5365316905, 406542.0224471837, 4453450.847271131, 4479300.099031694)

四、创建shapefile

• 自定义创建
• 拷贝

4.1 基本流程

0. 删除数据源： driver.DeleteDataSource：
   • 判断filename是否存在（若存在，则删除）
1. 创建数据源： driver.CreateDataSource(filename):
   • filename不能存在
2. 创建图层： datasource.CreateLayer(filename, srs=None, geom_type, options=None)
   • srs: SpatialReference,
   • geom_type: OGRwkbGeometryType
   • options: Options are driver specific
3. 创建字段结构： layer.CreateField(fiedlDefn) | CreateFields([fiedlDefn1, fiedlDefn2])
4. 创建要素：
   • ogr.Featrue() $\Rightarrow$ ogr.Geometry() $\Rightarrow$ geometry.SetGeometry() $\Rightarrow$ feature.SetField() $\Rightarrow$ layer.CreateFeature()

import os
from osgeo import ogr

# 如果shape文件存在，则通过DeleteSource删除
driver = ogr.GetDriverByName('ESRI Shapefile')
path2 = './osgeo_data/Test_new.shp'
if os.path.exists(path2):
    driver.DeleteDataSource(path2)

ds = driver.CreateDataSource(path2)

layer = ds.CreateLayer('test_new', geom_type=ogr.wkbPoint)

# 创建字段结构
fieldDefn = ogr.FieldDefn('id', ogr.OFTString)
fieldDefn.SetWidth(4)
layer.CreateField(fieldDefn)  # CreateFields([fieldDefn1, fieldDefn2])

# 从layer中读取相应的feature类型，并创建feature
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
# 设定几何形状
point = ogr.Geometry(ogr.wkbPoint)
point.SetPoint(0, 123, 123)  # SetPoint(int point, double x, double y, double z=0)
feature.SetGeometry(point)
# 设定某字段数值
feature.SetField('id', 23)
# 将feature写入layer
layer.CreateFeature(feature)
# 最后，从内存中清除 ds，将数据写入磁盘中
ds.Destroy()

4.2 创建要素

4.2.1 创建点要素（Point）

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10, 20)
print(point)

POINT (10 20 0)

4.2.2 创建线要素（Line|Ring）

# 创建线要素
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10, 10), line.AddPoint(20, 20)
print(line)
# 创建线环要素Ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(0, 0), ring.AddPoint(100, 0)
ring.AddPoint(100, 100), ring.AddPoint(0, 100)
ring.CloseRings()
print(ring)

LINESTRING (10 10 0,20 20 0)
LINEARRING (0 0 0,100 0 0,100 100 0,0 100 0,0 0 0)

4.2.3 创建面要素（Polygon）

# 创建多边形要素
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.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)
print(polygon)
# polygon.GetGeometryCount(), polygon.GetGeometryRef(0)

POLYGON ((0 0 0,100 0 0,100 100 0,0 100 0,0 0 0),(25 25 0,75 25 0,75 75 0,25 75 0,25 25 0))

4.2.4 创建复合几何形状（MultiGeometry）

multipoint = ogr.Geometry(ogr.wkbMultiPoint)
point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(10, 10)
multipoint.AddGeometry(point1)
point1.AddPoint(20, 20)
multipoint.AddGeometry(point1)
multipoint.GetGeometryCount()

2

4.3 拷贝

• 数据源层次：driver.CopyDataSource(ds, outputfile) $\Rightarrow$ out.Release()
• 图层层次：driver.CreateDataSource(outfile).CopyLayer(layer, outputfile) $\Rightarrow$ out.Destroy()
• 要素层次：driver.CreateDataSoure().CreateLayer().

4.3.1 数据源层次

from osgeo import ogr
import os, math
inshp = './osgeo_data/Test.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = './osgeo_data/Test_copy_ds.shp'
# access:检查用户是否具有某项权限（如打开文件的权限）,mode=os.F_OK/R_OK/W_OK/X_OK(存在性、可读性、可写性和可执行性)
if os.access(outputfile, os.F_OK):
    driver.DeleteDataSource(outputfile)
ds2 = driver.CopyDataSource(ds, outputfile)
ds2.Release()

4.3.2 图层层次

inshp = './osgeo_data/Test.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = './osgeo_data/Test_copy_layer.shp'
if os.access(outputfile, os.F_OK):
    driver.DeleteDataSource(outputfile)
layer = ds.GetLayer()
ds2 = driver.CreateDataSource(outputfile)
pt_layer = ds2.CopyLayer(layer, 'abcd')
ds2.Destroy()

4.3.3 要素层次

inshp = './osgeo_data/Test.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = './osgeo_data/Test_copy_feature.shp'
if os.access(outputfile, os.F_OK):
    driver.DeleteDataSource(outputfile)
newds = driver.CreateDataSource(outputfile)
layer = ds.GetLayer()
layernew = newds.CreateLayer('copy', srs=layer.GetSpatialRef(), geom_type=ogr.wkbLineString)
feature = layer.GetNextFeature()
while feature is not None:
    layernew.CreateFeature(feature)
    feature = layer.GetNextFeature()
newds.Destroy()

五、过滤器

• SQL查询：ExecuteSQL(DataSource self, char const * statement, Geometry spatialFilter=None, char const * dialect)

5.1 空间过滤（Spatial filters）

• 几何过滤：layer.SetSpatialFilter(geom)
• 坐标过滤：layer.SetSpatialFilterRect(minx, miny, maxx, maxy)

import os
from osgeo import ogr, osr

# 将叠合的要素生成一个shp文件
def create_shp_by_layer(shp, layer):
    outputfile = shp
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputfile):
        driver.DeleteDataSource(outputfile)
    ds = driver.CreateDataSource(outputfile)
    ds.CopyLayer(layer, '')
    ds.Destroy()

    
world_shp = './osgeo_data/Test.shp'
cover_shp = './osgeo_data/Test_copy_ds.shp'
ds_world = ogr.Open(world_shp)
ds_cover = ogr.Open(cover_shp)
layer_world = ds_world.GetLayer(0)
layer_cover = ds_cover.GetLayer(0)

cover_feats = layer_cover.GetNextFeature()
poly = cover_feats.GetGeometryRef()  # 获取几何参考
layer_world.SetSpatialFilter(poly)
out_shp = './osgeo_data/world_cover.shp'
create_shp_by_layer(out_shp, layer_world)

5.2 属性过滤

• SetAttributeFilter(Layer self, filter_string)
  • filter_string: char, eg.“AREA > 800000”

ds = ogr.Open('./osgeo_data/Test.shp')
layer = ds.GetLayer(0)
layer.SetAttributeFilter('id==0')
# 根据属性条件生成要素
driver = ogr.GetDriverByName('ESRI Shapefile')
outfile = './osgeo_data/filter_out.shp'
if os.access(outfile, os.F_OK):
    driver.DeleteDataSource(outfile)
newds = driver.CreateDataSource(outfile)
newlayer = newds.CreateLayer('rect', srs=None, geom_type=ogr.wkbPolygon)
feat = layer.GetNextFeature()
while feat is not None:
    newlayer.CreateFeature(feat)
    feat = layer.GetNextFeature()
newds.Destroy()

六、空间计算

多边形（Polygon）：

• 交：poly3.Intersection(poly2)
• 并：poly3.Union(poly2)
• 差：poly3.Difference(poly2)
• 补：poly3.SymmetricDifference(poly2)

几何形状(Geometry)：

• 缓冲区：<geom>.Buffer(<distance>)
• 判断相等：<geom1>.Equal(<geom2>)
• 最短距离：<geom1>.Distance(<geom2>)
• 相接矩形4角坐标：<geom>.GetEnvelope()
• 面积：<geom>.GetArea()

判断两个对象的关系:

• Intersect判断两个要素是否相交：poly2.Intersect(poly1)
• Disjoint是判断两个要素是否不相交: poly2.Disjoint(poly1)
• Touch表示相邻：poly2.Touches(poly1)
• Crosses表示横穿：poly2.Crosses(line)
• Within表示包含：ptB.Within(poly1)
• Contains也表示包含，跟Within正好相反：poly1.Contains(ptB)
• Overlaps重叠（只有两个多边形之间才能算作overlap）：poly2.Overlaps(poly3)