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

RS迷途小书童

已于 2024-08-08 09:48:59 修改

阅读量1k

点赞数 6

分类专栏： Python&Remote Sensing 文章标签： python 开发语言遥感影像栅格裁剪栅格数据处理

于 2024-07-10 17:23:23 首次发布

本文链接：https://blog.csdn.net/m0_56729804/article/details/140328692

版权

Python&Remote Sensing 专栏收录该内容

26 篇文章

订阅专栏

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 代码部分

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

我这里计算的是NDWI，并对计算后的结果进行二值化。原理是一样的，读的时候按照滑块的大小读取，写入的时候同样。

# -*- 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在处理大数据集时提供了实用参考，助力提升工作效率与数据处理能力。

20240807更新......

上面的思路是通过分块然后进行整块的影像计算，今天更新一下分块之后读取每个像素的波段值，有助于进行深度学习训练，填补上面代码无法分波段读取的空隙。这两个一起基本解决了所有分块处理大型影像的问题，无论是做指数计算、分类、训练模型都能用上。

如果大家感兴趣，可以留言交流，整体思路逻辑并不复杂，逐行分解就很好读懂了！

"""
@Time ： 2023/6/25 16:28
@Auth ： RS迷途小书童
@File ：RF_Predict_Personal.py
@IDE ：PyCharm
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""


    for y in range(0, ds_height, block_size):
        for x in range(0, ds_width, block_size):
            project_memory, system_memory = Memory()
            print(r"正在处理:第%s行,第%s列  程序内存占用:%smb, 系统剩余:%smb......" % (y+1, x+1, project_memory, system_memory))
            del project_memory, system_memory
            x_size = min(block_size, ds_width - x)  # 判断块和剩余的大小，避免右侧超出范围
            y_size = min(block_size, ds_height - y)
            array = ds.ReadAsArray(x, y, x_size, y_size).astype(np.float64)
            array_block = list()
            for i in range(0, y_size):
                for j in range(0, x_size):
                    array_band = list()
                    for b in range(0, ds_bands):
                        array_band.append(array[b][i][j])
                    array_block.append(array_band)
                    del array_band
            #  result = model.predict(array_block)
            predicted_image = result.reshape(y_size, x_size)
            del array_block, array, i, j, b, result, y_size, x_size
            ds_result.GetRasterBand(1).WriteArray(predicted_image, x, y)  # 写入数据
            del predicted_image