# Create mean raster stack from raster stack

I am trying to create a mean stack image(1 panchromatic layer) from a stacked tiff image(13 panchromatic layers). How do i go about doing this?

My code:

``````from osgeo import gdal
import numpy as np

ndvi_stack = gdal.Open('NDVI_stack/ndvi_stack.tif')
bands = ndvi_stack.RasterCount

for band in range(1, bands+1):
mean = np.mean(data[data != 0]) #calculate mean without value 0
``````

In your code, you are calculating a mean value for each band. However you likely want to get a mean value for each set of values in the same pixel in your `ndvi_stack` image.

One way to accomplish this is to create a 3D `numpy array` using `np.stack()` and then calculating the mean by specifying the axis.

``````from osgeo import gdal
import numpy as np

ndvi_stack = gdal.Open('NDVI_stack/ndvi_stack.tif')
bands = ndvi_stack.RasterCount

ndvi_array = []  # list to store all the bands
for band in range(1, bands+1):
data = ndvi_stack.GetRasterBand(band).ReadAsArray().astype('float')  # (n rows by n cols array)
ndvi_array.append(data)

ndvi_array = np.stack(ndvi_array)  # (n bands by n rows by n cols array)
result = ndvi_array.mean(axis=0)  # (n rows by n cols array)
``````

I see you are trying to ignore zeros when calculating the mean. Following the same principle from the code above, you can use a masked numpy array to accomplish this. The only differences is that you mask the array to ignore zeros and calculate the mean using `np.ma.mean()` rather than `np.mean()`.

``````from osgeo import gdal
import numpy as np

ndvi_stack = gdal.Open('NDVI_stack/ndvi_stack.tif')
bands = ndvi_stack.RasterCount

ndvi_array = []  # list to store all the bands
for band in range(1, bands+1):
data = ndvi_stack.GetRasterBand(band).ReadAsArray().astype('float')  # (n rows by n cols array)
ndvi_array.append(data)

ndvi_array = np.stack(ndvi_array)  # (n bands by n rows by n cols array)
Furthermore, if there are pixels where all the bands have a value of 0, you will end up with masked pixels, represented with a `--` when you print the array. You can easily convert all the masked pixels to an arbitrary value (e.g. a NoData value) using: `result = result.filled(-99)`, replacing `-99` with the value of your choice.
After you have the result, which is a `numpy array`, you will have to create a raster dataset, write the array to the first band of the raster and save it.