GDAL-栅格数据和矢量数据的创建和修改

1.栅格数据的创建

顺序：选择驱动gdal.GetDriverByName("GTiff")→创建空数据ds = driver.Create(output_raster, cols, rows, bands, dtype)→写入数据ds.GetRasterBand(1).WriteArray(data)→写入地理变换参数ds.SetGeoTransform(geotransform)→写入投影信息ds.SetProjection(srs.ExportToWkt())

import numpy as np
from osgeo import gdal, osr

# ========== 1. 创建栅格文件 ==========
output_raster = r"E:\DroughtAnalysis\data\SPEI_raster\gdal_spei12\created_raster.tif"

# 参数设置
cols = 100      # 列数（X方向像元数）
rows = 80       # 行数（Y方向像元数）
bands = 1       # 波段数
dtype = gdal.GDT_Float32  # 数据类型

# 创建空文件
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create(output_raster, cols, rows, bands, dtype)

# ========== 2. 写入数据（NumPy数组） ==========
# 创建示例数据：一个渐变的值数组
data = np.zeros((rows, cols), dtype=np.float32)
for i in range(rows):
    data[i, :] = i  # 每行值等于行号

# 写入第1波段
ds.GetRasterBand(1).WriteArray(data)
ds.GetRasterBand(1).SetDescription("elevation")  # 波段描述

# ========== 3. 设置地理变换参数 ==========
# 左上角坐标、像元宽度、旋转(0)、左上角Y、旋转(0)、像元高度(负值)
x_min = 110.0
x_res = 0.1  # 像元宽度（度）
y_max = 35.0
y_res = -0.1  # 像元高度（负值，表示Y轴向下）

geotransform = (x_min, x_res, 0, y_max, 0, y_res)
ds.SetGeoTransform(geotransform)

# ========== 4. 设置投影（WGS84） ==========
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # WGS84
ds.SetProjection(srs.ExportToWkt())

# ========== 5. 关闭文件 ==========
ds = None

print(f" 栅格文件创建完成: {output_raster}")

2.矢量数据的创建和修改

顺序：选择驱动→创建数据源→设置坐标系→创建图层→添加字段→创建要素并写入图层

from osgeo import ogr, osr

# ========== 1. 创建矢量文件 ==========
output_shp = r"E:\DroughtAnalysis\data\SPEI_raster\gdal_spei12\created_sites.shp"

# 获取驱动
driver = ogr.GetDriverByName("ESRI Shapefile")

# 创建数据源
ds = driver.CreateDataSource(output_shp)

# ========== 2. 设置坐标系（WGS84） ==========
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# ========== 3. 创建图层 ==========
layer = ds.CreateLayer("sites", srs, ogr.wkbPoint)

# ========== 4. 添加字段 ==========
# ID字段（整数）
field_id = ogr.FieldDefn("id", ogr.OFTInteger)
layer.CreateField(field_id)

# 名称字段（字符串，长度50）
field_name = ogr.FieldDefn("name", ogr.OFTString)
field_name.SetWidth(50)
layer.CreateField(field_name)

# 值字段（浮点数）
field_value = ogr.FieldDefn("value", ogr.OFTReal)
layer.CreateField(field_value)

# ========== 5. 创建要素（添加点） ==========
# 定义三个点的坐标和属性
points = [
    {"lon": 110.5, "lat": 35.2, "id": 1, "name": "站点A", "value": 12.5},
    {"lon": 112.3, "lat": 34.8, "id": 2, "name": "站点B", "value": 8.3},
    {"lon": 108.7, "lat": 33.5, "id": 3, "name": "站点C", "value": 15.2},
]

for p in points:
    # 创建几何点
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(p["lon"], p["lat"])

    # 创建要素
    feature = ogr.Feature(layer.GetLayerDefn())
    feature.SetGeometry(point)
    feature.SetField("id", p["id"])
    feature.SetField("name", p["name"])
    feature.SetField("value", p["value"])

    # 写入图层
    layer.CreateFeature(feature)

    print(f"  已添加: {p['name']} ({p['lon']}, {p['lat']})")

# ========== 6. 关闭文件 ==========
ds = None

print(f"\n 矢量文件创建完成: {output_shp}")