tif数据转为numpy进行栅格计算

# -*- coding: utf-8 -*-
"""
tif数据转为numpy进行栅格计算
2024.01.12, 22:52, YMJ
"""

import rasterio
import numpy as np

def raster_to_numpy(input_raster):
    # 打开栅格数据
    with rasterio.open(input_raster) as src:
        # 读取栅格数据为 NumPy 数组
        raster_array = src.read(1)

    return raster_array

def raster_arithmetic_numpy(input_array1, input_array2, operation):
    # 执行加减运算
    if operation == 'add':
        result_array = np.add(input_array1, input_array2)
    elif operation == 'subtract':
        result_array = np.subtract(input_array1, input_array2)
    else:
        raise ValueError("无效的运算操作")

    return result_array

def numpy_to_raster(input_array, output_raster, reference_raster):
    # 使用 reference_raster 的信息创建输出栅格文件
    with rasterio.open(reference_raster) as src:
        profile = src.profile
        profile.update(count=1)

        # 将 NumPy 数组写入栅格文件
        with rasterio.open(output_raster, 'w', **profile) as dst:
            dst.write(input_array, 1)

raster1_path = "D:\\resample_test\\2\\2022_01_01_hhh_clip.tif"
raster2_path = "D:\\resample_test\\output_19790116_clip.tif"
output_path = "D:\\resample_test\\output1.tif"

# 转换为 NumPy 数组
array1 = raster_to_numpy(raster1_path)
array2 = raster_to_numpy(raster2_path)

# 执行加法和减法运算
result_array_add = raster_arithmetic_numpy(array1, array2, 'add')
result_array_subtract = raster_arithmetic_numpy(array1, array2, 'subtract')

# 将 NumPy 数组写入栅格文件
numpy_to_raster(result_array_add, output_path, raster1_path)
# numpy_to_raster(result_array_subtract, output_raster_subtract, input_raster1)