This script takes the BQA band of all of the Nyriagongo and Kawa Ijen imagery and creates a binary mask where 1 is a non-clouded pixel and 0 is a clouded pixel. Then, it goes through every band of the source imagery (except BQA) and cuts out the pixels that are clouded. All results are rendered to the output directory.

cloud-hole-all.py

# cloud-hole-all.py

# Nicholas Hallahan nick@theoutpost.io
# Sun Nov 24 2013

# This script takes the BQA band of all of the Nyriagongo and Kawa Ijen imagery and
# creates a binary mask where 1 is a non-clouded pixel and 0 is a clouded pixel.

# Then, it goes through every band of the source imagery (except BQA) and cuts out
# the pixels that are clouded. All results are rendered to the output directory.

import os
import arcpy

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


cwd = os.getcwd()
nyriagongoPath = r"C:\vbox-shared\remote-sensing\nyiragongo"
kawahIjenPath = r"C:\vbox-shared\remote-sensing\kawah-ijen"


######################## create output paths #################################
outputPath = cwd + '\\output'
if not os.path.exists(outputPath): os.makedirs(outputPath)
outputNyiragongoPath = outputPath + '\\nyiragongo'
if not os.path.exists(outputNyiragongoPath): os.makedirs(outputNyiragongoPath)
outputKawahIjenPath = outputPath + "\\kawah-ijen"
if not os.path.exists(outputKawahIjenPath): os.makedirs(outputKawahIjenPath)

# Input either nyriagongoPath or kawahIjenPath
def createOutputLCDirectories(inputBasePath, outputBasePath):
	for f in os.listdir(inputBasePath):
		path = inputBasePath + "\\" + f
		if os.path.isdir(path) and f[:2] == 'LC':
			outputLCPath = outputBasePath + '\\' + f
			if not os.path.exists(outputLCPath): os.makedirs(outputLCPath)

createOutputLCDirectories(nyriagongoPath, outputNyiragongoPath)
createOutputLCDirectories(kawahIjenPath, outputKawahIjenPath)
##############################################################################


def cutHoleInBand(inputBandFile, bqaFile, outputBandFile):
	print("Cutting holes in " + inputBandFile)

	bqa = Raster(bqaFile)

	r61440 = bqa == 61440
	r59424 = bqa == 59424
	r57344 = bqa == 57344
	r56320 = bqa == 56320
	r53248 = bqa == 53248
	r52256 = bqa == 52256
	r52224 = bqa == 52224
	r49184 = bqa == 49184
	r49152 = bqa == 49152
	r48128 = bqa == 48128
	r45056 = bqa == 45056
	r43040 = bqa == 43040
	r39936 = bqa == 39936
	r36896 = bqa == 36896
	r36864 = bqa == 36864
	r32768 = bqa == 32768
	r31744 = bqa == 31744
	r28672 = bqa == 28672
	r16380 = bqa == 16380
	r13246 = bqa == 13246

	rSum = r61440 + r59424 + r57344 + r56320 + r53248 + r52256 + r52224 + r49184 + r49152 + r48128 + r45056 + r43040 + r39936 + r36896 + r36864 + r32768 + r31744 + r28672 + r16380 + r13246

	# we want all the pixels that are NOT cloud to be flagged as 1
	rSum = 1 - rSum
	
	holyBand = Raster(inputBandFile) * rSum
	holyBand.save(outputBandFile)
	print("Holy band saved to " + outputBandFile)


# We give this the main directory for Nyriagongo / Kawah Ijen for both
# input and output
def processBands(inputPath, outputPath):
	# Here we are iterating through all of the LC directories and executing 
	# cutHoleInBand
	lcDirs = os.listdir(outputPath) # we know (hope) the output dirs dont have extra crap

	for lcDir in lcDirs:
		inputLCPath = inputPath + '\\' + lcDir
		outputLCPath = outputPath + '\\' + lcDir
		bands = os.listdir(inputLCPath)

		# find the BQA band in the directory
		for band in bands:
			if band[-7:] == 'BQA.TIF':
				bqaPath = inputLCPath + '\\' + band
				break

		# now we can process all the bands that are not BQA
		for band in bands:
			if band.split('.')[-1] == 'TIF' and band[-7:] != 'BQA.TIF':
				inputBandPath = inputLCPath + '\\' + band
				outputBandPath = outputLCPath + '\\' + band
				
				cutHoleInBand(inputBandPath, bqaPath, outputBandPath)


# Ready... Set... GO!!!!
processBands(nyriagongoPath, outputNyiragongoPath)
processBands(kawahIjenPath, outputKawahIjenPath)