;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Finds the value with n percent of values above it and 100-(n percent) below it
;;		percentage is used to determine the cuttoff
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function n_percent, hist, percentage
	
	if percentage ge 1 then percentage=percentage/100.
	
	full = total(hist)
	max=0
	i=n_elements(hist)-1
	
	while max lt percentage do begin
		max=max+hist(i)/full
		i=i-1
	endwhile
	
	max=i+1
	return, max
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Takes the name of an ENVI header file and finds the boundaries
;;	in map coordinates.  Returns upperleft, lowerright corners
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function headerBounds, headerfile
	s=strsplit(headerfile, '.', /extract)
	name=headerfile
	if s(n_elements(s)-1) eq 'hdr' then begin
		name=s(0)
		if n_elements(s) gt 2 then $
			for i=1, n_elements(s)-2 do begin
				name=name+'.'+s(i)
			endfor
	endif
		
	envi_open_file, name, r_fid=fid
	if fid eq -1 then begin print, 'Bad Header File' & return, -1 & endif
	envi_file_query, fid, h_map=Hmapinfo, ns=ns, nl=nl, xstart=xstart, ystart=ystart
	if Hmapinfo eq -1 then begin print, 'Bad Header File' & return, -1 & endif
	Handle_Value, Hmapinfo, mapinfo

	pixsz=mapinfo.ps

;	left=(mapinfo.mc(0) + xstart )*pixsz(0)*(-1) + mapinfo.mc(2)
;	top=(mapinfo.mc(1) + ystart )*pixsz(1) + mapinfo.mc(3)
	left=mapinfo.mc(2)
	top=mapinfo.mc(3)
	right=left+(ns*pixsz(0))
	bottom=top-(nl*pixsz(1))

	envi_file_mng, /remove, id=fid

	return, [left, top, right, bottom, pixsz(0), pixsz(1)]
end



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Takes a series of lat long coordinates and converts them
;;	to UTM Zone 14 coordinates with the builtin envi routine
;;	ENVI_CONVERT_PROJECTION_COORDINATES, takes (nx2) array (lat,lon)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function convertLLtoUTM, locations
	s=size(locations)
	locations(*,1)=locations(*,1)*(-1)
	units=envi_translate_projection_units('Degrees')
	iproj=envi_proj_create(/geographic, units=units, datum='WGS-84')
	units=envi_translate_projection_units('Meters')
	oproj=envi_proj_create(/utm, units=units, datum='WGS-84', zone=14)

	envi_convert_projection_coordinates, locations(*,1), locations(*,0), iproj, $
						newXmap, newYmap, oproj
	map=fltarr(2, n_elements(newXmap))
	map(0,*)=newXmap
	map(1,*)=newYmap

	return, map
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	read the next line from the stnlist file and parse it into a station ID,
;;	a state ID and a lat deg, lat min, lon deg, lon min array
;;
;;	to get around the inability to store characters and numbers in the same array
;;	without making a structure, we will convert the state ID string into bytes first
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getnextline, unit
	s=''
	readf, unit, s
	s=strsplit(s, /extract)
	last = n_elements(s) -1

	stnid=0l
	reads, s(0), stnid
	
	stateid = byte(s(last-4))	;;two letters
	
	lat = strsplit(s(last-2), 'N', /extract)
	lon = strsplit(s(last-1), 'W', /extract)

	latn = 0
	lonn = 0
	reads, lat, latn
	reads, lon, lonn
	
;;	[stationID, stateID[2],  latitude,         	  longitude           ]
;;					       degrees, minutes      degrees, minutes
	return, [stnid, stateid, lat/100, lat mod 100, lon/100, lon mod 100]
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Determine whether or not the current position falls within the image boundary
;;	this would probably be much more efficient if it was performed on all points
;;	at once...
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;	
function IsItIn, curline, bounds

	lat = curline(3) + curline(4)/60.
	lon = curline(5) + curline(6)/60.
	
	utmloc = convertLLtoUTM([[lat],[lon]])

	return, utmloc(0) lt bounds(2) and utmloc(0) gt bounds(0) and $
			utmloc(1) lt bounds(1) and utmloc(1) gt bounds(3)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;	increase the area covered by the image boundary by factor meters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function enlargebounds, bounds, factor
	bounds(0) = bounds(0) - factor
	bounds(1) = bounds(1) + factor
	bounds(2) = bounds(2) + factor
	bounds(3) = bounds(3) - factor
	return, bounds
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	find all of the stations listed in the file stnlist that fall within 10km
;;	of the boundaries specified by bounds
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function pickfileIDs, bounds
	
	openr, sun, /get, 'STNLIST.TXT'
	
	bounds = enlargebounds(bounds, 100000)	;;increase boundary by 100Km on all sides
	fileIDlist = long([999999, 99,99,99,99,99,99]) ;;dummy first element
	
	while not eof(sun) do begin
		nextline = getnextline(sun)

;; determine if this station falls within the specified image boundaries, then add it
;; to the current file ID list if it does
		if IsItIn(nextline, bounds) then $
			fileIDlist = [[fileIDlist],[nextline]]
	endwhile
	close, sun
	free_lun, sun
	
;; remove the 0th element (the dummy initializer) and return the list of fileIDs
	return, fileIDlist[*, 1:n_elements(fileIDlist[0,*])-1]

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	finds all the wind files available for each station ID, max possible eq 16
;;	returns file paths relative to basedir.  this includes state 2-letter code
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getfilenames, fileIDs
	
;; again, this is an initial dummy file list
	filelist = strarr(2,16)			;;16 = 1982-1997
	filelist(*,*) = 'NA'

;; check for filenames of all stationIDs
	for i=0, n_elements(fileIDs(0,*))-1 do begin
		curbase = string(byte(fileIDs(1:2,i)))
		cd, curbase, current=olddir
		cd, 'DATA'
		tmpfile = strcompress(string(fileIDs(0,i))+'*', /remove_all)
		tmplist = findfile(tmpfile)

;;initialize the file list for this station
		curlist = strarr(2,16)
		curlist(*,*) = 'NA'
		
		if n_elements(tmplist) gt 1 then begin
			tmp=0
			
			for j=0, n_elements(tmplist)-1 do begin
				tmpname=tmplist(j)
				sep = strsplit(tmpname, '.', /extract)
;;if this file is not gzipped, add it to the file list
				if sep(n_elements(sep)-1) ne 'GZ' then begin
					curlist(0,tmp) = string(byte(fileIDs(1:2,i)))
					curlist(1,tmp) = tmpname
					tmp=tmp+1
				endif
			endfor
		endif
		cd, olddir

;;append the current list to the full list
		filelist = [[[filelist]],[[curlist]]]
	endfor


;; remove the first dummy element and return the list of files
	return, filelist[*,*,1:n_elements(filelist[1,1,*])-1]

end

function thisHistogram, fname
	openr, un, /get, fname
	s=''
	
	if eof(un) then begin
		close, un
		free_lun, un
		return, -1
	endif
	
	time=ulong64(0)
	speed = 0
	n=0
	tot=0.
	curlist=fltarr(10000)
	j=0l
	count = 0ul

;;	initialize month, and speed for this file
	readf, un, s
	s=strsplit(s,/extract)
	reads, s(1), time
	reads, s(3), speed
	if speed ne 9999 then begin
;		curmonth = time/1000000 ;;allows splitting by month
;		curmonth = curmonth mod 100
		curlist[0] = speed
		count=count+1
	endif
	
;; continue searching months until we reach the end of this file	
	while not eof(un) do begin
		s=''
		readf, un, s
		s=strsplit(s,/extract)
		if n_elements(s) gt 4 then begin
			reads, s(1), time
			reads, s(3), speed
			if speed ne 9999 then begin
;				curmonth = time/1000000 ;;allows splitting by month
;				curmonth = curmonth mod 100
				if count eq n_elements(curlist) then $
					curlist = [curlist,fltarr(10000)]
				
				curlist[count] = speed
				count=count+1
			endif
		endif
	endwhile
	close, un
	free_lun, un

	curlist = curlist(0:count-1)/10.

	return, histogram(curlist, binsize=0.5, min=0., max=100.)
end
	
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	open a single weather file and return the mean monthly max windspeed
;;	for that file
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function stat, fname
	openr, un, /get, fname
	s=''
	if eof(un) then begin
		close, un
		free_lun, un
		return, -1
	endif
	
	time=ulong64(0)
	speed = 0
	n=0
	tot=0.
	curmax=0
	j=0l

;;	initialize month, and speed for this file
	readf, un, s
	s=strsplit(s,/extract)
	reads, s(1), time
	reads, s(3), speed
	if speed eq 9999 then speed = 0 else j=1
	curmonth = time/1000000
	curmonth = curmonth mod 100
	curmax = max([curmax,speed])

;; continue searching months until we reach the end of this file	
	while not eof(un) do begin
		s=''
		readf, un, s
		s=strsplit(s,/extract)
		if n_elements(s) gt 4 then begin
			reads, s(1), time
			reads, s(3), speed
			if speed eq 9999 then speed = 0 else j=j+1; count missiing values as 0
			tmpmonth = time/1000000
			tmpmonth = tmpmonth mod 100
			
			if tmpmonth eq curmonth then $
				curmax = max([curmax,speed]) $
			else begin
				tot = tot+curmax
				if curmax eq 0 then begin
					print, fname,' has a max of ', curmax, ' for month', curmonth
					print, ' and there was at least', j, ' observation(s) with valid data'
				endif
				curmax = speed
				curmonth=tmpmonth
				n=n+(j gt 0)	;; n+1 if j gt 0, n otherwise
				if j lt 3 then print, 'Too few data points in ', fname, curmonth
				j=0
			endelse
		endif
	endwhile
	close, un
	free_lun, un

	if n gt 12 then print, n, ' months in ', fname
	return, tot/(10.*n)
end



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	takes a list of filenames (all of the same station) and outputs the mean
;;	monthly max windspeed in meters / second
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function findHistogram, fnames
	index = where(fnames(1,*) ne 'NA')
	if index(0) eq -1 then begin
		print, 'no valid files for this station ID'
		return, -1
	endif
	
	nfiles = n_elements(index)
	results = lonarr(202,nfiles)
	for i=0, nfiles-1 do begin
		cd, fnames(0,i),current=olddir
		cd, 'DATA'
		tmp=thisHistogram(fnames[1,index[i]])
		results(0:200,i) = tmp
		results(201,i) = stat(fnames[1,index[i]])
		cd, olddir
			
	endfor
	index = where(results eq -1)
	if index(0) ne -1 then $
		if (nfiles - n_elements(index)) lt 2 then return, -1 $
		else results = results[where(results ne -1)]
	return, results
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	assumes windbasdir is in the "standard"? format of the western US weather
;;	observations CD.  basedir contains 'DOCS' directory which contains
;;	'STNLIST.TXT'  and basedir also contains state directories by 2-lettercode
;;
;;	img must be a standard ENVI image file with .hdr map information
;;
;;	dataout is the name of the output file.  A postscript file showing the histograms
;;		off all wind speeds during the period of record for each station
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro windDisplay, img, windbasedir, dataout

	bounds=headerbounds(img)
	if bounds(0) eq -1 then return

	oldDevice=!D.name
	set_plot, 'ps'
	oldMulti=!p.multi
	device,file= dataout ,xsize=7.5, ysize=10, xoff=.5, yoff=.5, /inches, $
		set_font='Times', /tt_font, font_size=14
;; output file
;	openw, oun, /get, dataout
	
	cd, windbasedir, current=olddir

;;fileIDs = n x 5 array n files, [stnID,latD,latM,lonD,lonM]
	cd, 'DOCS', current=basedir
	fileIDs=pickfileIDs(bounds)
	cd, basedir

;; actual file names must include the year too
	fnames = getfilenames(fileIDs(0:2,*))

	nrows=10
	npages=0
	ncols=2
	while nrows gt 5 do begin 
		npages=npages+1
		nrows=n_elements(fnames(0,0,*))/(npages*ncols)
	endwhile
	!p.multi = [0,ncols,nrows]
	;print, fnames
	vals=fltarr(n_elements(fnames(0,0,*)),4)
;; for all stations, find the windspeed histogram, then print it
;;	and the location of that station to the output file	
	for i=0, n_elements(fnames(0,0,*))-1 do begin
		 hist = findHistogram(fnames(*,*,i))
		 if hist(0) ne -1 then begin
		 	tmp=total(hist(0:200,*), 2)
		 	index = where(tmp gt 0)
		 	mmm = mean(hist(201,*))
		 	
		 	if index(0) ne -1 then begin
				maxx = index(n_elements(index)-1)
				maxy = max(tmp)
				maxx = 60
				offset=maxy/10.
				twopercent=(n_percent(tmp, 2)/2.)
				onepercent=(n_percent(tmp, 1)/2.)
				vals(i,0) = mmm
				vals(i,1) = onepercent
				vals(i,2) = twopercent
				
				plot, tmp[0:index(n_elements(index)-1)], xr=[0,60],$
					title=string(fnames(0,0,i)+fnames(1,0,i))
				xyouts, 25, maxy-offset*2, 'Mon='+strmid(strtrim(mmm,1),0,5)+' m/s'
				xyouts, 25, maxy-offset*3, 'Max='+strmid(strtrim((index(n_elements(index)-1))/2.,1),0,5)+' m/s'
				xyouts, 25, maxy-offset*4, '2% ='+strmid(strtrim(twopercent,1),0,5)+' m/s'
				xyouts, 25, maxy-offset*5, '1% ='+strmid(strtrim(onepercent,1),0,5)+' m/s'
			endif
		endif else begin
			print, 'ERROR at : ',i
			print, fnames(*,i)
			print, fileIDs(i,*)
		endelse
		 print, byte(100*float(i)/(n_elements(fnames(0,0,*))-1)), ' % done'
	endfor
	cd, olddir
	
;	close, oun
;	free_lun, oun
	device, /close
	
	device,file= 'compare.ps' ,xsize=7.5, ysize=10, xoff=.5, yoff=.5, /inches, $
		set_font='Times', /tt_font, font_size=14
	!p.multi=[0,1,2]
	
	index = where(vals(*,0) ne 0 and vals(*,1) ne 0 and vals(*,2) ne 0)
	if index(0) ne -1 then begin
		plot, vals(index,0),vals(index,1), psym=3, xr=[0,40], $
			title='Mean Monthly Max vs. Top One Percent'
		plot, vals(index,0),vals(index,2), psym=3, xr=[0,40], $
			title='Mean Monthly Max vs. Top Two Percent'
	endif
	device, /close
	
	openw, oun, /get, 'mmmvspercent.data'
	writeu, oun, long(n_elements(index))
	writeu, oun, vals(index,0), vals(index,1), vals(index,2)
	close, oun
	free_lun, oun
	!p.multi=oldMulti
	set_plot, oldDevice
end