; This routine reads a 2B-GEOPROF file specified by infile and writes
; png images that match the quicklook segments found at
; http://www.cloudsat.cira.colostate.edu/dpcstatusQL.php
;
; Make sure you have a subdirectory called "images" under the data_directory.
; Also, to get a good range of scales for radar reflectivity the min and max
; values were set to -40 and 50.  Just look for "MINMAX" in the code below
; to change it.
;
; Phil Partain (partain@cira.colostate.edu)
; CloudSat Data Processing Center
; Cooperative Institute for Research in the Atmosphere
; Colorado State University
; Date: 9/27/06

pro geoprof_png

	infile = '2006190054715_01047_CS_2B-GEOPROF_GRANULE_P_R02_E01.hdf'
	data_directory = 'C:\geoprof_test\'
	swathname = '2B-GEOPROF'

	fieldname1 = 'Radar_Reflectivity'
	missing1 = -88.88
	factor1 = 100
	offset1 = 0
	missop1 = 'le'

	;fieldname2 = 'CPR_Cloud_mask'
	;missing2 = 247
	;factor2 = 1
	;offset2 = 0
	;missop2 = 'eq'

	print,infile

	prefix = strmid(infile,0,23)

	;--------------------------------------------------------------------

	fid = eos_sw_open(data_directory+infile,/READ)
	print,'file',fid
	swathid = EOS_SW_ATTACH(fid, swathname)
	print,'swath',swathid
	status = EOS_SW_READFIELD(swathid, 'Profile_time', time)
	print,'status',status
	status = EOS_SW_READFIELD(swathid, fieldname1, field1)
	print,'status',status
	;status = EOS_SW_READFIELD(swathid, fieldname2, field2)
	;print,'status',status
	status = EOS_SW_DETACH(swathid)
	print,'status',status
	status = EOS_SW_CLOSE(fid)
	print,'status',status

	plotfield,data_directory,prefix,swathname,field1,fieldname1,missing1,missop1,factor1,offset1
	;plotfield,data_directory,prefix,swathname,field2,fieldname2,missing2,missop2,factor2,offset2

end

;========================================================================
pro plotfield,data_directory,prefix,swathname,field,fieldname,missing,missop,factor,offset

	LOADCT, 39

	; set the missing data to a color that doesn't really show up in color table 39: purple
	TVLCT, red, green, blue, /GET
	red(0) = 200
	green(0) = 0
	blue(0) = 200
	tvlct,red,green,blue

	; scale the field
	field = (float(field)-offset)/factor

	; set up the window size
	win_xsize = n_elements(field(0,*))
	win_ysize = n_elements(field(*,0))

	; find out where the good and missing data are
	CASE missop OF
	    'eq': whgood = where(field NE missing)
	    'le': whgood = where(field GT missing)
	    'ge': whgood = where(field LT missing)
	    'lt': whgood = where(field GE missing)
	    'gt': whgood = where(field LE missing)
	ENDCASE

	CASE missop OF
	    'eq': whbad = where(field eq missing)
	    'le': whbad = where(field le missing)
	    'ge': whbad = where(field ge missing)
	    'lt': whbad = where(field lt missing)
	    'gt': whbad = where(field gt missing)
	ENDCASE

	; MINMAX
	;mx = max(field(whgood))
	mx = 50
	;mn = min(field(whgood))
	mn = -40

	print,"min field",mn
	print,"max field",mx

	; scale the image
	image = BYTSCL(field,MIN=mn,MAX=mx,top=254) + 1

	; set the missing data
	print,'n missing:',n_elements(whbad)
	if (n_elements(whbad) gt 1) then image(whbad) = 0

	image = reverse(transpose(image),2)

	;--------------------------------------------------------------------
	; set up the segments to match the standard cloudsat quicklooks
	NumSegments = ceil(float(win_xsize)/1200.)
	SegmentInc = win_xsize/NumSegments

	IndexBeg = 0L
	IndexEnd = 0L

	; loop over the segments
	for i = 0, NumSegments-1 do begin
	    IndexBeg = IndexEnd
	    IndexEnd = IndexEnd + SegmentInc

	    IntIndexBeg = long(IndexBeg)
	    IntIndexEnd = long(IndexEnd)

	    ;print,IntIndexBeg,IntIndexEnd-1

		; satellite travels from right to left
	    smallimg = reverse(image(IntIndexBeg:IntIndexEnd-1,*))

	    ; make the image
	    s = SIZE(smallimg)
	    image3d = BYTARR(3, s(1), s(2))
	    image3d(0, *, *) = red(smallimg)
	    image3d(1, *, *) = green(smallimg)
	    image3d(2, *, *) = blue(smallimg)

		; set up the filename
	    filename = data_directory+'images\'+prefix+swathname+'_'+fieldname+'_'+strtrim(i+1,2)+'.png'

		; write the png file
		WRITE_PNG, filename, image3d, r,g,b

		; want jpg?
		;filename = data_directory+'images\'+prefix+swathname+'_'+fieldname+'_'+strtrim(i+1,2)+'.jpg'
	    ;WRITE_JPEG, filename, image3d, TRUE=1

	endfor

	print, "Finished!"

end

;==============================================================================