How to read a MODIS HDF4 file using fortran ?

Following previous article on How to read a MODIS HDF file using python let's see how to do the same thing in fortran:
Note: HDF libraries have to be installed in your system first (go to the link: obtaining the hdf software )

For example, let's consider the following MODIS HDF file : "MYD06_L2.A2008001.2040.006.2013341205855.hdf".

How to print all SDS names ?

Example of fortran source code (called here how_to_read_modis_hdf4.f90) to get all Scientific Datasets (SDS) names stored in the HDF file (see HDF user guide for more details):

program how_to_read_modis_hdf4

!----------------------------------------------------------------------------------------!

implicit none

integer i

character*300 file_name

integer DFACC_READ, DFNT_INT32

parameter ( DFACC_READ = 1, DFNT_INT32 = 24 )

integer sfstart, sfselect, sfrdata, sfendacc, sfend

integer sffinfo , sfginfo

integer sd_id, sds_id, sds_index

integer status, n_attrs

integer n_datasets, n_file_attrs , index

integer rank, data_type

integer MAX_NC_NAME , MAX_VAR_DIMS

parameter ( MAX_NC_NAME = 256 , MAX_VAR_DIMS = 32 )

character*( MAX_NC_NAME ) name

integer dim_sizes( MAX_VAR_DIMS )

!----------------------------------------------------------------------------------------!

call getarg(1,file_name)

write(6,*) 'file_name : ' , file_name

!----------------------------------------------------------------------------------------!
! Print all SDS names

sd_id = sfstart(file_name, DFACC_READ)

status = sffinfo( sd_id , n_datasets , n_file_attrs )

write(6,*) '! Number of data sets in the file and Number of file attributes :'
write(6,*) 'sffinfo :', status , n_datasets , n_file_attrs
write(6,*)

do index = 0, n_datasets - 1

    sds_id = sfselect(sd_id, index)
    status = sfginfo(sds_id, name, rank, dim_sizes, data_type,n_attrs)

    write(6,*) "name : ", name(1:100), "rank : ", rank, &
             & "index : ", index, "dimension sizes are : ", (dim_sizes(i), i=1, rank), & 
             &  "data type is : ", data_type, "number of attributes is : ", n_attrs

end do

!----------------------------------------------------------------------------------------!
! Close hdf file

status = sfendacc(sds_id)
status = sfend(sd_id)

!----------------------------------------------------------------------------------------!

end program how_to_read_modis_hdf4

!----------------------------------------------------------------------------------------!

Explanations: sfstart function open the HDF file; sffinfo function gives the number of SDS (n_datasets) and attributes (n_file_attrs) associated with the HDF file; sfselect function select a SDS using its index; sfginfo function returns information about the selected SDS and finally sfend close the file. To compile the code:

gfortran how_to_read_modis_hdf4.f90 -I/path_to_hdf/include/ -L/path_to_hdf/lib/ -lmfhdf -ldf -ljpeg -lz -o how_to_read_modis_hdf4

(replace "path_to_hdf" by the full path where HDF libraries are installed on your system) and to execute the code

how_to_read_modis_hdf4 MYD06_L2.A2008001.2040.006.2013341205855.hdf

The above code returns:

 ! Number of data sets in the file and Number of file attributes :
 sffinfo :           0         127          14

 name : Latitude                                                                                            rank :            2 index :            0 dimension sizes are :          270         408 data type is :            5 number of attributes is :           10
 -------------------------
 name : Longitude                                                                                           rank :            2 index :            1 dimension sizes are :          270         408 data type is :            5 number of attributes is :           10
 -------------------------
 name : Scan_Start_Time                                                                                     rank :            2 index :            2 dimension sizes are :          270         408 data type is :            6 number of attributes is :           10
 -------------------------
 name : Solar_Zenith                                                                                        rank :            2 index :            3 dimension sizes are :          270         408 data type is :           22 number of attributes is :           10
 -------------------------
 name : Solar_Zenith_Day                                                                                    rank :            2 index :            4 dimension sizes are :          270         408 data type is :           22 number of attributes is :           10
 -------------------------
 name : Solar_Zenith_Night   
.
.
.

sffinfo function returns 127 (meaning that MODIS C6 cloud optical products HDF files MYD/MOD06 contain 127 Scientific Datasets (SDS)

How to get the data ?

Suppose we know the name of the SDS we are interested in (for exemple "cloud_top_temperature_1km"). To find the index of the SDS "cloud_top_temperature_1km":

index = sfn2index(sd_id, 'cloud_top_temperature_1km')

To select the SDS

sds_id = sfselect(sd_id, index)

To get info about the selected SDS:

status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )

write(6,*) "index : ", index, "dimension sizes are : ", (dim_sizes(i), i=1, rank), & 
    &   "data type is : ", data_type, "number of attributes is : ", n_attrs

we got:

index :           60 dimension sizes are :         1354        2040 data type is :           22 number of attributes is :           10

meaning than cloud_top_temperature_1km SDS is a (1354,2040) array, data type is 22 corresponding to a 16-bit integer type in fortran (to see how to translate data type HDF Data Types and Fortran 90 kind parameter ) and has 10 attributes (which contain additional information about the SDS such as the valid range, the unit, etc).

Now to read the SDS, using sfrdata function:

integer*2, dimension(:,:), allocatable :: cloud_top_temperature_1km

allocate(cloud_top_temperature_1km(dim_sizes(1),dim_sizes(2)),STAT=AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
status = sfrdata( sds_id, start, stride, edges, cloud_top_temperature_1km )

Print cloud_top_temperature_1km values:

do i = 1, 1354
do j = 1, 2040
write(6,*) cloud_top_temperature_1km(i,j)
end do
end do

gives

Note: Don't forget to deallocate the array if you don't want to go out of memory pretty quickly:

deallocate(cloud_top_temperature_1km, STAT = DeAllocateStatus)
if(DeAllocateStatus/=0) STOP

How to get the SDS attributes ?

cloud_top_temperature_1km SDS has 10 attributes, to get information about each attribute, there is the sfgainfo function:

do attr_index = 0, n_attrs - 1
status = sfgainfo(sds_id, attr_index, attr_name, data_type, n_values)
write(6,*) attr_name, data_type, n_values
end do

results

 valid_range                   22           2
 _FillValue                    22           1
 long_name                      4         123
 units                          4           1
 scale_factor                   6           1
 add_offset                     6           1
 Parameter_Type                 4           6
 Cell_Along_Swath_Sam          24           3
 Cell_Across_Swath_Sa          24           3
 Geolocation_Pointer            4          34

Let's say we want the 'scale_factor' and 'add_offset', the data type is 6 for both corresponding to a 64-bit floating-point type in fortran (see hdf data types and Fortran 90 kind parameter ):

real(KIND=8) scale_factor
real(KIND=8) add_offset

to select the attribute use the sffattr function and sfrattr to get selected attribute value. Example for 'scale_factor' attribute:

attr_index = sffattr(sds_id, 'scale_factor')
write(6,*) attr_index
status = sfrattr(sds_id, attr_index, scale_factor)
write(6,*) scale_factor

same for 'add_offset':

attr_index = sffattr(sds_id, 'add_offset')
write(6,*) attr_index
status = sfrattr(sds_id, attr_index, add_offset)
write(6,*) add_offset

Print cloud top temperature at 1km:

do i = 1, 1354
do j = 1, 2040
write(6,*) ( cloud_top_temperature_1km(i,j) - add_offset) * scale_factor 
end do
end do

results

242.45999458059669     
242.45999458059669     
242.45999458059669     
244.43999453634024     
244.43999453634024     
244.43999453634024     
244.43999453634024     
244.43999453634024

How to extract info from a byte ?

To extract info from a byte there is the fortran function ibits. Example: how to read Cloud_Mask_1km SDS stored in MOD/MYD06 cloud optical products (see MODIS C6 cloud optical products user guide page 115). The Cloud Mask 1km SDS has two bytes (1byte = 8 bits, so 16 bits for Cloud Mask 1km). Lets consider we want to get the "Cloud Mask Status Flag" and "Cloud Mask Cloudiness Flag". Both are stored in the first byte (see page 115 of MODIS C6 cloud optical products user guide ):

Cloud_Mask_1km(1,i,j)

where i,j are the pixel indexes of along track and cross track respectively. To extract "Cloud Mask Status Flag" and "Cloud Mask Cloudiness Flag" using ibits function:

cloud_mask_status_flag =  ibits(Cloud_Mask_1km(1,i,j),0,1)
cloud_mask_cloudiness_flag = ibits(Cloud_Mask_1km(1,i,j),1,2)

ibits function needs the position of the first bit and the number of bits used to code the information (0-based indexing):

Byte 1 ( Cloud_Mask_1km(1,i,j) )

Flag Name	First bit index	Number of Bits
Cloud Mask Status Flag	0	1
Cloud Mask Cloudiness Flag	1	2
Day/Night Flag	3	1
Sunglint Flag	4	1
Snow/Ice Flag	5	1
Surface Type Flag	6	2

Byte 2 ( Cloud_Mask_1km(2,i,j) )

Flag Name	First bit index	Number of Bits
Heavy Aerosol Flag	0	1
Thin Cirrus Flag	1	1
Shadow Flag	2	1

In summary: Cloud Mask Cloudiness Flag is stored in the first byte of Cloud_Mask_1km SDS (Cloud_Mask_1km(1,i,j)), first bit index = 1 and number of bits = 2: ibits(Cloud_Mask_1km(1,i,j),1,2).

See full code below.

Example of source code :

program how_to_read_modis_hdf4

!----------------------------------------------------------------------------------------!

implicit none

integer i

character*300 file_name

integer DFACC_READ, DFNT_INT32

parameter ( DFACC_READ = 1, DFNT_INT32 = 24 )

integer sfstart, sfselect, sfrdata, sfendacc, sfend

integer sffinfo , sfginfo

integer sd_id, sds_id, sds_index

integer status, n_attrs

integer n_datasets, n_file_attrs , index

integer rank, data_type

integer MAX_NC_NAME , MAX_VAR_DIMS

parameter ( MAX_NC_NAME = 256 , MAX_VAR_DIMS = 32 )

character*( MAX_NC_NAME ) name

integer dim_sizes( MAX_VAR_DIMS )

character * 100 :: sds_name

integer start(3), edges(3), stride(3)

integer sfn2index

integer :: AllocateStatus, DeAllocateStatus

! data type 22 => integer 16 bits

integer*2, dimension(:,:), allocatable :: cloud_top_temperature_1km

integer sffattr, sfrattr, sfgainfo

integer      attr_index, n_values

character*20 attr_name

real(KIND=8) scale_factor

real(KIND=8) add_offset

! data type 20 => integer 8 bits

integer*1, dimension(:,:,:), allocatable :: Cloud_Mask_1km

integer cloud_mask_status_flag, cloud_mask

!----------------------------------------------------------------------------------------!

call getarg(1,file_name)

write(6,*) 'file_name : ' , file_name

!----------------------------------------------------------------------------------------!
! Get all SDS names

sd_id = sfstart(file_name, DFACC_READ)

status = sffinfo( sd_id , n_datasets , n_file_attrs )

write(6,*) '! Number of data sets in the file and Number of file attributes :'
write(6,*) 'sffinfo :', status , n_datasets , n_file_attrs
write(6,*)

do index = 0, n_datasets - 1

    sds_id = sfselect(sd_id, index)
    status = sfginfo(sds_id, name, rank, dim_sizes, data_type,n_attrs)

    write(6,*) "name : ", name(1:100), "rank : ", rank, &
             & "index : ", index, "dimension sizes are : ", (dim_sizes(i), i=1, rank), & 
             &  "data type is : ", data_type, "number of attributes is : ", n_attrs

end do

!----------------------------------------------------------------------------------------!
! Example of how to extract the data

index = sfn2index(sd_id, 'cloud_top_temperature_1km')

write(6,*) 'index', index

sds_id = sfselect(sd_id, index)
status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )

sds_name = name(1:len(name))

write(6,*) "index : ", index, "dimension sizes are : ", (dim_sizes(i), i=1, rank), & 
        &   "data type is : ", data_type, "number of attributes is : ", n_attrs

do i = 1 , rank
start(i) = 0
edges(i) = dim_sizes(i)
stride(i) = 1
end do

allocate(cloud_top_temperature_1km(dim_sizes(1),dim_sizes(2)),STAT=AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
status = sfrdata( sds_id, start, stride, edges, cloud_top_temperature_1km )

do i = 1, 100
    write(6,*) cloud_top_temperature_1km(i,100)
end do

!----------------------------------------------------------------------------------------!
! Example of how to get SDS attributes

do attr_index = 0, n_attrs - 1
    status = sfgainfo(sds_id, attr_index, attr_name, data_type, n_values)
    write(6,*) attr_name, data_type, n_values
end do

attr_index = sffattr(sds_id, 'scale_factor')
write(6,*) attr_index
status = sfrattr(sds_id, attr_index, scale_factor)
write(6,*) scale_factor

attr_index = sffattr(sds_id, 'add_offset')
write(6,*) attr_index
status = sfrattr(sds_id, attr_index, add_offset)
write(6,*) add_offset

do i = 1, 10
    write(6,*) ( cloud_top_temperature_1km(i,100) - add_offset) * scale_factor
end do

!----------------------------------------------------------------------------------------!
! Example of how to extract information from a Byte

index = sfn2index(sd_id, 'Cloud_Mask_1km')

write(6,*) 'Cloud_Mask_1km index', index

sds_id = sfselect(sd_id, index)
status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )

sds_name = name(1:len(name))

do i = 1 , rank
start(i) = 0
edges(i) = dim_sizes(i)
stride(i) = 1
end do

allocate(Cloud_Mask_1km(dim_sizes(1),dim_sizes(2),dim_sizes(3)),STAT=AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory Cloud_Mask_1km***"
status = sfrdata( sds_id, start, stride, edges, Cloud_Mask_1km )

write(6,*) 'READ MODIS Cloud Mask Example : '

do i = 1, 100

    cloud_mask_status_flag = ibits(Cloud_Mask_1km(1,i,1),0,1)
    cloud_mask = ibits(Cloud_Mask_1km(1,i,1),1,2)

    write(6,*) cloud_mask_status_flag, cloud_mask

end do

!----------------------------------------------------------------------------------------!
! deallocate

deallocate(cloud_top_temperature_1km, STAT = DeAllocateStatus)
if(DeAllocateStatus/=0) STOP

deallocate(Cloud_Mask_1km, STAT = DeAllocateStatus)
if(DeAllocateStatus/=0) STOP

!----------------------------------------------------------------------------------------!
! Close hdf file

status = sfendacc(sds_id)
status = sfend(sd_id)

!----------------------------------------------------------------------------------------!

end program how_to_read_modis_hdf4

!----------------------------------------------------------------------------------------!

References

Link	WebSite
MODIS C6 cloud optical products user guide	NASA
HDF user guide	HDF group
HDF Data Types	HDF group
Fortran 90 kind parameter	stackoverflow
Example code	HDF Group
IBITS — Bit extraction	gcc
Language Reference	tamu.edu

Developer in my spare time, I’m interested in Earth observation and satellite-based fire detection. Passionate about Python, machine learning, and open science.