This script averages together all of the RGB values from each image. Because the clouded pixels are set to 0, Those values are ignored in the average. We divide the sum of values by another raster that is the sum of valid pixels based on BQA classification.

average-rgb-nyira.py

# average-rgb2.py

# Nicholas Hallahan nick@theoutpost.io
# Sat Nov 30 2013

import os
import arcpy

from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True


cwd = os.getcwd()
inputPath = r"C:\vbox-shared\remote-sensing\code\cloud-hole-all\output\nyiragongo"
validPixelSum = r"C:\vbox-shared\remote-sensing\code\bqa\output-sum-not-cloud\LC81730612013322LGN00_BINSUM.TIF"
outputPath = r"C:\vbox-shared\remote-sensing\code\average-rgb\output\RGB_NYIRA_MEAN.TIF"
rOutputPath = r"C:\vbox-shared\remote-sensing\code\average-rgb\output\R_NYIRA_MEAN.TIF"
gOutputPath = r"C:\vbox-shared\remote-sensing\code\average-rgb\output\G_NYIRA_MEAN.TIF"
bOutputPath = r"C:\vbox-shared\remote-sensing\code\average-rgb\output\B_NYIRA_MEAN.TIF"

######################## list of input reds, greens, blues ##################
ls = os.listdir(inputPath)

reds = []
greens = []
blues = []
for lcDir in ls:
	lcDir = inputPath + '\\' + lcDir
	bands = os.listdir(lcDir)
	for band in bands:
		band = lcDir + '\\' + band
		if band[-6:] == 'B4.TIF':
			reds.append(band)
		elif band[-6:] == 'B3.TIF':
			greens.append(band)
		elif band[-6:] == 'B2.TIF':
			blues.append(band)

##############################################################################


print("Mean Averaging Nyiragongo RGBs")

rSum = Raster( reds.pop(0)   ) / 65000.0
gSum = Raster( greens.pop(0) ) / 65000.0
bSum = Raster( blues.pop(0)  ) / 65000.0

for r in reds:
	rSum += Raster(r) / 65000.0 # normalizing signed 32 bit int into float
for g in greens:
	gSum += Raster(g) / 65000.0
for b in blues:
	bSum += Raster(b) / 65000.0
	
# mean is based on the number of valid pixels for a given pixel
rMean = rSum / Raster(validPixelSum)
gMean = gSum / Raster(validPixelSum)
bMean = bSum / Raster(validPixelSum)

rMean.save(rOutputPath)
print("Red mean composite saved to " + rOutputPath)
gMean.save(gOutputPath)
print("Green mean composite saved to " + gOutputPath)
bMean.save(bOutputPath)
print("Blue mean composite saved to " + bOutputPath)

print("Compositing RGB bands into single RGB TIF")
arcpy.CompositeBands_management(rOutputPath+';'+gOutputPath+';'+bOutputPath, outputPath)
print("Composite RGB saved to " + outputPath)