RSer工作时不可避免会用到大型的遥感影像,由于分辨率过高、区域过大、波段信息过多等原因,都会导致数据非常的大。这个时候我们在进行一些简单的操作,如计算NDVI、二值化、分类等时,计算机的内存都会溢出。因此今天跟大家分享一下我平时分块的方法,中间如何计算就按照自己的需求来即可。
原创作者:RS迷途小书童
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'
![点击并拖拽以移动](/i/l/?n=24&i=blog/3204676/202407/3204676-20240710172359817-1017749863.gif)
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)
![点击并拖拽以移动](/i/l/?n=24&i=blog/3204676/202407/3204676-20240710172359817-1017749863.gif)
在遥感影像处理中,面对高分辨率、大区域及多波段数据带来的巨大挑战,内存溢出成为常见难题。本博客旨在分享一种高效应对策略——影像分块处理法。通过将大型遥感影像分割成多个小块,独立处理后再整合结果,有效缓解了内存压力,使得NDVI计算、二值化、分类等操作得以顺畅进行。此方法灵活易行,为RSer在处理大数据集时提供了实用参考,助力提升工作效率与数据处理能力。
标签:block,分块,Python,RS,result,array,ds,gdal,size From: https://www.cnblogs.com/RSran/p/18294633