GDAL——栅格数据的裁剪，矢量数据的裁剪和按属性导出

1. 栅格数据的裁剪

需要注意的小细节：

• 核心函数是gdal.Warp()
• Warp()的属性中output_raster和input_raster可以使路径也可以是Dataset,而cutlineDSName只能是路径。
• srcBand属性：表示只裁剪第几个波段，不设置则裁剪全部波段。
• cutlineWhere和cutlineSQL属性：可以对矢量数据进行简单的筛选

from osgeo import gdal
import numpy as np

gdal.UseExceptions()


if __name__=="__main__":
    input_raster=r"待裁剪文件的路径"
    output_raster=r"输出路径"
    shp_path=r"矢量边界路径"


    gdal.Warp(
        output_raster,              # 路径 or Dataset
        input_raster,               # 路径 or Dataset
        srcBand=1,                  # 只取第1个波段（波段编号从 1 开始）

        # 对裁剪矢量的要素进行筛选
        # cutlineWhere="名称 = '上游'",
        # cutlineSQL="SELECT * FROM '1黄河流域分区范围' WHERE 名称 = '上游'",

        cutlineDSName=shp_path,     # 裁剪的矢量文件路径
        cropToCutline=True,         # 裁剪到矢量边界精确范围
        dstNodata=np.nan,           # 空缺值设置
        format="GTiff"
    )

2.矢量数据的按属性导出

实质：创建一个新文件 → 筛选符合条件的要素 → 写入新文件
要点：

• 中文字段的问题难以解决，特别是在SetAttributeFilter()函数中
• 设置驱动：driver = ogr.GetDriverByName("ESRI Shapefile")
• 创建源文件：driver.CreateDataSource()
• 创建图层：dst_layer = dst_ds.CreateLayer()
• 设置图层字段：dst_layer.CreateField(layer_defn.GetFieldDefn(i))
• 添加元素：dst_layer.CreateFeature(feature)

from osgeo import ogr
from osgeo import gdal

gdal.UseExceptions()

input_shp = r"E:\DroughtAnalysis\data\研究区\2020年黄河流域矢量范围数据\Data\1黄河流域分区\1黄河流域分区范围.shp"
output_shp = r"E:\DroughtAnalysis\data\SPEI_raster\gdal_spei12\1黄河流域分区范围select.shp"


# ========== 按属性筛选 ==========
*实质*：读取源文件 → 判断条件 → 将符合条件的要素复制 → 写入新文件


# 1. 打开源文件
src_ds = ogr.Open(input_shp)
src_layer = src_ds.GetLayer()   # .shp文件一个图层；.gpkg文件多个图层；.gdb文件多个图层

# 2. 查看属性字段（先确认字段名）
layer_defn = src_layer.GetLayerDefn()
# print("属性字段列表:")
# for i in range(layer_defn.GetFieldCount()):
#     print(f"  [{i}] {layer_defn.GetFieldDefn(i).GetName()}")

# 3. 设置属性过滤（假设字段名为“分区”或“名称”，请根据实际修改）
# src_layer.SetAttributeFilter("名称 = '上游'")      # 中文字段名，会报错
# src_layer.SetAttributeFilter("name = 'upstream'")  # 英文字段名
src_layer.SetAttributeFilter("name = '上游'")  # 示例

# 4. 创建输出文件
driver = ogr.GetDriverByName("ESRI Shapefile")  #"ESRI Shapefile"驱动
if driver.TestCapability(ogr.ODrCCreateDataSource):     # 有些 Driver 只能读取，不能创建新文件；需要判断是否可写新文件
    # 关键：添加编码选项
    dst_ds = driver.CreateDataSource(output_shp, options=["ENCODING=UTF-8"])    # options=["ENCODING=UTF-8"]为设置中文编码
    dst_layer = dst_ds.CreateLayer(
        "upstream",                 # 图层名称，如果是.shp文件里添加layer，则无意义
        src_layer.GetSpatialRef(),  # 坐标系信息
        src_layer.GetGeomType(),     # 几何类型（点、线、面、多面）
        options = ["ENCODING=UTF-8"]    # 会警告，但是不加会因为中文字段名报错，必须在创建源文件和图层是都添加
    )

    # 复制字段定义
    for i in range(layer_defn.GetFieldCount()):
        dst_layer.CreateField(layer_defn.GetFieldDefn(i))

    # 复制过滤后的要素
    for feature in src_layer:
        dst_layer.CreateFeature(feature)
        print(f"已复制要素 ID: {feature.GetFID()}")

    dst_ds = None
    src_ds = None

    print(f"\n 按属性筛选完成！输出文件: {output_shp}")
else:
    print("无法创建输出文件")

3. 矢量数据的裁剪

from osgeo import gdal

# 设置中文编码（避免乱码）
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")

# 路径设置
input_sites = r""
clip_boundary = r""
output_clipped = r""


print("正在裁剪站点数据...")

# 一行代码完成裁剪
gdal.VectorTranslate(
    output_clipped,  # 输出文件
    input_sites,  # 输入文件（被裁剪的）
    format="ESRI Shapefile",  # 输出格式
    clipSrc=clip_boundary,  # 裁剪边界(只接受路径，不接受图层对象)
    layerCreationOptions=["ENCODING=UTF-8"]  # 避免中文乱码
)


print(f" 裁剪完成！输出文件: {output_clipped}")