Update - Added additional method to compare with.
Using a mixture of arcpy
and numpy
the methods below should do the trick:
import arcpy, numpy
##Get Rasters bands as well as raster height, width, and no data value
in_path = r'' #enter the fullpath to your raster stack here
arcpy.env.workspace = in_path
bands = arcpy.ListRasters()
in_ras = arcpy.Raster(in_path)
raster_rows = in_ras.width
raster_columns = in_ras.height
no_data = in_ras.noDataValue # this property can be tricky so you may need to set the value manually
#create an empty numpy array of zeros the same size as input raster to store the results
out_raster_data = numpy.zeros((raster_columns, raster_rows), numpy.int)
From here I will present two different ways you can go.
The first method (shown below) relies on native numpy functions. It uses native methods from numpy to mask data based on the provided no data value. This method is much faster:
bands_list = []
for i in bands:
temp_array = arcpy.RasterToNumPyArray(i).astype(numpy.float32)
masked_array = numpy.ma.masked_values(temp_array, no_data)
bands_list.append(masked_array.mask) #use the mask from masked_array that is an array of True/False values depending on if the value in the array was no_data or not
for i in xrange(0, len(bands)):
out_raster_data += bands_list[i] #perform addition on the Boolean mask arrays (1 for False (no_data), 0 for True (data))
The second method is a looping method which loops over the pixels/array elements individually. This method is much slower.
#Create a list that houses each raster array
bands_list = []
for i in bands:
bands_list.append(arcpy.RasterToNumPyArray(i).astype(numpy.float32))
#loop over all columns, rows, then bands
#loop over bands last to get a value from each band and increase the output pixel value if it is equal to no_data
for y in xrange(0, raster_columns):
for x in xrange(0, raster_rows):
out_pixel_data = 0
for j in range(0,len(bands)):
sample_value = bands_list[j][y,x]
if sample_value == no_data:
out_pixel_data += 1
out_raster_data[y,x] = out_pixel_data
To save the output in raster format for both methods, you can finish with this:
#convert numpy array to raster and save the results
result_raster = arcpy.NumPyArrayToRaster(out_raster_data)
result_raster.save(r'') #enter output raster path here
If you are working with a spatial reference you will need to project the resulting raster either when converting to raster or after the fact.
numpy.ndarray
or do you first need to read it from disk? How is your data structured?