关键词：gdal

graph LR
  Dataset --(包含)--> Layer
  Layer --(包含)--> Feature
  1. Dataset（数据集）：是 GDAL/OGR 中数据的最高容器，代表一个完整的数据集（可以是矢量或栅格数据）。
    • 对于矢量数据：包含一个或多个Layer（图层），并管理这些图层的空间参考、元数据等全局信息。
    • 对于栅格数据：包含栅格波段、投影、地理变换等信息。
  2. Layer（图层）：是Dataset中包含的矢量数据层，用于组织同一类型的空间要素。一个Dataset可以包含多个Layer
    • 所有要素具有相同的几何类型（如点、线、面）。
    • 共享同一套属性字段结构（如 “名称”“面积” 等字段）。
    • 共享相同的空间参考（投影坐标系）。
  3. Feature（要素）：是Layer中的基本空间单元，代表现实世界中的一个地理实体（如一个城市、一条道路）。
    • 由几何信息和属性信息组成。
      • 几何信息（Geometry）：点、线、面等空间位置和形状。
      • 属性信息：键值对形式的非空间数据（如名称、人口、面积等）。

打开关闭数据源

1
2
3
4
5
6
7
8
9
# 打开影像文件
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly, GA_Update

src_ds = gdal.Open(file_path, GA_ReadOnly)  # 只读模式
src_ds = gdal.Open(file_path, GA_Update)  # 可写模式

# 关闭数据集
del dataset  # 释放资源

栅格数据操作

获取基本信息

1
2
3
4
5
6
7
src_ds.RasterXSize       # 宽度
src_ds.RasterYSize       # 高度
src_ds.RasterCount       # 波段数
src_ds.GetProjection()   # 投影
src_ds.GetGeoTransform() # 地理变换
src_ds.GetMetadata()     # 元数据

波段操作

1
2
3
4
5
6
n_bands = src_ds.RasterCount

# 波段从1开始计数
band = src_ds.GetRasterBand(idx)

band_data = band.ReadAsArray()

创建数据集

1
2
3
4
5
6
7
8
9
10
11
12
13
# TIFF格式
driver = gdal.GetDriverByName(
    'GTiff'
)
    
# 创建数据集
dst_ds = driver.Create(
    output_path,
    width,
    height,
    bands,
    dtype
) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Shapefile格式
driver = ogr.GetDriverByName(
    "ESRI Shapefile"
)

# 判断是否已经存在相同的shp
if os.path.exists(file_path):
    driver.DeleteDataSource(
        file_path
    )
    
# 创建数据源
ds = driver.CreateDataSource(
    file_path
)
1
2
3
4
5
6
7
8
# dtype:
gdal.GDT_Byte    # 无符号8位整数。
gdal.GDT_UInt16  # 无符号16位整型
gdal.GDT_Int16   # 有符号16位整型
gdal.GDT_UInt32  # 无符号32位整型
gdal.GDT_Int32   # 有符号32位整型
gdal.GDT_Float32 # 32位浮点型
gdal.GDT_Float64 # 64位浮点型

设置投影和地理变换

1
2
3
4
5
projection = src_ds.GetProjection()
dst_ds.SetProjection(projection)

geotransform = src_ds.GetGeoTransform()
dst_ds.SetGeoTransform(geotransform)

读写数据

1
2
3
4
src_data = src_ds.ReadAsArray() # 读取图像像素
src_data = src_ds.GetRasterBand(idx).ReadAsArray()

dst_ds.GetRasterBand(idx).WriteArray(data) # 写图像像素

Shapefile操作

打开Shapefile

1
2
3
4
5
6
7
from osgeo import gdal, ogr, osr
from osgeo.gdalconst import GA_ReadOnly, GA_Update

# 打开Shapefile

driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(file_path, GA_Update)

获取Shapefile信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
layer = src_ds.GetLayer(0)
srs = layer.GetSpatialRef()
geom_type = layer.GetGeomType()
feature_count = layer.GetFeatureCount()
    
# 获取属性字段信息
fields = []
layer_defn = layer.GetLayerDefn()
for i in range(layer_defn.GetFieldCount()):
    field_defn = layer_defn.GetFieldDefn(i)
    fields.append({
        "name": field_defn.GetName(),
        "type": ogr.GetFieldTypeName(field_defn.GetType())
    })
     
info = {
    "空间参考": srs.ExportToWkt() if srs else "未知",
    "几何类型": ogr.GeometryTypeToName(geom_type),
    "要素数量": feature_count,
    "属性字段": fields
}

# 向图层添加属性字段
fd_defn = ogr.FieldDefn(field_name, field_type)
layer.CreateField(field_defn)

Layer操作

创建图层

1
2
3
4
5
6
7
8
9
10
11
# 设置空间参考(默认WGS84)
srs = osr.SpatialReference()
srs.ImportFromWkt(projection)
#srs.ImportFromEPSG(4326)  # WGS84坐标系

layer = ds.CreateLayer(
    layer_name,
    srs=srs,
    geom_type=geom_type  # 如ogr.wkbPoint, ogr.wkbPolygon等
)

从数据源中获取图层

1
2
layer = ds.GetLayerByName(layer_name)  # 通过图层名字
layer = ds.GetLayer(layer_index)       # 通过图层索引

获取图层基本信息

1
2
3
4
5
6
7
8
layer.GetName()         # 名称
layer.GetFeatureCount() # 要素数量
layer.GetSpatialRef().ExportToProj4() if layer.GetSpatialRef() else None # 空间参考
layer.GetExtent()       # 范围: (min_x, max_x, min_y, max_y)

layer_defn = layer.GetLayerDefn()
ogr.GeometryTypeToName(layer_defn.GetGeomType()) # 几何类型
layer_defn.GetFieldCount()                       # 字段数量

获取字段

1
2
for i in range(layer_defn.GetFieldCount()):
    field_defn = layer_defn.GetFieldDefn(i)

添加字段

1
2
3
4
5
6
7
8
9
10
field_defn = ogr.FieldDefn(field_name, field_type)
field_defn.SetWidth(width)
field_defn.SetPrecision(precision)
    
# 检查字段是否已存在
if layer.FindFieldIndex(field_name, True) != -1:
    raise Exception("字段已存在")
    
if layer.CreateField(field_defn) != 0:
    raise Exception(f"无法添加字段: {field_name}")

删除字段

1
2
3
4
5
6
7
8
index = layer.FindFieldIndex(field_name_or_index, True)
if index == -1:
    raise Exception("字段不存在")
else:
    index = field_name_or_index
    
if layer.DeleteField(index) != 0:
    raise Exception("无法删除字段")

要素操作

创建要素

1
2
layer_defn = layer.GetLayerDefn()
feature = ogr.Feature(layer_defn)

读取要素

1
2
3
4
5
for feature in layer:
    ...

# 通过ID获取单个要素
layer.GetFeature(fid)

设置要素属性

1
feature.SetField(key, value)

将要素添加到图层

1
layer.CreateFeature(feature)

设置要素几何

1
feature.SetGeometry(geometry)

删除要素

1
2
fid = feature.GetFID()
layer.DeleteFeature(fid) # -2表示所有

获取要素属性

1
2
3
4
5
attributes = {}
layer_defn = layer.GetLayerDefn()
for i in range(layer_defn.GetFieldCount()):
		field_name = layer_defn.GetFieldDefn(i).GetName()
		attributes[field_name] = feature.GetField(i)

释放资源

1
feature.Destroy()

创建几何

点

1
2
3
4
5
6
7
8
9
10
11
12
13
from osgeo import ogr

# 创建点几何（WGS84坐标：经度116.4，纬度39.9）
point = ogr.Geometry(ogr.wkbPoint)  # 初始化点对象
point.AddPoint(116.4, 39.9)  # 添加坐标（X, Y）

# 可选：添加Z值（三维点）
point_3d = ogr.Geometry(ogr.wkbPoint25D)
point_3d.AddPoint(116.4, 39.9, 500)  # （X, Y, Z）

# 查看几何信息
print("点坐标：", point.GetX(), point.GetY())
print("点WKT格式：", point.ExportToWkt())  # WKT格式便于查看和存储

线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 创建线几何（由多个点组成）
line = ogr.Geometry(ogr.wkbLineString)

# 依次添加线的顶点（顺序决定线的走向）
line.AddPoint(116.3, 39.8)  # 起点
line.AddPoint(116.4, 39.9)
line.AddPoint(116.5, 40.0)  # 终点

# 可选：三维线
line_3d = ogr.Geometry(ogr.wkbLineString25D)
line_3d.AddPoint(116.3, 39.8, 100)
line_3d.AddPoint(116.4, 39.9, 200)

# 查看线信息
print("线顶点数量：", line.GetPointCount())
print("线长度：", line.Length())  # 计算线的长度（单位取决于空间参考）

面

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 创建面几何（单个外环）
polygon = ogr.Geometry(ogr.wkbPolygon)

# 创建外环（线性环）
outer_ring = ogr.Geometry(ogr.wkbLinearRing)
# 依次添加外环顶点（最后一点需与第一点重合，确保闭合）
outer_ring.AddPoint(116.3, 39.8)
outer_ring.AddPoint(116.5, 39.8)
outer_ring.AddPoint(116.5, 40.0)
outer_ring.AddPoint(116.3, 40.0)
outer_ring.AddPoint(116.3, 39.8)  # 闭合环

# 将外环添加到面
polygon.AddGeometry(outer_ring)

# 可选：添加内环（孔洞）
inner_ring = ogr.Geometry(ogr.wkbLinearRing)
inner_ring.AddPoint(116.35, 39.85)
inner_ring.AddPoint(116.45, 39.85)
inner_ring.AddPoint(116.45, 39.95)
inner_ring.AddPoint(116.35, 39.95)
inner_ring.AddPoint(116.35, 39.85)  # 闭合

polygon.AddGeometry(inner_ring)  # 面可以包含多个环（1个外环+N个内环）

# 查看面信息
print("面面积：", polygon.Area())  # 计算面积（单位取决于空间参考）
print("面环数量：", polygon.GetGeometryCount())  # 外环+内环数量

多部件几何

1
2
3
4
5
6
7
8
9
10
# 多面点（MultiPoint）
multipoint = ogr.Geometry(
    ogr.wkbMultiPoint
)
p1 = ogr.Geometry(ogr.wkbPoint)
p1.AddPoint(116.4, 39.9)
p2 = ogr.Geometry(ogr.wkbPoint)
p2.AddPoint(121.4, 31.2)
multipoint.AddGeometry(p1)
multipoint.AddGeometry(p2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 多线（MultiLineString）
multiline = ogr.Geometry(
    ogr.wkbMultiLineString
)
line1 = ogr.Geometry(
		ogr.wkbLineString
)
line1.AddPoint(116.3, 39.8)
line1.AddPoint(116.4, 39.9)
line2 = ogr.Geometry(
		ogr.wkbLineString
)
line2.AddPoint(116.5, 40.0)
line2.AddPoint(116.6, 40.1)
multiline.AddGeometry(line1)
multiline.AddGeometry(line2)
1
2
3
4
5
6
7
8
9
10
# 多面（MultiPolygon）
multipolygon = ogr.Geometry(
	ogr.wkbMultiPolygon
)
poly1 = ogr.Geometry(ogr.wkbPolygon)
# （省略poly1的外环定义，同单个面）
poly2 = ogr.Geometry(ogr.wkbPolygon)
# （省略poly2的外环定义）
multipolygon.AddGeometry(poly1)
multipolygon.AddGeometry(poly2)

将几何写入Shapefile

1
2
3
4
feature = ogr.Feature(layer.GetLayerDefn())
geom = ogr.Geometry(ogr.wkbPoint)
geom.AddPoint(116.4, 39.9)
feature.SetGeometry(geom)
许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。