前言

在遥感数据处理、空间分析以及指数计算等应用中，栅格计算是常用且必不可少的工具。根据不同的分析需求和基础数据源，栅格计算可以被细分为几种类型：首先是影像不同波段之间的栅格运算（例如NDVI计算等）；其次是不同影像之间的栅格运算，涉及坐标系、像元大小、图像范围的统一处理；最后是分区栅格计算，它特别适用于需要在不同区域应用不同计算系数的情形，比如在同一栅格计算公式下，不同区域可能需要乘以不同的系数。本文将提供实现上述三种栅格计算方法的代码和详细解析。

1.栅格计算

栅格计算包括基本的数学运算（加、减、乘、除等）和其他函数运算。进行栅格计算时，首先需要确保栅格的坐标系、范围和分辨率一致。如果我们对同一影像的不同波段进行计算，这些参数通常已经统一，无需额外考虑。然而，当进行不同影像之间的栅格计算时（例如两个不同的 tif 数据），就必须确保它们的坐标系、范围和分辨率相同，以保证计算结果的准确性。

1.1影像的不同波段之间的运算（以计算植被指数为例子）

在生态遥感中，同一影像的不同波段之间的运算是非常常见的，尤其是计算各类植被指数（如NDVI、NDWI、EVI等）。这些植被指数能够有效反映植被状况，广泛应用于环境监测、农作物生长评估等领域。以下代码示例以Landsat影像为例，展示如何计算多个植被指数，并将其结果叠加保存。

import rasterio
import numpy as np

file=r"C:\Users\HP\Desktop\repro\stack01.tif.tif"
output_file=r"C:\Users\HP\Desktop\repro\cal01.tif"
with rasterio.open(file) as src:
    profile = src.profile
    #更新即将输出结果的元数据,2个波段，防止输出结果取整数
    profile.update(count=2, dtype=rasterio.float32)
    # print(profile)
    red=src.read(3).astype(np.float32)#转为浮点型，防止计算的为整数
    NIR = src.read(4).astype(np.float32)
    NDVI=(NIR-red)/(NIR+red)
    NDVI[np.isnan(NDVI)] = -99#设置nodata值
    blue=src.read(1).astype(np.float32)
    evi = 2.5 * (NIR - red) / (NIR + 6 * red - 7.5 * blue + 1)
    evi[np.isnan(evi)] = -99
    with rasterio.open(output_file, "w", **profile) as dst:
        dst.write(NDVI,1)
        dst.write(evi, 2)#叠加输出，数字指的是第几波段

1.2不同影像之间的栅格运算

当进行不同来源影像之间的栅格运算时，我们首先需要确保它们具有一致的坐标系、像元大小和空间范围。这是因为不同影像的分辨率和坐标系可能存在差异，若不进行统一，会导致运算结果不准确。其大致步骤就是先考虑影像的坐标系是否一致，不一致先重投影，在对应裁剪至相同范围，在对影像重采样到同一像元大小。以符合计算要求。

'''
来源不同影像的栅格计算
'''
import rasterio
import geopandas as gpd
from rasterio.mask import mask
import numpy as np
from rasterio.io import MemoryFile

tif_ref = r"C:\Users\HP\Desktop\image_process\origin\st_rgb01.tif"
tif_src = r"C:\Users\HP\Desktop\image_process\origin\ldst_rgb01.tif"
shp_file= r"C:\Users\HP\Desktop\image_process\origin\range.shp"
#影像裁剪
def maskRaster(raster_file,vecror_file):
    with rasterio.open(raster_file) as src:
        # raster_data=src.read(1)
        gdf = gpd.read_file(vecror_file)
        shapes = [feature["geometry"] for feature in gdf.to_dict("records")]#循环geometry
        out_image, out_transform = mask(src, shapes, crop=True)
        profile=src.profile
        profile.update({
            "height": out_image.shape[1],
            "width": out_image.shape[2],
            "transform": out_transform
        })

        memfile = MemoryFile()#这里我是想到envi的memory file，就是文件并未实际保存
        with memfile.open(**profile) as mem_dst:
            mem_dst.write(out_image)
        return memfile
#影像重采样
def resize_raster(src_path, ref_path):
    with rasterio.open(ref_path) as ref:
        new_height, new_width = ref.height, ref.width

    with rasterio.open(src_path) as src:

        data = src.read(out_shape=(src.count, new_height, new_width), resampling=rasterio.enums.Resampling.bilinear)
        profile = src.profile
        profile.update({"height": new_height, "width": new_width})

        memfile = MemoryFile()
        with memfile.open(**profile) as mem_dst:
            mem_dst.write(data)
        return memfile

#栅格计算
def compute_index(img1_path, img2_path,output_path):
    with rasterio.open(img1_path) as img1, rasterio.open(img2_path) as img2:
        nir1, red1 = img1.read(1), img2.read(1)  # 读取波段

        index = nir1 + red1
        index[np.isnan(index)] = -99
        profile = img1.profile.copy()
        profile.update(dtype=rasterio.float32, count=1)

        with rasterio.open(output_path, "w", **profile) as dst:
            dst.write(index, 1)
clip_ref=maskRaster(tif_ref,shp_file)
clip_src=maskRaster(tif_src,shp_file)
resize_src=resize_raster(clip_src,clip_ref)

output_path=r"C:\Users\HP\Desktop\image_process\origin\test_index01.tif"
cal_index=compute_index(clip_ref,resize_src,output_path)

2.分区栅格计算

分区栅格计算是一种针对不同区域应用不同栅格计算方法的技术。在实际应用中，我们可能需要在同一栅格影像中，对不同区域应用不同的计算系数。例如，在计算某一省份的植被指数时，各市可能需要乘以不同的系数进行调整。在此过程中，我们首先需要获得各市的矢量数据，并为每个市添加相应的系数值。然后，通过将矢量数据栅格化，将系数值转换为栅格值，从而实现分区栅格计算。

import rasterio
import geopandas as gpd
import numpy as np
from rasterio.mask import mask
from rasterio.features import rasterize


index_raster_path = r"C:\Users\HP\Desktop\image_process\origin\st_rgb01.tif"  # 省级指数栅格
city_shapefile = r"C:\Users\HP\Desktop\image_process\origin\coef.shp"  # 市级矢量数据（包含 coefficient 字段）
output_raster_path = r"C:\Users\HP\Desktop\image_process\origin\adjusted_index.tif"  # 输出文件
#读取数据
with rasterio.open(index_raster_path) as src:
    profile = src.profile
    index_data = src.read(1)  # 读取指数数据
    transform = src.transform  # 获取地理变换信息

gdf = gpd.read_file(city_shapefile)

# 生成与指数栅格同尺寸的空矩阵用来存系数，这一步也是为了使矢量外的值设置为1做铺垫，防止最后得到栅格矢量外的数据为nan
coefficient_raster = np.full(index_data.shape, np.nan, dtype=np.float32)

# 遍历市级矢量，将 `coefficient` 栅格化
for index, row in gdf.iterrows():
    shape = row['geometry'] # 获取市的边界
    coef = row["coef"]  # 获取该市的系数

    # 创建栅格化掩膜（与指数栅格一致）
    mask_array = rasterize(
        [(shape, coef)],
        out_shape=index_data.shape,
        transform=transform,
        fill=np.nan,
        dtype=np.float32)

#rasterio.features.rasterize(
#    shapes,         # 需要转换的矢量数据 [(geometry, value)]
#    out_shape,      # 目标栅格的形状 (height, width)
#    transform,      # 目标栅格的地理变换参数（坐标与像素的对应关系）
#    fill,           # 栅格的默认填充值（没有匹配到矢量区域的像素赋值）
#    dtype           # 输出数据类型)

# 赋值给系数栅格
    coefficient_raster[~np.isnan(mask_array)] = mask_array[~np.isnan(mask_array)]
    #这个代码的意思市coefficient_raster与mask_array对应的非NAN区域等于mask_array的值，np.isnan(mask_array)是一个bool索引

# 处理无数据值（如果某些区域没有系数，默认设为 1），保证矢量外的区域是原栅格值，若直接使用mask_array矢量边界外的都为NAN
coefficient_raster[np.isnan(coefficient_raster)] = 1.0

# 计算调整后的指数栅格
adjusted_index = index_data * coefficient_raster

# 更新栅格文件元数据
profile.update(dtype=rasterio.float32, count=1, nodata=-99)

with rasterio.open(output_raster_path, "w", **profile) as dst:
    dst.write(adjusted_index, 1)