;+
; NAME:
;	MOD35Display.pro
;
; PURPOSE:
;	To display a bit- and spatial-subset of an HDF file (useful for interpreting MODIS cloud mask individual spectral tests)
;	Note: look at the console output to see what each greyscale value might mean. Values shown here are scaled relative to the start of each byte, so you will have to convert in your head
;
; AUTHOR:
;	Alex Fraser
;
; CATEGORY:
;	Cloud Mask Utility
;
;
; CALLING SEQUENCE:
;	see example
;
;
; INPUTS:
;	cloudMaskFile: The location and filename of the cloud mask HDF file (will prompt if not specified)
;	bitStart: The starting bit number of the data to be displayed
;	numBits: The number of bits following the start bit to be scaled and displayed. Defaults to 1.
;	xrange: a 2-element vector describing the x-range of the output plot
;	yrange: a 2-element vector describing the y-range of the output plot
;
; KEYWORDS:
;	/confidence: If set, will use the "Quality_Assurance" SDS instead of the "Cloud_Mask" one. 
;	Google "Strabala Cloud Mask User's Guide" for more info on what each SDS contains.
;
; OUTPUTS:
;	Output to display
;
; RESTRICTIONS:
;	Heaps! 
;	Note: Required a compiled version of HDFreturnSDS.
;
; EXAMPLE:
;	This one will show the results of the simple IR threshold test.
;	MOD35Display, cloudMaskFile='c:\whatever.hdf', bitStart=13, xrange=[1173, 1223], yrange=[1354, 1404]
;
;	this one will display the 2 bits (4 greyscale levels including black and white) showing the number of spectral tests used to generate the cloud mask
;	MOD35Display, cloudMaskFile='c:\whatever.hdf', bitStart=50, numBits=2, xrange=[1173,1223], yrange=[1354,1404], /conf
;
; MODIFICATION HISTORY:
;
;	Written by Alex Fraser, December 2007.
;-

PRO MOD35Display, cloudMaskFile=CldMskFile, bitStart=bF, numBits=nB, xrange=xrng, yrange=yrng, confidenceSDS=confidence

if n_elements(CldMskFile) eq 0 then begin
	print, "Pick the MODIS cloud mask granule of interest"
	CldMskFile = dialog_pickfile(PATH='/u/adfraser/IDL/HDFtesting/HDFdisplay/')
endif

print, CldMskFile

if n_elements(confidence) eq 0 then begin
	;CloudData is a 3D array of size [1354, 2030, 6]
	CloudData = HDFreturnSDS(CldMskFile, 'Cloud_Mask')
endif else begin
	;CloudData is a 3D array of size [1354, 2030, 10]
	CloudData = HDFreturnSDS(CldMskFile, 'Quality_Assurance')
endelse

if n_elements(nB) eq 0 then begin
	;set numBits to 1 if not defined.
	nB = 1
endif

;break up the bitField into a byte and a bit offset
byte=bF/8
print, byte

;assemble the bit mask
;bitmask=2^(bF MOD 8)
bFstart = bF MOD 8
bitmask=0
for j=0, nB-1 DO BEGIN
	bitmask = bitmask + 2^(bFstart+j)
ENDFOR
print, bitmask

;extract the relevant byte of data from CloudData
if n_elements(confidence) eq 0 then begin
	data = CloudData[*,*,byte]
endif else begin
	data = CloudData[byte,*,*]
	data = REFORM(data)
endelse

;mask that array
data = data AND bitmask

;pull out the subset we want
data = data[xrng[0]:xrng[1],yrng[0]:yrng[1]]

;reverse the row order to correct for IDL drawing it upside down
numrows = yrng[1]-yrng[0]
numcols = xrng[1]-xrng[0]

data2 = data

for j=0, numrows-1 DO BEGIN
	data2[*,j]=data[*,numrows-j]
ENDFOR

;tvscl handles the normalization for us, so no need to do that.
window, 0, xsize=(xrng[1]-xrng[0]), ysize=(yrng[1]-yrng[0]) & tvscl, data2

print, "Max of this array:"
print, max(data2)
print, "min of this array:"
print, min(data2)




END