【Python&RS】基于Python分块处理大型遥感影像的方法

​

        RSer工作时不可避免会用到大型的遥感影像，由于分辨率过高、区域过大、波段信息过多等原因，都会导致数据非常的大。这个时候我们在进行一些简单的操作，如计算NDVI、二值化、分类等时，计算机的内存都会溢出。因此今天跟大家分享一下我平时分块的方法，中间如何计算就按照自己的需求来即可。

原创作者：RS迷途小书童

博客地址：https://blog.csdn.net/m0_56729804?type=blog

1 忽略警告信息

        由于GDAL版本的更新，现在运行代码都会有一条警告Warning，但不影响程序，我有强迫症所以必须没有警告。另外就是指定Proj库的路径，以免坐标系进行警告。

gdal.DontUseExceptions()
os.environ['PROJ_LIB'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data/proj'
os.environ['GDAL_DATA'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data'

2 代码部分

        代码的思路就是按照设定好的块滑动处理影像，并按照原始的排列顺序将处理后的重新写入即可。中间计算部分可以按照自己需求来。

# -*- coding: utf-8 -*-
"""
@Time ： 2023/6/25 16:28
@Auth ： RS迷途小书童
@File ：Raster Data Block Processing.py
@IDE ：PyCharm
@Purpose：栅格数据分块处理
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import os
import numpy as np
from osgeo import gdal
gdal.DontUseExceptions()
os.environ['PROJ_LIB'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data/proj'
os.environ['GDAL_DATA'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data'


def block_big_image(input_path: str, output_path: str, block_size: int) -> None:
    """
    :param input_path: 输入需要处理的影像
    :param output_path: 输出路径
    :param block_size: 输入块的大小
    :return: None
    """
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)
    if not ds:  # 判断文件是否能打开
        print(f"无法打开文件 {input_path}")
        return
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    # ds_bands = ds.RasterCount  # 获取波段数
    # ds_type = gdal.GetDataTypeName(ds.GetRasterBand(1).DataType).lower()  # 获取位数信息
    ds_prj = ds.GetProjection()  # 获取投影信息
    driver = gdal.GetDriverByName('GTiff')  # 创建输出文件
    ds_result = driver.Create(output_path, ds_width, ds_height, 1, gdal.GDT_Float64,
                              options=["TILED=YES", "COMPRESS=LZW", "BIGTIFF=YES"])
    ds_result.SetGeoTransform(ds_geo)  # 设置仿射地理变化参数
    ds_result.SetProjection(ds_prj)  # 设置投影坐标信息

    # -----------------------------------------------逐块读取和写入数据----------------------------------------------------
    for y in range(0, ds_height, block_size):
        for x in range(0, ds_width, block_size):
            x_size = min(block_size, ds_width - x)  # 判断块和剩余的大小，避免右侧超出范围
            y_size = min(block_size, ds_height - y)  # 判断块和剩余的大小，避免右侧超出范围
            array_green = ds.GetRasterBand(2).ReadAsArray(x, y, x_size, y_size).astype(np.float64)
            array_nir = ds.GetRasterBand(4).ReadAsArray(x, y, x_size, y_size).astype(np.float64)
            b1 = array_green - array_nir
            b2 = array_green + array_nir
            array_result = np.divide(b1, b2, out=np.zeros_like(b2), where=b2 != 0)  # 相除，排除被除数为0的情况
            array_result = np.where(array_result > 0.2, 1, 0)  # 筛选array_result大于0.2的部分
            # array = np.where(array_result == 1, array, 0)  # 掩膜筛选，将array_result等于1的array部分保留，其余为0
            ds_result.GetRasterBand(1).WriteArray(array_result, x, y)  # 写入数据
            del array_green, array_nir, array_result, b1, b2, x_size, y_size
    ds_result.FlushCache()  # 清理缓存
    del ds, ds_width, ds_height, ds_prj


if __name__ == "__main__":
    # 定义输入和输出文件夹路径
    input_path1 = r'Y:\彭俊喜\1.tif'
    output_path1 = r'Y:\彭俊喜\2.tif'
    block_big_image(input_path1, output_path1, 256)

        在遥感影像处理中，面对高分辨率、大区域及多波段数据带来的巨大挑战，内存溢出成为常见难题。本博客旨在分享一种高效应对策略——影像分块处理法。通过将大型遥感影像分割成多个小块，独立处理后再整合结果，有效缓解了内存压力，使得NDVI计算、二值化、分类等操作得以顺畅进行。此方法灵活易行，为RSer在处理大数据集时提供了实用参考，助力提升工作效率与数据处理能力。

​