Calculating the mean value of a raster stack using numpy

Question

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.