0

I am calculating the mean value of a stack of images. I created a virtual raster which is the stack of 68 sentinel-1 images. I want to calculate the mean value for every pixel.

Here is my code:

import os
import rasterio 
import numpy as np
import time 


outvrt = "stacked_s1.vrt"

in_ds = rasterio.open(outvrt)
print(in_ds.count)

array_mean = np.zeros(120560400, dtype=float)

for i in range(in_ds.count):
    start = time.time()
    index = i + 1
    band = in_ds.read(index)
    print(f"reading band index {i + 1}")
    flatten_band = band.flatten()
    for r in np.nditer(flatten_band):

        array_mean[r] = (array_mean[r] + r)/2
    end_time = time.time()
    hours, rem = divmod(end_time-start, 3600)
    minutes, seconds = divmod(rem,60)
    print("Finished visit of array in:")
    print(f"{int(hours)}:{int(minutes)}:{int(seconds)}")

print(len(array_mean))

As you can see, I created an zero array, that will hold the result of the calculation made in the iteration. Then I iterated over every band in the stack, and read the band as an numpy array. Afterwards, I flattened the array and used np.nditer to iterate over every value on the flattened array. I then, updated the values of the array_mean by adding the values of the array and divide it by 2. This is an error for the first iteration because the initial value in array_mean is 0. Therefore I inserted some if statements to check for this condition:

for r in np.nditer(flatten_band):
        if array_mean[r] == 0:
            array_mean[r] = r
        else:
            array_mean[r] = (array_mean[r] + r)/2

In the first code I had time completion of ~1.40 minutes per band. By introducing this extra condition checking I had a performance of ~2.30 minutes per band.

Should I get rid of calculating the mean during the iteration, and just adding the new values to create a cumulative sum array and then apply numpy.divide in the final array by the number of bands (68)?

This approach helps me to calculate the mean of a time series of data in my small laptop, and don't incur into Memory Errors by accessing values in a dense time series.

I did not write the code to reshape and write out the final result.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.