问题

I'm using Rasterio to work with a satellite image, and I need to iterate through the entire file. and applying the formula to each pixel. This process takes a long time and makes it difficult for me to try out different modifications because it takes a long time to see the results each time.

使用Rasterio处理卫星图像，我需要遍历整个文件并对每个像素应用公式。这个过程需要很长时间，使我难以尝试不同的修改，因为每次查看结果都需要很长时间。

any suggestions to improve time execution ?

有没有建议来改善执行时间？

And is it better to work on this project locally or via Jupiter, Google Colab, or other tools?

是在本地还是通过Jupyter、Google Colab或其他工具进行此项目更好？

def dn_to_radiance(data_array, band_number):
    # getting the G value
    channel_gain = float(Landsat8_mlt_dict['RADIANCE_MULT_BAND_' + str(band_number) + ' '])
    # Getting the B value
    channel_offset = float(Landsat8_mlt_dict['RADIANCE_ADD_BAND_' + str(band_number) + ' '])
    # creating a temp array to store the radiance value
    # np.empty_like Return a new array with the same shape and type as a given array.
    new_data_array = np.empty_like(data_array)

    # loooping through the image
    for i, row in enumerate(data_array):
        for j, col in enumerate(row):

            # checking if the pixel value is not nan, to avoid background correction
            if data_array[i][j].all() != np.nan:
                new_data_array[i][j] = data_array[i][j] * channel_gain + channel_offset
    print(f'Radiance calculated for band {band_number}')
    return new_data_array

Landsat8_mlt_dict = {}
with open('LC08_L2SP_190037_20190619_20200827_02_T1_MTL.txt', 'r') as _:
    for line in _:
        line = line.strip()
        if line != 'END':
            key, value = line.split('=')
            Landsat8_mlt_dict[key] = value

def radiance_to_reflectance(arr, ESUN, ):
    # getting the d value
    d = float(Landsat8_mlt_dict['EARTH_SUN_DISTANCE '])
    # calculating rh phi value from theta
    phi = 90 - float(Landsat8_mlt_dict['SUN_ELEVATION '])
    # creating the temp array
    new_data_array = np.empty_like(arr)
    # loop to finf the reflectance
    for i, row in enumerate(arr):
        for j, col in enumerate(row):
            if arr[i][j].all() != np.nan:
                new_data_array[i][j] = np.pi * arr[i][j] * d ** 2 / (ESUN * cos(phi * math.pi / 180))
    print("Reflectance of Band calculated")
    return new_data_array

def dn_to_radiance(data_array, band_number):
    # 获取G值
    channel_gain = float(Landsat8_mlt_dict['RADIANCE_MULT_BAND_' + str(band_number) + ' '])
    # 获取B值
    channel_offset = float(Landsat8_mlt_dict['RADIANCE_ADD_BAND_' + str(band_number) + ' '])
    # 创建一个临时数组来存储辐射度值
    # np.empty_like返回一个与给定数组具有相同形状和类型的新数组。
    new_data_array = np.empty_like(data_array)

    # 循环遍历图像
    for i, row in enumerate(data_array):
        for j, col in enumerate(row):

            # 检查像素值是否不是NaN，以避免背景校正
            if data_array[i][j].all() != np.nan:
                new_data_array[i][j] = data_array[i][j] * channel_gain + channel_offset
    print(f'已计算波段{band_number}的辐射度')
    return new_data_array

Landsat8_mlt_dict = {}
with open('LC08_L2SP_190037_20190619_20200827_02_T1_MTL.txt', 'r') as _:
    for line in _:
        line = line.strip()
        if line != 'END':
            key, value = line.split('=')
            Landsat8_mlt_dict[key] = value

def radiance_to_reflectance(arr, ESUN, ):
    # 获取d值
    d = float(Landsat8_mlt_dict['EARTH_SUN_DISTANCE '])
    # 从theta计算rh phi值
    phi = 90 - float(Landsat8_mlt_dict['SUN_ELEVATION '])
    # 创建临时数组
    new_data_array = np.empty_like(arr)
    # 循环查找反射率
    for i, row in enumerate(arr):
        for j, col in enumerate(row):
            if arr[i][j].all() != np.nan:
                new_data_array[i][j] = np.pi * arr[i][j] * d ** 2 / (ESUN * cos(phi * math.pi / 180))
    print("已计算波段的反射率")
    return new_data_array

def dn_to_radiance(data_array, band_number):
    # getting the G value
    channel_gain = float(Landsat8_mlt_dict['RADIANCE_MULT_BAND_' + str(band_number) + ' '])
    # Getting the B value
    channel_offset = float(Landsat8_mlt_dict['RADIANCE_ADD_BAND_' + str(band_number) + ' '])
    # creating a temp array to store the radiance value
    # np.empty_like Return a new array with the same shape and type as a given array.
    new_data_array = np.empty_like(data_array)

    # loooping through the image
    for i, row in enumerate(data_array):
        for j, col in enumerate(row):

            # checking if the pixel value is not nan, to avoid background correction
            if data_array[i][j].all() != np.nan:
                new_data_array[i][j] = data_array[i][j] * channel_gain + channel_offset
    print(f'Radiance calculated for band {band_number}')
    return new_data_array
Landsat8_mlt_dict = {}
with open('LC08_L2SP_190037_20190619_20200827_02_T1_MTL.txt', 'r') as _:
    # print(type(_))
    for line in _:
        line = line.strip()
        if line != 'END':
            key, value = line.split('=')
            Landsat8_mlt_dict[key] = value


# print(Landsat8_mlt_dict)
def radiance_to_reflectance(arr, ESUN, ):
    # getting the d value
    d = float(Landsat8_mlt_dict['EARTH_SUN_DISTANCE '])
    # calculating rh phi value from theta
    phi = 90 - float(Landsat8_mlt_dict['SUN_ELEVATION '])
    # creating the temp array
    new_data_array = np.empty_like(arr)
    # loop to finf the reflectance
    for i, row in enumerate(arr):
        for j, col in enumerate(row):
            if arr[i][j].all() != np.nan:
                new_data_array[i][j] = np.pi * arr[i][j] * d ** 2 / (ESUN * cos(phi * math.pi / 180))
    print(f"Reflectance of Band calculated")
    return new_data_array

答案1

得分: 0

你可以使用第三方库，比如EOReader，来为你将Landsat波段转换为反射率。

from eoreader.reader import Reader
from eoreader.bands import RED, GREEN

prod = Reader().open(r"LC08_L1TP_200030_20201220_20210310_02_T1.tar")

# 将这些波段加载为xarray.DataArray字典
bands = prod.load([RED, GREEN])
green = bands[GREEN]
red = bands[RED] 

声明：我是EOReader的维护者

如果你想自己做，你应该学习如何在Python中处理数组，绝对不要使用循环！相反，你应该矢量化你的计算，这样速度会快得多！

from eoreader.reader import Reader
from eoreader.bands import RED, GREEN

prod = Reader().open(r"LC08_L1TP_200030_20201220_20210310_02_T1.tar")

# Load those bands as a dict of xarray.DataArray
bands = prod.load([RED, GREEN])
green = bands[GREEN]
red = bands[RED]