; IDL_read.pro
;
; DISCLAIMER: This software is provided for users of TRMM data.
; TSDIS DOES NOT PROVIDE SUPPORT FOR THIS SOFTWARE.
;
; Sample code for reading some TRMM 2A25 fields into interactive IDL
; Minimal error checking, be sure you are really reading a TRMM 
; 2A25 HDF file.
;
; Some of this source is hard-coded for 2A25 for the sake of simplicity.
; The following fields are read from the file and are accessable from the
; IDL prompt after execution:
;
;    SDSs:
;    geolocation(nray,nscan)
;    rain(ncell,nray,nscan) [nbin=80,nray=49,nscan= ~9150] This array is BIG.
;                           This sample only reads in all bins for all rays
;                           in first 10 scans.
; 
;    Vdata:
;    PR_scan_status(nscan) -> Fractional Orbit Number
; 
;
; All the above fields are documented in the TSDIS file specifications.
; READ THE FILE SPECIFICATIONS!  They are described in the ICS Vols 3,4,6
; on the world wide web at http://tsdis.gsfc.nasa.gov/tsdis/tsdis.html
;
; Also see the TSDIS Toolkit include files for the variable names
;
; No warranty expressed or implied. Some code copied from an RSI
; IDL/HDF example
; 
; The Vdata code has been taken and modified from the TSDIS OrbitViewer.
;
; Requirements: X-windows and lots of memory
;   Uses an X dialog box for filename input.  If your not using X alter 
;   the code and pass the filename as an argument.
;   This code also displays a sample plot of geolocation data. Comment
;   that out in the main routine if your not using X.
;
; Usage:
;  
; IDL> .run IDL_read_HDF
;
;
;------------------------------------------------------------------------------


;------------------------------------------------------------------------------
; Procedure to read the SDS and Vdata fields

  pro  read2a25, geolocation, rain, frac_orbit_num, file_only


    print, '-------- Select a 2A25 granule  --------'
; Use built-in IDL to pop up a filename dialog box 
    filen = dialog_pickfile()
   

; Do some error checking on the filename
; First 4 chars should be '2A25'
; Call remove_path routine which is included below
    remove_path, filen,file_only
    first_four = ''
    reads, file_only,first_four, format='(A4)' 
    if (first_four ne '2A25')then begin
	print, 'Are you sure ',file_only,' is a 2A25 HDF product?'
        print, 'If so, first 4 chars should be 2A25'
        retall
        stop
    endif

    print, 'This may take some time if reading a full 250MB granule'
    print,''

; Initialize the scientific data set interface (hdf_sd* routines) 
    sdsfileid = hdf_sd_start(filen,/read)

; Get info on the SDSs, number of SDSs and number of global attributes 
    hdf_sd_fileinfo,sdsfileid,numsds,ngatt
    names = strarr(numsds)
    ndims = lonarr(numsds)
    dtype = strarr(numsds)

; Print out a table of the name, number of dimensions and type of each SDS
; in file    
    for i = 0, numsds - 1 do begin
	sds_id = hdf_sd_select(sdsfileid, i)
        hdf_sd_getinfo, sds_id, name = na, ndim = nd,type= typ
        names( i ) = na 
        ndims( i ) = nd
        dtype(i) = typ  
    endfor
    F1='(" ",A,I4)'
    print,'List of SDS names'
    print,' # of SDSs = ',numsds,FORMAT=F1
    print,''
    if numsds gt 0 then begin
        print,"     Label       Dims   Type"
        print,"---------------- ---- --------"
        for i=0,numsds-1 do begin
	    print,names(i),ndims(i),dtype(i),FORMAT='(A14,"   ",I4," ",A8," ")'
        endfor
	print,"---------------- ---- --------"
    endif

; Read the geolocation SDS. This is kind of cheating since I know ahead of
; time that the name of the SDS is 'geolocation'.  However, I could run
; the code above to find out the name of the SDS containing geolocation.

; Get ID of SDS 
; This line calls another built in HDF function inside hdf_sd_select to 
; translate the name of the SDS into an index.
    sds_id = hdf_sd_select(sdsfileid,hdf_sd_nametoindex(sdsfileid, 'geolocation'))

; Read the SDS data into a variable. This reads the entire geolocation array
    hdf_sd_getdata, sds_id, geolocation

; Read 3D rainfall array 
    sds_id = hdf_sd_select(sdsfileid,hdf_sd_nametoindex(sdsfileid, 'rain'))
    hdf_sd_getdata, sds_id, rain,count=[80,49,10]
; If count is omitted above, the entire SDS will be read in. This line
; reads in all 80 vertical bins for all rays in the first 10 scans.
; For some reason using the count keyword also makes a temporary hidden
; variable that generates a message when returning from the procedure.
; Take the count keyword off and you won't get the message. No harm done.

; We are done with SDSs so we can close the interface
    hdf_sd_end, sdsfileid

; Read in the Scan Status Vdata
; Scan status contains various flags for each scan. The fields are:
;
; Name                   Size Type Position
; ----------------------------------------
; Missing                1byte int   0 
; Validity               1byte int   1
; QAC                    1byte int   2
; Geolocation Quality    1byte int   3
; Data Quality           1byte int   4
; Current SC Orientation 1byte int   5
; Current ACS mode       1byte int   6
; Yaw Update Status      1byte int   7
; PR mode                1byte int   8
; PR Status 1            1byte int   9
; PR Status 2            1byte int   10
; Fractional Orbit Num.  4byte float 11
;
; How did I know this?  READ THE FILE SPECIFICATIONS
;
 
; Open the file for read and initialize the Vdata interface
    file_handle = hdf_open(filen,/read)

; Get the ID of the first Vdata in the file 
    vdata_ID = hdf_vd_getid(  file_handle, -1 )
    is_NOT_fakeDim = 1
    num_vdata = 0

; Loop over Vdata 
    while (vdata_ID ne -1) and (is_NOT_fakeDim) do begin

; Attach to the vdata 
	vdata_H = hdf_vd_attach(file_handle,vdata_ID)

; Get vdata name
	hdf_vd_get, vdata_H, name=name,size= size, nfields = nfields
        
; Detach vdata
	hdf_vd_detach, vdata_H
       
; Check to see if this is a dummy
; Can't really explain why this happens but sometimes a dummy dimension
; gets returned as a Vdata name depending on the HDF file.  
	is_NOT_fakeDim = strpos(name,'fakeDim') eq -1

; Build up the list of Vdata names,sizes and number of fields 
	if (num_vdata eq 0) then begin
	    Vdata_name = name 
	    Vdata_size = size
	    Vdata_nfields = nfields
	    num_vdata = 1
	endif else if is_NOT_fakeDim then begin
	    Vdata_name = [Vdata_name,name]
            Vdata_size = [Vdata_size,size]
            Vdata_nfields = [Vdata_nfields,nfields]
	    num_vdata = num_vdata + 1
	endif
       
; Get ID of next Vdata
	vdata_ID = hdf_vd_getid( file_handle, vdata_ID )

    endwhile 

; Print out the list of names   
    print,''   
    print, 'List of Vdata names    Size (bytes)   Num. Fields'
    print, '-------------------------------------------------'
    for i = 0,num_vdata-1  do begin
	print, Vdata_name(i),Vdata_size(i),Vdata_nfields(i),$
	       format='(A18,I10,I14)'
    endfor
    print, '-------------------------------------------------'

; Find the Scan status Vdata
    vdata_ID = hdf_vd_find(file_handle,'pr_scan_status')

; Attach to this Vdata
    vdata_H = hdf_vd_attach(file_handle,vdata_ID)

; Get the Vdata stats
    hdf_vd_get,vdata_H,name=name,fields=raw_field

; Separate the fields
    fields = str_sep(raw_field,',')

; Read the Vdata, returns the number of records
; The data for all records is returned in a BYTE ARRAY of (record_size,nscans)
; IDL will issue a warning to remind you there are mixed data types in
; the array
    nscan = hdf_vd_read(vdata_h,data)
; Could have just read in the fractional orbit number with the
; fields keyword but this shows you how to extract the data from the 
; full record BYTE array.

; Make up an array for the fractional orbit number
    frac_orbit_num = fltarr(nscan)

; Loop over the records and pull out the fractional orbit number   
    for i = 0,nscan-1 do begin

; We know that the frac_orbit_number starts at position 11 in the byte array
      frac_orbit_num(i) = float(data(*,i),11)
    endfor

; Detach from the Vdata
    hdf_vd_detach, Vdata_H
   
; Close the hdf file
    hdf_close,file_handle   
    return
    end
;------------------------------------------------------------------------------

;------------------------------------------------------------------------------
; Procedure to remove path from file name
; Taken from OrbitViewer by Owen Kelley
       
   pro remove_path, file, file_no_path
;------------------------------------
;
;------------------------------------
        file_no_path    = file

        while( strpos( file_no_path , '/' ) ne -1 ) do $
          file_no_path = strmid( file_no_path , $
            strpos( file_no_path , '/' ) + 1, 1000 )

   end


;------------------------------------------------------------------------------
; Main routine
;
; All this does is call the read2a25 procedure and print out some values
; The main routine returns you to $MAIN$ with the three 
; variables defined.  Use 'help' to see them.
;

    read2A25, geolocation,rain,frac_orbit_num,filen
    
    print,'All data has been read, you can ignore the above warning'
    print,''
    print,'Hit Enter to output data to screen'
    dummy = get_kbrd(1)
    
; UNPACK THE RAIN VARIABLE
; The numbers in the file are mm/hr*10 
; The missing flag is not unscaled
  rain_out=fltarr(80,49,10)
  for i=0,79 do begin
    for j=0,48 do begin
      for k=0,9 do begin
        case rain(i,j,k) of
           -9999:  rain_out(i,j,k) = -9999.0
           else:   rain_out(i,j,k) = rain(i,j,k)/10.0
        endcase
      endfor
    endfor
  endfor

; Print out some of the values
    print,''
    print, '-------------- OUTPUT for file ',filen
    print,''
; Rain rate profile for first ray in first scan
    print, 'Rainfall rate in mm/hr, profile for first scan, first ray'
    print, 'Index 79 is the Earth ellipsoid, bins have 250meter resolution'
    print, 'Index  Rainrate      Index  Rainrate      Index  Rainrate '
    print, '------------------------------------------------------------'
    for i = 0,79 do begin
    	print, i, rain_out(i,0,0),format='(I4,2X,F7.2,TR8,$)'
        if ((((i+1) mod 3) eq 0) and (i ne 0 )) then print,format='(/,$)'  
    endfor
    print,''
    print,''

; Fractional orbit number for first 10 scans
    print, 'Fractional Orbit Number, first 10 scans'
    print, 'Scan Number    Fractional Orbit Number'
    print, '--------------------------------------'
    for i = 0,9 do begin
         print,i, frac_orbit_num(i),format='(I5,F25.8)'
    endfor

; Plot the latitude for the nadir ray (24)
; Comment this out if you are not connected to an X device
    plot, geolocation(0,24,*),/ynozero,title='Nadir Ray Latitude',$
          xtitle = 'Scan Number', ytitle = 'Latitude in degrees'
    print, ''
    print,'Type "help" to get a list of variables and sizes'

  end