0

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

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

for band in range(1, bands+1):
    data = ndvi_stack.GetRasterBand(band).ReadAsArray().astype('float')
    mean = np.mean(data[data != 0]) #calculate mean without value 0
1

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

#read ndvi_stack layers
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

#read ndvi_stack layers
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)
ndvi_array = np.ma.array(ndvi_array, mask=(ndvi_array == 0))  # masked array
result = np.ma.mean(ndvi_array, axis=0)  # (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.

  • Thank you so much it worked a charm!! – Marlin May 17 at 13:13

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.