;+
; NAME:
;	wso_read
; PURPOSE:
;	Read file with Stanford synoptic map data into 2D array
; CATEGORY:
;	gen/idl/util
; CALLING SEQUENCE:
	FUNCTION wso_read, wsofile, Y,  $
		_extra		= _extra	, $
		rotation	= rotation	, $
		longitude	= longitude , $
		cvariable	= cvariable , $
		align		= align 	, $
		header		= header	, $
		todegrees	= todegrees , $
		padto90 	= padto90	, $
		dropfirst	= dropfirst , $
		droplast	= droplast	, $
		fill		= fill		, $
		badfillvalue= badfillvalue, $
		errormessage= errormessage
; INPUTS:
;	wsofile			scalar; type: string
;						name of WSO data file
; OPTIONAL INPUTS:
;	align=align 	scalar; type: integer; default: not set
;						longitude at start of array.
;						Used to rearrange the first (longitudinal) dimension of
;						the array to put longitude 'align' at the start of the
;						array.
;	todegrees=todegrees
;					if 30 numbers per longitude are read the grid is assumed to
;					have equal steps of sine(latitude) from -14.5/15 to +14.5/15
;					If /todegrees is set this is converted to 29 equal steps in
;					latitude from -70 to +70 degrees in steps of 5 degrees
;					(using spline interpolation).
;	padto90=padto90
;					If 29 numbers per longitude are read (or 30 with /todegrees
;					set) the grid is assumed to cover -70 to +70 degrees in
;					steps of 5 degrees. If /padto90 is set then 4 rows of bad
;					values are added to the top and bottom of the array so that
;					the grid covers -90 to +90 degrees.
;
;	/dropfirst		remove the first longitude (as stored in input file)
;	/droplast		remove the last longitude (as stored in input file)
;						the output array will have only 72 longitudes.
;	/fill			fill bad values read from file by a call to gridfill or set
;						them to badfillvalue (if set).
;						This only applied to bad values read from the data file;
;						NOT for the padding to +/- 90 degrees if /padto90 is set.
;	badfillvalue=badfillvalue
;					scalar; type: float; default: none
;						used to fill bad values (NaNs) if keyword /fill is set
; OUTPUTS:
;	status			scalar; type: integer
;						return status (0=Failure,1=Success)
;	map				2D-array[N,M]; type: float
;						map=-1 if status=0
;						2D magnetic field map (in units stored on file)
;						N=73, M=30: equal steps in sine(latitude) from -14.5/15 to +14.5/15
;						N=73, M=29: equal steps in latitude from -70 to +70 degrees
;						N=73, M=37: equal steps in latitude from -90 to +90 degrees
;									(with 4 rows of bad data at top and bottom)
;						The first dimension stores data in order of increasing
;						heliographic longitude (i.e. later data first).
;						The second dimension stores data in order of increasing
;						latitude.
; OPTIONAL OUTPUTS:
;	rotation=rotation
;					array[N]; type: integer
;						Carrington rotation for each column in Z
;	longitude=longitude
;					array[N]; type: integer
;						heliographic longitude for each column in Z
;	cvariable=cvariable
;					array[N]; type: double
;						Carrington variable = rotation+(1-longitude/360)
;						If keyword align is used then this array is NOT in
;						monotonic decreasing order anymore !!!
;	header=header	scalar or array; type: string
;						records preceeding the magnetic data.
;						(May not exist on return).
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; RESTRICTIONS:
;	The most recent synoptic map on the WSO map sites could be incomplete.
;	In that case the returned array Z will contain less than 73 longitudes
; CALLS:
;	InitVar, gunzip_file, flt_read, BadValue, gridgen, IsType, boost, gridfill
;	do_file, TimeSet, TimeGet
; PROCEDURE:
;	The output array covers 360 degrees in longitude (not necessarily starting
;	at 0 degrees) and -90 to +90 in latitude in steps of 5 degrees.
;	The data only extend up to +/-70 degrees latitude, i.e the top and bottom
;	four rows are filled with bad values.
;
;	The input map maps may be in equal steps of sine(latitude). These are
;	converted to equal steps in latitude by cubic spline interpolation.
; MODIFICATION HISTORY:
;	NOV-2002, Paul Hick (UCSD/CASS)
;	JAN-2003, Paul Hick (UCSD/CASS)
;		Added /fill and bad fillvalue keywords to fill bad values (NaNs).
;		These were set to zero before.
;	OCT-2003, Paul Hick (UCSD/CASS)
;		Added check for bad values of -9999.90 in the WSO files.
;		The NOAA files occasionally use this value to fill missing data. These
;		are set to BadValue() in the output array by flt_read.
;	MAY-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		If keyword 'OBSTIME' is not found then the time assigned to the file
;		is extracted from keyword 'UTLAST' (in Feb 2006 the Fits header were
;		changed and OBSTIME is now store in UTLAST, it seems.
;-

InitVar, todegrees, /key
InitVar, padto90  , /key
InitVar, dropfirst, /key
InitVar, droplast , /key
InitVar, fill	  , /key

; This should always produce an 8 x 292 array.

type = GetFileSpec(wsofile,part='type')
IF type EQ '.gz' THEN type = GetFileSpec(GetFileSpec(part='name'),part='type')

SWITCH type OF

'.dat': BEGIN		; WSO ascii file

	status = flt_read(wsofile, Y, atleast=7, errormessage=errormessage, _extra=_extra,	$
		padvalue=BadValue(0.0), comments=comments, bad=[-9999.90,9999.90])

	IF status THEN BEGIN

		; Comments contains lines from the start of the file with less
		; than 7 numbers. These should be the lines preceeding the first
		; record starting with CT.

		IF IsType(comments,/defined) THEN BEGIN
			n = where(comments NE '')
			IF n[0] NE -1 THEN header = comments[n]
		ENDIF

		; Rearrange array from 8 x 292 to 32 x 73

		n = size(Y,/dim)
		Y = reform(Y,n[0]*4,n[1]/4, /overwrite)
		n = size(Y,/dim)

		; For files storing 29 latitudes in equal steps in degrees the last
		; column now contains bad values. Check and remove column, if present.

		IF total( finite(Y[n[0]-1,*]) ) EQ 0 THEN BEGIN
			Y = Y[0:n[0]-2,*]
			todegrees = 0				; Latitudinal grid already in degrees
		ENDIF

		; Extract Carrington rotation (1st column) and longitude (2nd column)
		; and construct Carrington variable

		cvariable = reform( Y[0,*]+(360d0-Y[1,*])/360d0 )

		; Drop first 2 columns and rearrange from (29 or 30) x 73 to 73 x (29 or 30)

		Y = Y[2:*,*]					; Extract magnetic field data
		Y = transpose(Y)				; Longitude now is first dimension

	ENDIF

	break

END

'.fits':
'.fts' : BEGIN

	message, /info, hide_env(wsofile)

	Y = readfits(wsofile, hdr)
	status = IsType(Y,/array)

	IF status THEN BEGIN

		; n[0] is always 72?

		n = size(Y,/dim)
		even = n[0]/2*2 EQ n[0]

		IF even THEN Y = [Y,Y[0,*]]

		n = size(Y,/dim)
		even = n[0]/2*2 EQ n[0]

		; Leading Carrington rotation and heliographic longitude
		; Combine into leading Carrington variable

		carrot = fxpar(hdr,'CRROTL')
		IF !err NE 0 THEN carrot = fxpar(hdr,'CARROT')
		IF !err NE 0 THEN message, 'could not determine leading Carrington rotation, '+hide_env(wso_file)

		carrlong = fxpar(hdr,'SNLONGL')
		IF !err NE 0 THEN carrlong = fxpar(hdr,'CARRLONG')
		IF !err NE 0 THEN message, 'could not determine leading Carringto longitude, '+hide_env(wso_file)

		cvariable = carrot+(360.0d0-carrlong)/360.0d0

		even = n[0]/2*2 EQ n[0]
		cvariable = gridgen(n[0]+even, range=cvariable-[0,1])
		cvariable = cvariable[0:n[0]-1]

		; The Fits file store the map with increasing longitude in the first
		; and increasing latitude in the second dimension.
		; The ASCII files do the exact reverse. Rearrange here to become consistent
		; with the ASCII files

		cvariable = reverse(cvariable)
		Y = reverse(reverse(Y,1),2)

		header = fxpar(hdr,'OBSTIME')
		IF !err GE 0 THEN BEGIN		; If keyword OBSTIME present
			; For at least one wso_noaa file the keyword OBSTIME is specified as
			;	OBSTIME = 'Observation time: 2005:06:28_22h:59m:45s'
			; Strip 'Observation time: ' if present
			IF strpos(header,'Observation time: ') EQ 0 THEN BEGIN
				message, /info, 'fix: '+header
				header = strmid(header,strlen('Observation time: '))
			ENDIF
		ENDIF ELSE BEGIN		; If keyword UTLAST is present (starting Feb 2006)
			header = fxpar(hdr,'UTLAST')
		ENDELSE

		CASE strlen(header) EQ strlen('YYYY:MN:DD_hhH:mmM:ssS') OF
		0: BEGIN
			header = TimeSet(header)
			header = ['OBSERVAT='+fxpar(hdr,'OBSERVAT'),'OBSTIME='+TimeGet(header,/_ydoy,upto=TimeUnit(/sec))]
		END
		1: BEGIN
			header = TimeSet(header,format='YYYY:MN:DD_hhH:mmM:ssS')
			header = ['OBSERVAT=WSO-NOAA','OBSTIME='+TimeGet(header,/_ydoy,upto=TimeUnit(/sec))]
		END
		ENDCASE

	ENDIF

	break

END

ENDSWITCH


IF status THEN BEGIN

	IF IsType(align, /defined) THEN BEGIN

		longitude = round((1-(cvariable-floor(cvariable)))*360) mod 360

		; Rearrange the first (longitudinal) dimension to put longitude 'align'
		; in the first column.

		tmp = min( abs(longitude-align), m )
		message, /info, 'aligning at'+strcompress(longitude[m])+' degrees'

		n = size(Y, /dim)
		n = indgen(n[0])

		CASE m OF
		0   : 
		1   : n = [ n[m:*], n[m] ]
		ELSE: n = [ n[m:*], n[1:m-1], n[m] ]
		ENDCASE

		; The Carrington variables are not monotonic anymore !!!

		Y = Y[n,*]
		cvariable = cvariable[n]

	ENDIF

	CASE 1 OF
	dropfirst: BEGIN 					; Drop first record in file = earliest data
		Y = Y[1:*,*]
		cvariable = cvariable[1:*]
	END
	droplast: BEGIN						; Drop last record in file = most recent data
		n = (size(Y, /dim))[0]
		Y = Y[0:n-2,*]
		cvariable = cvariable[0:n-2]
	END
	ELSE:								; Nothing to do
	ENDCASE

	; The file stores data from high to low latitude in each row and high
	; to low longitude in each column. So both dimensions need to be reversed.

	cvariable = reverse(cvariable)
	Y = reverse(reverse(Y,1),2)

	rotation  = floor(cvariable)		; Carrington rotation
	longitude = round((1-(cvariable-rotation))*360)
										; Heliographic longitude (degrees)
	; Check for bad values and fill them if requested
	; (XuePus program sometimes produces isolated bad values).

	IF fill THEN BEGIN
		bad = where(1-finite(Y),n)
		IF n NE 0 THEN BEGIN
			message, /info, strcompress(n,/rem)+' bad values filled'
			CASE IsType(badfillvalue, /defined) OF
			0: Y = gridfill(Y)
			1: Y[bad] = badfillvalue
			ENDCASE
		ENDIF
	ENDIF

	IF todegrees THEN BEGIN

		message, /info, 'converting to equal steps in latitude'

		n = size(Y,/dim)

		; Convert to equal steps in latitude by cubic spline interpolation

		;x0 = gridgen(n[1], range=14.5/15*[-1,1])
		x0 = gridgen(n[1],range=[-1,1],/open)
		x1 = sin( gridgen(n[1]-1, range=70/!radeg*[-1,1]) )

		R = make_array(type=IsType(Y), dim=[n[0],n[1]-1])
		for i=0,n[0]-1 do R[i,*] = spl_interp(x0,Y[i,*],spl_init(x0,Y[i,*]),x1)

		Y = R

	ENDIF


	IF padto90 THEN BEGIN

		message, /info, 'padding to +/- 90 degrees latitude'

		n = size(Y,/dim)

		; Add couple of rows at top and bottom of map to fill map up to +/- 90 degrees

		n = (37-n[1])/2 		; Should be 4
		Bad = BadValue(Y)
		boost, Y, pushdim=1, plus=n, /lead, value=Bad
		boost, Y, pushdim=1, plus=n, /top , value=Bad

	ENDIF

ENDIF

RETURN, status  &  END