; =============================================================
;
; Example of reading HDF 4.0 data files of the UV index or the UV dose
; as produced by TEMIS; for information, see
; http://www.temis.nl/uvradiation/info/uvhdfread.html
;
; (c) TEMIS/KNMI, last modified: 11 February 2004
;
; =============================================================
;
; Usage within IDL to compile
;
;    IDL> .run uvhdfread
;
; And to run:
;
;    IDL> uvhdfread
;
; ====================================================================

pro uvhdfread

; Read the name of the HDF file to be read and remove
; leading and trailing spaces
    print,''
    print,'>>> Give the name of the HDF file to read'
    filename=''
    read,filename
    filename=strtrim(filename,2)
    print,''


; Open the HDF file for reading; the 'sd_id' is an identifier
; linked to the file and used below when reading from the file.
    sd_id=HDF_SD_START(filename,/READ)
    if sd_id le 1 then begin
       print,' *** Error opening HDF file'
       goto,endprog
    endif


; To retrieve the Value of a certain Global Attribute, this Attribute
; must be located first and can then be read.
; Retrieve, for example, the product name.
    pt_id=HDF_SD_ATTRFIND(sd_id,"Product")
    HDF_SD_ATTRINFO,sd_id,pt_id,data=plottitle
    print,'    Product: ',plottitle


; The HDF file contains a set of Global Attributes and their Value,
; as well as a number of data sets.
    HDF_SD_FILEINFO,sd_id,datasets,attributes
    nattr=strcompress(attributes)
    nsets=strcompress(datasets)
    print,'    No. of global attributes:',nattr
    print,'    No. of datasets:',nsets


; But doing it this way requires knowing the name of the Attribute.
; The following loop lists the names of the Attributes, with their
; index number (0,1,...,nattr-1). For example the first four.
    print,''
    print,'      id  Global Attribute (first five)'
    print,'     -----------------------------------'
    for id=0,4 do begin
        HDF_SD_ATTRINFO,sd_id,id,NAME=attr_name,DATA=attr_data
        print,id,'  ',attr_name
    endfor


; In order to read the Value of a given Attribute, one needs to know
; the structure of the Value: is it a string, an array of value, etc.
; The routine HDF_SD_ATTRINFO provides the TYPE keyword for this,
; and the COUNT keyword says how many values there are.
; For example the "UV_field_date", being the 3rd Attribute, which
; has index 2 (as the previous list showed).
    print,''
    print,'    Details of the UV_field_date attribute:'
    print,''
    HDF_SD_ATTRINFO,sd_id,2,TYPE=date_type,COUNT=date_count,NAME=date_name
    print,FORMAT="('       ',a,' is of type ',a,'; no. of values: ',i2)", $
          date_name,date_type,date_count


; In other words, it is a bit fiddling around to get the Attribute Values
; read. The date can thus be read like this, assuming the NAME is known,
; but not the index:
    months=strarr(12)
    months=["January","February","March","April","May","June","July", $
            "August","September","October","November","December"]

    pt_id=HDF_SD_ATTRFIND(sd_id,"UV_field_date")
    HDF_SD_ATTRINFO,sd_id,pt_id,data=pd
    year=pd(0)
    nmon=pd(1)
    month=months(nmon-1)
    day=pd(2)
;    hour=pd(3)   ; this is "0" as the UV value is for the whole day
;    minute=pd(4) ; this is "0"
;    second=pd(5) ; this is "0"
    plotdate=string(day,month,year,format="(i2,1x,a,1x,i4)")
    print,FORMAT="('       UV_field_date values stored : ',6i5)",pd
    print,FORMAT="('       UV_field_date as a nice date: ',a)",plotdate
    print,''


; The Global Attributes also give information on the latitude and
; longitude grid of the data sets. For the UV index or dose this is:
;   Number_of_longitudes     720
;   Longitude_range          -179.75  179.75
;   Longitude_step           0.5
;   Number_of_latitudes      360
;   Latitude_range            -89.75   89.75
;   Latitude_step            0.5
;
; As an example, read the number of grid cells.
    lonattr=HDF_SD_ATTRFIND(sd_id,"Number_of_longitudes")
    latattr=HDF_SD_ATTRFIND(sd_id,"Number_of_latitudes")
    HDF_SD_ATTRINFO,sd_id,lonattr,data=nlon
    HDF_SD_ATTRINFO,sd_id,latattr,data=nlat


; Create the coordinates of the grid cells, using that the
; number of grid cells and their size is known.
    dlon=0.5
    dlat=0.5
    nlon=720
    nlat=360
    rlon=-179.75+indgen(nlon)*dlon
    rlat=-89.75+indgen(nlat)*dlat


; Now read the data in the data sets.
; There are two data sets, as was determined above.

; Take the first one, with index 0 (zero). Then find its name
; with HDF_SD_GETINFO, which has several additional keywords,
; e.g. to find the data type, but here we know in advance
; all that there is to know.
    sds_id=HDF_SD_SELECT(sd_id,0)
    HDF_SD_GETINFO,sds_id,NAME=var
    print,'    Reading data set ',var,' ...'

; The first one is the UV index or dose, which is stored as integer after
; multiplication by 100. Store the real UV dose in the 'uvval' array.
    HDF_SD_GETDATA,sds_id,iuvval
    uvval=iuvval/100.

; Do the same for the error in the UV index or dose (this error is the
; error due to the error in the ozone column; errors in the other terms
; of the algorithm computing the UV value are NOT taken into account).
    sds_id=HDF_SD_SELECT(sd_id,1)
    HDF_SD_GETINFO,sds_id,NAME=var
    print,'    Reading data set ',var,' ...'
    HDF_SD_GETDATA,sds_id,iuverr
    uverr=iuverr/100.


; Close the HDF file.
    HDF_SD_END,sd_id

; Ask for
    print,''
    print,'>>> Give a "return" to continue'
    inp=""
    read,inp


; From here on the data can be used for plotting or exporting.
; Show, for example, the values in the range
;   +14.0 < longitude < +16.0
;   +44.0 < latitude  < +44.5
; using a loop over all grid cells (not very smart programming, as the
; indices of the range can be found easily, but it's an example).
    print,''
    print,'    Part of the data as example:'
    print,''
    print,'       ilat  latitude    ilon  longitude   UV_value'
    print,'      ----------------------------------------------'
    for ilat=0,nlat-1 do begin
       lat=rlat(ilat)
       if (lat gt +44.0) and (lat lt +45.5) then begin
          for ilon=0,nlon-1 do begin
             lon=rlon(ilon)
             if (lon gt +14.0) and (lon lt +16.0) then begin
                print,FORMAT="(2x,2(5x,i4,3x,f6.2),6x,f6.3)", $
                      ilat,lat,ilon,lon,uvval(ilon,ilat)
             endif
          endfor
       endif
    endfor
    print,''


; Now extract the value at a user-specified coordinate, or to be more
; precise: the value of the grid cell within which the user-specified
; coordinate lies.

nextpoint:

    print,'>>> Extract the UV value at a certain longitude, latitude coordinate'
    print,'>>> Give the longitude, range: [-180:180] degrees'
    print,'    use -999 to end'
    read,lon
    if lon lt -998 then goto,endprog
    if lon lt -180 or lon gt 180 then begin
       print,'*** Error: value outside range'
       goto,nextpoint
    endif
    print,'>>> Give the latitude, range: [-90:90] degrees'
    read,lat
    if lat lt -90 or lat gt 90 then begin
       print,'*** Error: value outside range'
       goto,nextpoint
    endif

    for ilon=0,nlon-1 do begin
       if abs(lon-rlon(ilon)) le (dlon/2.) then goto,gotilon
    endfor

gotilon:

    for ilat=0,nlat-1 do begin
       if abs(lat-rlat(ilat)) le (dlat/2.) then goto,gotilat
    endfor

gotilat:

    print,FORMAT="('    Given coordinates fall in the grid cell with centre:')"
    print,FORMAT="('      longitude = ',f6.2)",rlon(ilon)
    print,FORMAT="('      latitude  = ',f6.2)",rlat(ilat)
    if uvval(ilon,ilat) gt -0.5 then begin
       print,FORMAT="('    The UV value for this cell is: ',f6.3)",uvval(ilon,ilat)
    endif else begin
       print,FORMAT="('    The UV value for this cell is -1, meaning that')"
       print,FORMAT="('    there is no data available for this grid cell')"
    endelse
    print,''

    goto,nextpoint


; End of the program.
endprog:
    print,''
    end

; ====================================================================