GDAL-栅格数据重采样，波段运算

1. 栅格数据重采样

from osgeo import gdal

gdal.UseExceptions()

input_raster = r""
output_raster = r""

# 目标分辨率（度）
target_xres = 0.1
target_yres = 0.1

gdal.Warp(
    output_raster,
    input_raster,
    xRes=target_xres,
    yRes=target_yres,
    resampleAlg=gdal.GRA_Bilinear,  # 连续数据用双线性插值
    format="GTiff"
)

2. 波段运算

import os
import numpy as np
from osgeo import gdal

# 1. 读取多个时间的数据
path = r""
spei_1 = gdal.Open(os.path.join(path, "1.tif")).ReadAsArray()
spei_2 = gdal.Open(os.path.join(path, "2.tif")).ReadAsArray()
spei_3 = gdal.Open(os.path.join(path, "3.tif")).ReadAsArray()

# 2. 计算平均
spring_mean = (spei_1 + spei_2 + spei_3) / 3

# 3. 计算像元值 < -1）
drought_mask = spring_mean < -1
drought_count = np.sum(drought_mask)
print(f"筛选像元数: {drought_count}")

# 4. 将结果写入新文件
driver = gdal.GetDriverByName("GTiff")  # GeoTIFF 格式的驱动
out = driver.Create("spring_spei.tif", spring_mean.shape[1], spring_mean.shape[0], 1, gdal.GDT_Float32) # 创建一个空的 GeoTIFF 文件
out.GetRasterBand(1).WriteArray(spring_mean)    # 将 NumPy 数组写入文件的第 1 个波段（从1开始）。

# 复制地理信息
ds = gdal.Open("spei_200001.tif")   # 打开一个已有的参考栅格文件（用于获取地理信息）
out.SetGeoTransform(ds.GetGeoTransform())   # 地理变换参数（左上角坐标、像元大小、旋转参数）
out.SetProjection(ds.GetProjection())   # 投影信息
ds = None
out = None