;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Takes a series of lat long coordinates and converts them
;;	to UTM Zone 14 coordinates with the builtin envi routine
;;	ENVI_CONVERT_PROJECTION_COORDINATES
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function convertLLtoUTM, locations
	s=size(locations)
	locations(*,1)=locations(*,1)*(-1)	;;to convert westings into eastings
	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



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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Gets envi header info (map and ns, nl, nb, type, interleave)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function getImageInfo, name
	 
     envi_open_file, name, r_fid=id
     if id eq -1 then return, -1
     
	envi_file_query, id, ns=ns, nl=nl, nb=nb, h_map=maph,	$
		interleave=interleave, data_type=type, descrip=desc

	handle_value, maph, mapinfo
	
	envi_file_mng, id=id, /remove
	
	return, {fileInfo, ns:ns, nl:nl, nb:nb, map:mapinfo,		$
			interleave:interleave, type:type, desc:desc}

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Sets envi header info (map and ns, nl, nb, type, interleave)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro setWindHdr, info, fname, ps, windSize
	info.map.ps = ps
	info.nb=1
	info.ns=windSize[1]
	info.nl=windSize[2]
	info.type=windSize[3] 
	
	envi_setup_head, fname=fname, ns=info.ns, nl=info.nl, nb=info.nb, $
		interleave=info.interleave, data_type=info.type, $
		descrip=info.desc, map_info=info.map, /write
end	



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;	Takes a wind meta file and krigs the data into a map
;;	based on the hdr file of a landsat scene
;;
;;	metawind is a meta file that contains, 
;;		column 1 windspeed (monthly mean max m/s)
;;		column 2-5 lat degrees, lat minutes, lon degrees, lon minutes
;;
;;	imgFile is the name of an ENVI header fileto get map info from it
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro spatwind, windFile, imgFile

	openw, oun, /get, string(windFile+'.out')
	cols=load_cols(windFile, winddata)
	loc=float(winddata(1:4, *))
	loc(0,*)=loc(0,*) + loc(1,*)/60
	loc(1,*)=loc(2,*) + loc(3,*)/60 
	print, loc
	locations=transpose(loc(0:1,*))
	npoints=n_elements(winddata)/cols

	utmLoc=convertLLtoUTM(locations)

	print, utmLoc

	bounds=headerBounds(imgFile)

	z=fltarr(npoints)
	e=[abs(bounds(2)-bounds(0))/2.,0.5] 

	flip=0		  ;false
	if bounds(1) gt bounds(3) then begin
		realbounds=bounds
		bounds(1) = bounds(3)
		bounds(3) = realbounds(1)
		flip=1	  ;true
	endif

	print, bounds
	gs=bounds(4:5)*10

;	for j=0, npoints-1 do begin
;		z(j)=winddata(1,j,where(centdata(0,j,*) eq centdats(i)))
;	endfor

	z=transpose(winddata(0,*))
	index =where(z ne 0 and z ne 0/0.)
	if n_elements(index) ne n_elements(z) then begin 
		print, 'Error with wind values, ', z
		z = z(index)
		utmLoc = utmLoc(*,index)
	endif
	
	print, 'Wind values = ',z
	print, 'Locations = ', utmLoc
	print, 'Range and Nugget = ', e
	res= krig2d(z, transpose(utmLoc(0,*)), transpose(utmLoc(1,*)), gs=gs, bounds=bounds(0:3), spherical=e)
	if flip then res=rotate(res,7)	;;mirrors the image about the x-axis

	print, size(res)
	writeu, oun, res
	imgInfo = getImageInfo(imgFile)
	setWindHdr, imgInfo, string(windFile+'.out'), gs, size(res)
	
	res=[0b,0b]

	close, oun
	free_lun, oun
end