rasterio data - python - execution time - preprocessing

Question

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. any suggestions to improve time execution ? And is it better to work on this project locally or via Jupiter, Google Colab, or other tools?

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

Answer 1

You could use third-party libraries such as EOReader to convert Landsat bands to reflectance for you.

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] 

Disclaimer: I am the maintener of EOReader

---

If you want to do that yourself, you should do some tutorials on how to handle arrays in Python. Never ever loop over them! You should instead vectorize your computations: it will go way faster!