Making Enquiries Against a Scientific Data Set

But observe that we've been cheating, just a little, all along. We knew what was on the file and hence we also knew how to read it, e.g., how much space to reserve for the data to be read from the file, and we have hardwired this knowledge into our little program.

The following program still reads the data as before, but it makes some inquiries against both the file and against the data set and prints the result of these inquiries on standard output.

This program contains many new elements, which are all explained below.

Here's the program:

program read_SDS

  use hdf
  use dffunc
  use netcdf

  implicit none

  ! constants

  character(len=7), parameter :: sd_file_name = "SDS.hdf"
  integer, parameter :: number_of_data = 10, sd_set_index = 0
  integer, dimension(1), parameter :: start = (/ 0 /), stride = (/ 1 /), &
       edges = (/ number_of_data /)

  ! variables

  integer :: sd_file_id, sd_set_id, sd_read_status, sd_set_status, &
       sd_file_status, n_datasets, n_file_attributes, file_info_status, &
       rank, data_type, n_set_attributes, i, set_info_status
  integer(kind=4), dimension(number_of_data) :: data = 0
  integer(kind=4), dimension(MAXVDIMS) :: dim_sizes
  character(len=MAXNCNAM) :: sd_set_name

  sd_file_id = sfstart(sd_file_name, DFACC_READ)
     write (*, *) "File ", sd_file_name, " opened"
     write (*, *) "sd_file_id        = ", sd_file_id
  file_info_status = sffinfo(sd_file_id, n_datasets, n_file_attributes)
     write (*, *) "Obtained information about the file"
     write (*, *) "file_info_status  = ", file_info_status
     write (*, *) "n_datasets        = ", n_datasets
     write (*, *) "n_file_attributes = ", n_file_attributes
  sd_set_id = sfselect(sd_file_id, sd_set_index)
     write (*, *) "First scientific data set selected"
     write (*, *) "sd_set_id         = ", sd_set_id
  set_info_status = sfginfo(sd_set_id, sd_set_name, rank, dim_sizes, &
       data_type, n_set_attributes)
     write (*, *) "Obtained information about the data set"
     write (*, *) "sd_set_name       = ", sd_set_name(:10)
     write (*, *) "rank              = ", rank
     write (*, *) "dim_sizes         = ", (/ (dim_sizes(i), i = 1, rank) /)
     write (*, *) "data_type         = ", data_type
     write (*, *) "n_set_attributes  = ", n_set_attributes
     write (*, *) "data (before)     = ", data
  sd_read_status = sfrdata(sd_set_id, start, stride, edges, data)
     write (*, *) "sd_read_status    = ", sd_read_status
     write (*, *) "data (after)      = ", data
  sd_set_status = sfendacc(sd_set_id)
     write (*, *) "sd_set_status     = ", sd_set_status
  sd_file_status = sfend(sd_file_id)
     write (*, *) "sd_file_status    = ", sd_file_status

end program read_SDS

Before we discuss the program have a look at how it compiles, links, and runs:

gustav@blanc:../src 14:44:19 !526 $ make hdf-4
f90 -c -M/afs/ovpit.indiana.edu/@sys/HDF/modules hdf-4.f90
f90 -o hdf-4 hdf-4.o -L/afs/ovpit.indiana.edu/@sys/HDF/lib \
       -lmfhdf -lnsl -ldf -ljpeg -lz -lm
gustav@blanc:../src 14:44:23 !527 $ ./hdf-4
 File SDS.hdf opened
 sd_file_id        = 393216
 Obtained information about the file
 file_info_status  = 0
 n_datasets        = 1
 n_file_attributes = 0
 First scientific data set selected
 sd_set_id         = 262144
 Obtained information about the data set
 sd_set_name       = SDStmpl   
 rank              = 1
 dim_sizes         = 10
 data_type         = 24
 n_set_attributes  = 0
 data (before)     = 10*0
 sd_read_status    = 0
 data (after)      = 1,  2,  3,  4,  5,  6,  7,  8,  9,  10
 sd_set_status     = 0
 sd_file_status    = 0
gustav@blanc:../src 14:44:26 !528 $

Observe the output shorthand:

data (before)     = 10*0

This does not mean $10 \times 0$. It means 10 zeros in a row.

The first thing you should notice is that this time we have also used module netcdf. This module contains definitions of some constants, which are used in the program: most notably MAXVDIMS  and MAXNCNAM .

The next novelty I would like to draw your attention to is the initialisation of array data. This time we don't use an implicit DO loop. Instead we simply say:

integer(kind=4), dimension(number_of_data) :: data = 0

Whenever an ANSI Fortran compiler sees a constant assigned to an array, it results in assigning this constant to every element of the array.

The program begins as before with opening the file "SDS.hdf". The next step, however, is new: 

  file_info_status = sffinfo(sd_file_id, n_datasets, n_file_attributes)
     write (*, *) "Obtained information about the file"
     write (*, *) "file_info_status  = ", file_info_status
     write (*, *) "n_datasets        = ", n_datasets
     write (*, *) "n_file_attributes = ", n_file_attributes

Here we call function sffinfo, whose interface is defined in module dffunc, The function takes as it argument the file ID, and returns the total number of scientific data sets saved on the file (in our case this is just one) and the total number of the so called file attributes. These are simply annotations and we shall learn how to add them to the file soon enough. In HDF the file can have general annotations, and every scientific data sets can have its own annotations (or attributes) too. For the time being the number of file attributes returned by sffinfo is zero. The function returns exit status, which we save on file_info_status. When the exit status is zero, it means everything is fine. When it is -1 it means that an error has occurred. Module hdf defines the following: 

hdf.f90:      parameter(SUCCEED         = 0, FAIL       = -1)

After having made a call to sffinfo we again proceed as before, i.e., we select the first scientific data set. But then there is a new call to function sfginfo: 

  set_info_status = sfginfo(sd_set_id, sd_set_name, rank, dim_sizes, &
       data_type, n_set_attributes)
     write (*, *) "Obtained information about the data set"
     write (*, *) "sd_set_name       = ", sd_set_name(:10)
     write (*, *) "rank              = ", rank
     write (*, *) "dim_sizes         = ", (/ (dim_sizes(i), i = 1, rank) /)
     write (*, *) "data_type         = ", data_type
     write (*, *) "n_set_attributes  = ", n_set_attributes
     write (*, *) "data (before)     = ", data

This function takes as its argument the scientific set ID and returns

sd_set_name

name of the scientific data set, in this case it is "SDStmpl". In general a scientific data set name can be up to MAXNCNAM characters long, which is why we have defined

character(len=MAXNCNAM) :: sd_set_name

Observe that when we write that name on standard output, we select only its first 10 characters: sd_set_name(:10). In Fortran strings are not terminated with a null character as they are in C and C++. Instead they are padded with spaces. So if we were to write simply sd_set_name, we'd end up writing 7 characters followed by MAXNCNAM - 7 spaces, and module netcdf defines

parameter(MAXNCNAM = 128)

rank

is a rank of the data set, i.e., 2 for a 2-D matrix, 3 for a 3-D cube, and so on. In our case rank is 1, because we have written a 1-dimensional array on the data set.

dim_sizes

is an array of rank 1, whose length must be at least equal to the rank of the scientific data set. The maximum length of that array, MAXVDIMS, is defined on module netcdf:

parameter(MAXNCDIM = 32)
...
parameter(MAXVDIMS = MAXNCDIM)

That is, the scientific data sets can have up to 32 dimensions. On exit function sfginfo fills every slot of that array with a number that corresponds to the length of the scientific data set in an appropriate direction. For example, if the data set has dimensions $32\times 48$ then dim_sizes(1) = 32 and dim_sizes(2) = 48. In principle this array is going to be rank entries long. So when we write its contents on standard output we use the implicit DO loop:

write (*, *) "dim_sizes         = ", (/ (dim_sizes(i), i = 1, rank) /)

In our case this loop returns only one entry: 10.

data_type

this is the numerical value of a data type that fills the scientific data set. Our program returns number 24, which corresponds to

parameter(DFNT_INT32      = 24)

as defined on module hdf.

n_set_attributes

this is a number of the scientific data set attributes, which is like a number of the file attributes, i.e., a number of the data set annotations. For the time being this number remains zero.

The remainder of the program works as before.