;+
; NAME:
;   MODCldMskBat.pro
;
; PURPOSE:
;	This program is a lightweight utility to mask out cloudy pixels in a single MODIS L1B
;	HDF file. It is designed to be called multiple times from a unix script which will handle the file pre-duplication 
;	and output directory creation, etc...
;
; AUTHOR:
;	Alex Fraser
;
; CATEGORY:
;	Utility
;
;
; CALLING SEQUENCE:
;	see example
;
;
; INPUTS:
;	NONE: Automatically looks for files named ./output/L1B.hdf and ./output/cloud.hdf
;
; KEYWORDS:
;	/strict: if set, will mask if probably clear, probably cloudy and confident cloudy.
;	/smooth: if set, will smooth the cloud mask (default values are a 5*5 kernel, thresholding at 0.25)
;
;
; OUTPUTS:
;	Replaces the L1B SDSs within the HDF with the updated values masking cloudy pixels (masks all channels)
;
;
; RESTRICTIONS:
;	Currently only 1km L1B MODIS granules are supported
;
;
; EXAMPLE:
;	MODCldMskBat, /strict, /smooth, kernelSize=5, threshold=0.25
;
;
; MODIFICATION HISTORY:
;
;	Written by Alex Fraser, July/August and November/December 2007.
;-

PRO MODCldMskBat, strict=srct, smooth=smth, kernelSize=krnsz, threshold=thld

l1b='./output/L1B.hdf'
cloud='./output/cloud.hdf'

;read in cloud mask data - close cloud HDF
;size of this array is 1354*2030*6 bytes, reduced to 1354*2030*1 byte
;because we only need the first byte to summarize the presence/absence of cloud.
CloudData = HDFreturnSDS(cloud, 'Cloud_Mask')
CloudData = CloudData[*,*,0]


;process cloud mask info, editing L1B data as we go
; Algorithm from Ackerman1997:
; For no tolerance of cloud at all, the following algorithm is suitable:
; Read bit 0 to determine if cloud mask was calculated. If 0, mask not calculated.

; Generate mask bit pattern
if keyword_set(srct) then begin
	mask=6 ;Both need to be equal to 1 (confident cloud-free)
endif else begin
	mask=4 ;Only bit 2 needs to be 1 (probably cloud-free)
endelse 


; Make note of (in CloudDataCalculated) only those pixels for which the cloud mask was calculated (mask=1)
CloudDataCalculated = CloudData AND 1
CloudDataCloudFree = CloudData AND mask


;Now that we have the cloud free array (CloudDataCloudFree), smooth it then threshold
if n_elements(smth) ne 0 then begin
	CloudDataCloudFreeNorm = float(CloudDataCloudFree)/float(mask) ;(normalize it to a matrix of 1 (cloud free) and 0 (cloudy))

	;Apply a convolution to the matrix
	;First, assemble the kernel
	if n_elements(krnSz) eq 0 then krnSz=5
	if n_elements(thld) eq 0 then thld=0.25

	constant = 1/float(krnSz*krnSz)
	kernel = Replicate(constant,krnSz,krnSz) ;produces a kernelSize*kernelSize kernel filled with the value 1/kernelSize*kernelSize
	filtered = Convol(CloudDataCloudFreeNorm, kernel, Center=1)


	;now threshold this "filtered" array
	for g=0, 2029 DO BEGIN
		for h=0, 1353 DO BEGIN
			if (filtered[h,g] gt thld) then filtered[h,g] = mask else filtered[h,g] = 0
			;note: cloud free regions may now be labelled as cloudy after the thresholding! Following line fixes that.
			if (CloudDataCloudFree[h,g] eq mask) then filtered[h,g] = mask


		ENDFOR
	ENDFOR

	CloudDataCloudFree = filtered
endif 




;open ALL the SDSs in the L1B file (EV_250_Aggr1km_RefSB, EV_500_Aggr1km_RefSB, EV_1KM_RefSB and EV_1KM_Emissive)
print, "Trying to open file "+l1b+" ..."

L1B_250m = HDFreturnSDS(l1b, 'EV_250_Aggr1km_RefSB')
num_250 = n_elements(L1B_250m)
print, "Read 250m data"
L1B_500m = HDFreturnSDS(l1b, 'EV_500_Aggr1km_RefSB')
num_500 = n_elements(L1B_500m)
print, "Read 500m data"
L1B_1km = HDFreturnSDS(l1b, 'EV_1KM_RefSB')
num_1km_vis = n_elements(L1B_1km)
print, "Read 1km vis data"
L1B_emiss = HDFreturnSDS(l1b, 'EV_1KM_Emissive')
; all MODIS L1B granules include an 'EV_1KM_Emissive' SDS so no need to consider the case of not existing. num_1km_TIR = n_elements(L1B_emiss)
print, "Read 1km TIR data"



;loop through all pixels
for j=0, 2029 DO BEGIN
	for k=0, 1353 DO BEGIN

		;MASK if bits 1 and 2 NOT equal to mask (if, mask if not clear)
		;OR MASK if not calculated (being conservative)

		if (CloudDataCloudFree[k,j] ne mask) OR (CloudDataCalculated[k,j] ne 1)then begin
			L1B_250m[k,j,*]=0
			L1B_500m[k,j,*]=0
			L1B_1km[k,j,*]=0
			L1B_emiss[k,j,*]=0
		endif

		;only allow pixels with no shadow contamination to pass
		;if CloudDataShadowFree[k,j] ne 4 then L1Bdata[k,j,*]=0
		; note - this appears to do nothin because the shadow mask is always 1.
		; see Ackerman (1997), JGR, for an explaination.


	ENDFOR
ENDFOR

print, "Masked all SDSs within HDF"

;write back to all SDSs in HDF, then close
if (num_250 gt 0) THEN BEGIN
	HDFeditSDS, l1b, 'EV_250_Aggr1km_RefSB', L1B_250m
	Print, "Wrote back 250m data"
endif
if (num_500 gt 0) THEN BEGIN
	HDFeditSDS, l1b, 'EV_500_Aggr1km_RefSB', L1B_500m
	Print, "Wrote back 500m data"
endif
if (num_1km_vis gt 0) THEN BEGIN
	HDFeditSDS, l1b, 'EV_1KM_RefSB', L1B_1km
	Print, "Wrote back 1km vis data"
endif
HDFeditSDS, l1b, 'EV_1KM_Emissive', L1B_emiss
Print, "Wrote back 1km TIR data"

print, "Masking completed!"

END