; =============================================================
;
; Example of reading HDF 4.0 data files of the SO2 column densities
; as produced by TEMIS/PROMOTE using IDL.
;
; For information, see the on-line product information page
; http://www.oma.be/BIRA-IASB/Molecules/SO2archive/info/?kind=dataspecidlhdf
;
; (c) TEMIS / PROMOTE / BIRA-IASB
;     last modified: 17 August 2007
;
; =============================================================
;
; Usage within IDL to compile
;
;    IDL> .run so2hdfread
;
; And to run:
;
;    IDL> so2hdfread
;
; after which the routine will ask for
; the name of the HDF file to process.
;
; ====================================================================

pro so2hdfread

; Read the name of the HDF file to be read and remove
; leading and trailing spaces.
; Alternatively, comment out the line "read,filename"
; and give the filename hardcoded.

    print,''
    print,'>>> Give the name of the HDF file to read'

    filename=''
    read,filename

    ;filename='data/so2cd2007031921.hdf'

    filename=strtrim(filename,2)
    print,''


; Check whether the HDF file exists.

    fd_id=FILE_TEST(filename,/READ)
    if fd_id ne 1 then begin
       print,' *** Error: HDF file does not exist'
       goto,endprog
    endif


; 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 with HDF_SD_START'
       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 six.

    print,''
    print,'      id  Global Attribute (first six)'
    print,'     -----------------------------------'

    for id=0,5 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 "SO2_field_date_1" and "SO2_field_date_2", specifying
; the date range for which the plot is valid.
; In principle we know what attribute has which index, but let us
; assume we have to search for the index.

    print,''
    print,'    Details of the SO2_field_date_# attributes:'
    print,''

    fd1_id=HDF_SD_ATTRFIND(sd_id,"SO2_field_date_1")
    HDF_SD_ATTRINFO,sd_id,fd1_id,TYPE=fd1_type,COUNT=fd1_count,NAME=fd1_name
    print,FORMAT="('       ',a,' is of type ',a,'; no. of values: ',i2)", $
          fd1_name,fd1_type,fd1_count

    fd2_id=HDF_SD_ATTRFIND(sd_id,"SO2_field_date_2")
    HDF_SD_ATTRINFO,sd_id,fd2_id,TYPE=fd2_type,COUNT=fd2_count,NAME=fd2_name
    print,FORMAT="('       ',a,' is of type ',a,'; no. of values: ',i2)", $
          fd2_name,fd2_type,fd2_count


; Let us print the contents of the first of these two fields.

    months=strarr(12)
    months=["January","February","March","April","May","June","July", $
            "August","September","October","November","December"]

    HDF_SD_ATTRINFO,sd_id,fd1_id,data=fd1
    year=fd1(0)
    nmon=fd1(1)
    month=months(nmon-1)
    day=fd1(2)
    plotdate1=string(day,month,year,format="(i2,1x,a,1x,i4)")

    print,''
    print,FORMAT="('       ',a,' values stored : ',6i5)",fd1_name,fd1
    print,FORMAT="('       ',a,' as a nice date: ',a)",fd1_name,plotdate1


; The attributes "Data_begin" and "Data_end" have 6 values: the data and
; the time (in UTC). These are read in the same way, but the array read
; has 6 elements.
 

; The Global Attributes also give information on the latitude and
; longitude grid of the data sets.
; As an example, read the number of grid cells.

    lon_id=HDF_SD_ATTRFIND(sd_id,"Number_of_longitudes")
    lat_id=HDF_SD_ATTRFIND(sd_id,"Number_of_latitudes")
    HDF_SD_ATTRINFO,sd_id,lon_id,data=mlon
    HDF_SD_ATTRINFO,sd_id,lat_id,data=mlat

; Unfortunately 'mlon' and 'nlon' can now NOT be used directly to
; create the coordinate arrays of the grid. The reason for this is
; that 'mlon' and 'nlon' appear to be of type LONG, rather than
; being simple integers. But we want simple integers. The following
; does the trick.

    nlon=fix(mlon(0),type=2)
    nlat=fix(mlat(0),type=2)

    print,''
    print,'    Number of longitudes, latitudes:',nlon,nlat


; Setup the latitude and longitude arrays for the centres of
; the grid cells.
; In principle we could read the size of the grid cells and
; the value for the centre of the first grid cell from the
; attributes, but we might as well compute it here.

    dlon=360.0/float(nlon)
    dlat=180.0/float(nlat)
    rlon=(-180.0+dlon/2.0)+indgen(nlon)*dlon
    rlat=( -90.0+dlat/2.0)+indgen(nlat)*dlat

    print,''
    print,'       Grid cell size (lon,lat): ',dlon,dlat
    print,'       First longitude value   :',rlon(0)
    print,'       First latitude  value   :',rlat(0)
    print,''


; Now read the data of the data sets.
; There are 'nsets' data sets, as was determined above.
; While storing the fields, take into account that the values
; in the data sets have been multiplied by 1000.
; As an example, read the slant column and cloud cover data.

    for iset=0,nsets-1 do begin

       sds_id=HDF_SD_SELECT(sd_id,iset)
       HDF_SD_GETINFO,sds_id,name=var

       CASE var OF
       "Iscd_field" : begin
              HDF_SD_GETDATA,sds_id,infield
              print,'    Reading  ',var,' data set ...'
              scdfield=infield/1000.
           end
       "Ivcd_field_2" : begin
              HDF_SD_GETDATA,sds_id,infield
              print,'    Reading  ',var,' data set ...'
              vcdfield=infield/1000.
              break
           end
       "Iccf_field" : begin
              HDF_SD_GETDATA,sds_id,infield
              print,'    Reading  ',var,' data set ...'
              ccffield=infield/1000.
           end
       ELSE : begin
              ; HDF_SD_GETDATA,sds_id,infield
              print,'    Skipping ',var,' data set ...'
              break
           end
       ENDCASE

    end

; The data is successfully read, so close the HDF file.

    HDF_SD_END,sd_id


; Ask for a key stroke to continue.

    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
;  +114.0 < longitude < +115.0
;   +35.0 < latitude  <  +36.0
; using a loop over all grid cells (not very smart programming, as the
; indices of the range can be found easily, but it is just an example).

    print,''
    print,'    Part of the data as example:'
    print,''
    print,'     ilat  latitude  ilon  longitude   SCD [DU]  VCD [DU]  cloud frac.'
    print,'    ------------------------------------------------------------------'

    for ilat=0,nlat-1 do begin

       lat=rlat(ilat)
       if (lat gt 35.0) and (lat lt 36.0) then begin

          for ilon=0,nlon-1 do begin
             lon=rlon(ilon)
             if (lon gt 114.0) and (lon lt 115.0) then begin
                print,FORMAT="(2x,2(3x,i4,3x,f7.3),3(4x,f6.3))", $
                      ilat,lat,ilon,lon,scdfield(ilon,ilat),     $
                      vcdfield(ilon,ilat),ccffield(ilon,ilat)
             endif
          endfor

       endif

    endfor

    print,''


; To disable the next part, enable either
; of the following 'goto' statements.

;   goto,endprog 
;   goto,plotdata
 

; 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 field value at a user-specified location'

    print,'>>> Give the longitude, range: [-180:180] degrees;'
    print,'    use 999 to end this extracting for values.'
    read,lon

    if lon eq 999 then goto,plotdata

    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,''
    print,FORMAT="('    Given coordinates fall in the grid cell with centre:')"
    print,FORMAT="('      longitude = ',f8.3)",rlon(ilon)
    print,FORMAT="('      latitude  = ',f8.3)",rlat(ilat)

    if scdfield(ilon,ilat) gt -90.0 then begin
       print,FORMAT="('    The SCD field value for this cell is: ',f6.3)",scdfield(ilon,ilat)
    endif else begin
       print,FORMAT="('    The SCD field value for this cell is -99, meaning that')"
       print,FORMAT="('    there is no data available for this grid cell')"
    endelse

    print,''

    goto,nextpoint


; Make a simple world plot of the slant column data

plotdata:

    print,''
    print,'>>> Want to see a plot of the SO2 slant column field?'
    print,'    Give "1" (one) for "yes" and something else for "no"'
    read,np
    if np ne 1 then goto,endprog

    print,''
    print,'--- Plotting the SO2 slant column field, which takes a while ...'
    print,'    The colour scale runs from 0.5 DU to 3 DU.'

; Use as colour table "Rainbow+white". 
; The white is used for the background of the plot.
; The first colour of this table is black, which makes it not possible
; to use black as colour for any other line, such as the coast lines.
; Therefore it is better to start higher up in the colour table, e.g.
; at index 32 (some shade of blue) and go up to 254 (red). 

    device, decomposed=0, retain=2
    loadct, 39, /silent

    map_set, /cylindrical, /continents, con_color=0, $
             e_horizon={fill:1, color:255}

    xx=fltarr(5)
    yy=fltarr(5)

    for ilat=0,nlat-1 do begin
    for ilon=0,nlon-1 do begin

       quant = scdfield(ilon,ilat)

       ; Do not plot "no data" grid cells

       if quant gt -90. then begin

          ; Set the corner coordinates of the grid cell

          xx(0)=rlon(ilon)-dlon/2.
          xx(1)=rlon(ilon)+dlon/2.
          xx(2)=rlon(ilon)+dlon/2.
          xx(3)=rlon(ilon)-dlon/2.
          yy(0)=rlat(ilat)-dlat/2.
          yy(1)=rlat(ilat)-dlat/2.
          yy(2)=rlat(ilat)+dlat/2.
          yy(3)=rlat(ilat)+dlat/2.
          xx(4)=xx(0)
          yy(4)=yy(0)

          ; Plot values from 0.5 (blue) to 3 (red) DU

          icol=fix((quant-2.5)*(254-32)/2.5)+32
          if icol lt 32  then icol=32
          if icol gt 255 then icol=254

          polyfill,xx,yy,color=icol

       endif

    endfor
    endfor

    map_continents, /coasts, color = 0, mlinethick=2


; End of the program.

endprog:

    print,''
    end

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