;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Takes the name of an ENVI header file and finds the boundaries
;;	in map coordinates.  Returns upperleft, lowerright corners
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function headerBounds, headerfile, zone=zone
	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

        if keyword_set(zone) then zone=mapinfo.proj.params[0]

	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, zone
  FORWARD_FUNCTION envi_translate_projection_units, envi_proj_create
  
  if n_elements(zone) eq 0 then zone=14

	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=zone)

	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, zone

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

	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, zone
	
	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


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	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 stats, 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 = fltarr(nfiles)
	for i=0, nfiles-1 do begin
		cd, fnames(0,i),current=olddir
		cd, 'DATA'
		results(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, mean(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 krigged wind map will also be
;;		generated with the samename+'out'. This will be a floating point file
;;		with the same number of samples and lines as the input file, 
;;		co-registered with the image file.  
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro windstat, img, windbasedir, dataout
  envistart
	bounds=headerbounds(img, zone=zone)
	if bounds(0) eq -1 then return

;; 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, zone)
	cd, basedir

;; actual file names must include the year too
	fnames = getfilenames(fileIDs(0:2,*))
	
;; for all stations, find the mean monthly maximum windspeed, then print it
;;	and the location of that station to the output file	
	for i=0, n_elements(fnames(0,0,*))-1 do begin
		 monthlymeanmax = stats(fnames(*,*,i))
                 if (monthlymeanmax ne -1) and $
                   (monthlymeanmax gt 0 and monthlymeanmax lt 1000) then begin
                    printf, oun, monthlymeanmax, fileIDs(3,i), fileIDs(4,i), $
                      fileIDs(5,i), fileIDs(6,i)
                 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
	
	close, oun
	free_lun, oun
	
	cd, olddir
	
	spatwind, dataout, img

end