# Calculating the mean value of a raster stack using numpy

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.