Rasterio卫星图像预处理执行时间优化及工具选择咨询

一、优化执行时间的核心建议 #

1. 用Numpy向量化运算彻底替换Python嵌套循环 #

你的代码中最影响性能的是双重Python循环遍历像素——Python的解释型循环在处理百万级像素时效率极低，而Numpy的向量化运算基于C实现，能带来几十到上百倍的速度提升。

修改后的dn_to_radiance：

import numpy as np

def dn_to_radiance(data_array, band_number):
    # 注意：提前处理MTL字典的键，去掉末尾空格（见下文优化建议2）
    channel_gain = float(Landsat8_mlt_dict[f'RADIANCE_MULT_BAND_{band_number}'])
    channel_offset = float(Landsat8_mlt_dict[f'RADIANCE_ADD_BAND_{band_number}'])
    
    # 直接用数组向量化运算，自动遍历所有元素
    new_data_array = data_array * channel_gain + channel_offset
    # 保留原数组的NaN值（无需手动判断，Numpy会自动处理）
    new_data_array[np.isnan(data_array)] = np.nan
    
    print(f'Radiance calculated for band {band_number}')
    return new_data_array

修改后的radiance_to_reflectance：

import numpy as np
import math

def radiance_to_reflectance(arr, ESUN):
    d = float(Landsat8_mlt_dict['EARTH_SUN_DISTANCE'])
    sun_elevation = float(Landsat8_mlt_dict['SUN_ELEVATION'])
    phi_rad = math.radians(90 - sun_elevation)
    
    # 预计算常量，避免重复计算
    constant = (np.pi * (d ** 2)) / (ESUN * math.cos(phi_rad))
    new_data_array = arr * constant
    new_data_array[np.isnan(arr)] = np.nan
    
    print(f"Reflectance of Band calculated")
    return new_data_array

2. 预处理MTL参数，减少字典查找开销 #

原代码中MTL字典的键带有末尾空格（比如'RADIANCE_MULT_BAND_1 '），每次查找都要拼接带空格的字符串，既容易出错又浪费性能。读取MTL时直接清理键的空格：

Landsat8_mlt_dict = {}
with open('LC08_L2SP_190037_20190619_20200827_02_T1_MTL.txt', 'r') as f:
    for line in f:
        line = line.strip()
        if line == 'END':
            break
        key, value = line.split('=')
        # 清理键的前后空格，避免后续查找出错
        Landsat8_mlt_dict[key.strip()] = value.strip()

3. 用Rasterio分块处理超大图像 #

如果卫星图像分辨率极高（比如Landsat 8的单景图像超过1GB），一次性加载到内存会导致卡顿甚至崩溃。用Rasterio的block_windows分块读取处理，既能节省内存，又能利用CPU缓存提升效率：

import rasterio

with rasterio.open('your_image.tif') as src:
    with rasterio.open('output.tif', 'w', **src.profile) as dst:
        # 遍历图像的每个分块
        for ji, window in src.block_windows():
            # 读取分块数据
            data = src.read(window=window)
            # 对分块应用你的计算逻辑
            processed_data = dn_to_radiance(data, band_number=1)
            # 将处理后的分块写入输出文件
            dst.write(processed_data, window=window)

4. 优化数据类型减少计算开销 #

检查输入数组的数据类型：如果原数据是uint16（Landsat的DN值常用类型），计算时可以转换为float32而非默认的float64，在精度足够的前提下，减少一半内存占用和计算时间：

# 在读取数据时指定 dtype
data_array = src.read(1, dtype=np.float32)

二、运行环境选择 #

本地运行 #

• 适合场景：图像数据量大、需要频繁读写本地文件、对数据隐私要求高。
• 优势：如果你的本地机器有多核心CPU或独立GPU，可以配合Numba（CPU即时编译）或CuPy（GPU加速Numpy）进一步提升性能；无需上传数据，避免网络延迟。

Jupyter/Google Colab #

• 适合场景：快速调试代码、尝试不同参数、分享分析成果。
• 优势：
  • Jupyter的交互式环境能快速查看中间结果，适合频繁试错的阶段；
  • Google Colab提供免费GPU/TPU资源，能大幅加速大规模数组运算；
  • 无需配置本地环境，直接在线运行代码。
• 注意：如果图像超过Colab的内存限制，需要结合分块处理逻辑；数据上传到云端会有隐私风险，敏感数据不建议使用。

内容的提问来源于stack exchange，提问作者Salah Amani