c readPFV5_clim.f - Read Pathfinder V5 SST climatologies and display c ---------------------------------------------------------------------------- c c Example fortran program for using the fortran HDF API to read in c AVHRR Pathfinder Version 5 SST Climatologies created by NODC. c Each climatology contains the following layers (with some small differences c depending on the averaging period chosen): c c Clim_SST c contains the climatological SST without any extra gap filling attempted, c convert from pixel value to degrees C using scale=0.075 and offset of -3.0 c (SST in deg C = scale*pixel_Value + offset) c c Clim_StandardDeviation: c the standard deviation of those averages, convert to degrees C c using scale=0.075 and offset=0 c c Clim_Counts c the number of observations going into those averages, no conversion needed c c Longitude and Latitude c these are vectors containing the longtiude and latiude values, where c negative values are S and W, positive values N and E. c c Clim_SST_Filled c scale=0.075, offset=-3,0, a gap filled version of Clim_SST c derived using the methods of Casey and Cornillon, 1999 (Journal of c Climate, vol 12, 1848-1863). This layer is the one intended for most c users of the data. The Annual and seasonal climatologies do not c contain this layer c c For access to the climatologies: c c ftp://data.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.0_Climatologies c c For more information on the climatologies see the README file at: c c ftp://data.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.0_Climatologies/README.txt c c compile with: c f77 -o readPFV5_clim readPFV5_clim.f -L/home/kcasey/HDF4/lib -lmfhdf -ldf -ljpeg -lz c c where /home/kcasey/HDF4/lib is replaced with the location of your HDF c libraries. c c Then, run simply with: c c read_PFV5_clim c c There is a lot of output to the screen, so you may wish to run with: c c read_PFV5_clim > out.txt c c In either case, three ASCII text files will be created showing a c subset of the data contained in the HDF file. These three files are: c example_output_climsst.txt - contains every 512th pixel value from c the dataset "Clim_SST", which is the c mean SST climatology before gap-filling c example_output_climstd.txt - Same as above, but contains the c climatology standard deviation c example_output_climfill.txt - Same as first, but contains the c gap-filled version c c The assistance of Takaya Namba is gratefully acknowledged, c along with the the example HDF code found at: c http://hdf.ncsa.uiuc.edu/training/UG_Examples/SD c c By: c Kenneth S. Casey, NOAA National Oceanographic Data Center. c Kenneth.Casey@noaa.gov c August 2005 c c c ---------------------------------------------------------------------------- program readPFV5_clim c Variable declarations: integer sd_id,sds_id,sds_index,file_id integer sfselect, sfendacc, sffinfo integer sfginfo, sfstart, sfend, sfrdata, sfrcdata integer sfgainfo, sfgcal, sfrcatt, sfrnatt integer start(2),edges(2),stride(2) integer status integer n_file_attributes integer n_datasets integer dim_sizes(2) integer rank integer data_type, n_values integer n_set_attributes, attr_index integer n_dim_attributes integer dataset_number character*56 sd_set_name character*80 attr_name character*900 attr_val_char real*8 attr_val_num(20) character*20 dim_name integer*2 climsst(8192,4096) integer*2 climstd(8192,4096) integer*2 climfill(8192,4096) real*8 latitude(4096) real*8 longitude(8192) integer*2 value real*8 scale_factor, scale_factor_err real*8 add_offset, add_offset_err c Parameter declarations: character*30 FILE_NAME character*30 OUTPUT_FILE_CLIMSST character*30 OUTPUT_FILE_CLIMSTD character*30 OUTPUT_FILE_CLIMFILL parameter(FILE_NAME ='month01_combined.hdf') parameter(OUTPUT_FILE_CLIMSST='example_output_climsst.txt') parameter(OUTPUT_FILE_CLIMSTD='example_output_climstd.txt') parameter(OUTPUT_FILE_CLIMFILL='example_output_climfill.txt') parameter(DFACC_RDONLY = 1) parameter(DFNT_CHAR = 4) parameter(DFNT_NUM = 6) c Open access to the HDF-SDS file. After starting acces, print c the sd_id value to the screen. c If it opened properly, the value will be a positve integer. c If it failed to open properly, the value will be -1. write(*,*) " " write(*,*) "====================================================" sd_id=sfstart(FILE_NAME,DFACC_RDONLY) if (sd_id .eq. -1) then write(*,*) "Failed to open access to ", FILE_NAME write(*,*) " File may not exist. Exiting..." stop else write(*,*) "Successfully opened ", FILE_NAME endif c Now get some information from the HDF file: n_datasets is the c number of datasets in this HDF file, and n_file_attributes c indicates how many global metadata attributes are in the file. After c getting this information, print the status of the get information c request. A status value of 0 indicates success, while -1 is failure. status=sffinfo(sd_id,n_datasets,n_file_attributes) write(*,*) "Number of datasets in this file = ", n_datasets write(*,*) "Number of global attributes = ", n_file_attributes c Now print out global attributes of the HDF-SDS file: write(*,*) " " write(*,*) "-------------------------------------------------" write(*,*) " Attribute Name Length Value" write(*,*) "-------------------------------------------------" do 10 attr_index = 0, n_file_attributes-1 status = sfgainfo(sd_id, attr_index, attr_name, data_type, + n_values) if (data_type .eq. DFNT_CHAR) then attr_val_char="" status = sfrcatt(sd_id, attr_index, attr_val_char) if (n_values .LE. 30) then write(*,*) " ", attr_name, n_values, + " ", attr_val_char(1:30) else write(*,*) " ", attr_name, n_values, + " ", attr_val_char(1:30), "..." endif else c status = sfrnatt(sd_id, attr_index, attr_val_num) c I know all these have integer values, even though some are c stored as real, so it is ok to convert to int for display: c write(*,111) attr_name, n_values, c + (int(attr_val_num(i)),i=1,n_values) 111 format (2X,A,1X,I2,4X,20I4) endif 10 continue write(*,*) "-------------------------------------------------" write(*,*) " " c Now select the dataset we want to work with. HDF numbering begins c with 0, so the you can select a value from 0 to n_datasets. A positive value c of sds_id indicates the selection was successful, while -1 indicates failure. c Note that at this stage we are only pointing to the dataset... we c have not actually read it in yet. In the Pathfinder Version 5 c climatology files: c c Clim_SST = 0 c Clim_StandardDeviation = 1 c Clim_Counts = 2 c Longitude = 3 c Latitude = 4 c Clim_SST_Filled = 5 c c First, we use 0 to select the Clim_SST array, which contains the c climatological SST before any gap filling has been performed: dataset_number=0 sds_id=sfselect(sd_id,dataset_number) c Now we get information on the selected dataset. c sd_set_name = the name of the dataset c rank = number of dimensions in the dataset c dim_sizes = the number of elements in each dimension c data_type = the type of data stored in the dataset (integer, etc.) c n_set_attributes = the number of metadata attributes associated c with this dataset (this is different than the c number of global attributes identified above) c After getting this information, print some of it to the screen. status=sfginfo(sds_id,sd_set_name,rank,dim_sizes, + data_type,n_set_attributes) write(*,*) "-----------------------------------" write(*,*) " Dataset name: ", sd_set_name write(*,*) " Dataset number = ", dataset_number write(*,*) " Number of dimensions = ", rank write(*,*) " dim_sizes(1) = ", dim_sizes(1) write(*,*) " dim_sizes(2) = ", dim_sizes(2) write(*,*) " Number of attributes = ", n_set_attributes write(*,*) " " write(*,*) "-------------------------------" write(*,*) " Attribute Name Length" write(*,*) "-------------------------------" do 20 attr_index = 0, n_set_attributes-1 status = sfgainfo(sds_id, attr_index, attr_name, data_type, + n_values) write(*,*) " ",attr_name, n_values 20 continue write(*,*) "-----------------------------------" write(*,*) " " c Although the scale and offset needed to convert the pixel c values into degrees are contained as attributes for this c dataset, they are also included in a special HDF-defined c attribute which has its own access function. Let's use c this to get the scale and offset for this dataset: status = sfgcal(sds_id, scale_factor, scale_factor_err, + add_offset, add_offset_err, data_type) write(*,*) " scale_factor = ", scale_factor write(*,"(A,F4.1)") " add_offset = ", add_offset write(*,*) " " write(*,*) " Remember to multiply by scale_factor and " write(*,*) " add add_offset to convert from pixel to deg C!" write(*,*) " " c Now prepare to do the reading of the data. c xl=number of elements in x [dim_sizes(1)] c yl=number of elements in y [dim_sizes(2)] c start(1)=where in x to start reading the array (0 based) c start(2)=where in y to start reading the array (0 based) c stride(1)=whether to read all values in x (stride(1)=1) or to skip c and read only a subsampling of all the values (for example, a c stride(1)=2 will skip and read only every other value) c stride(2)=same as above, but for y c edges(1)=how far to read in x c edges(2)=how far to read in y c These elements are all much better explained in the HDF 4 Users Guide. xl = dim_sizes(1) yl = dim_sizes(2) start(1)=0 start(2)=0 stride(1)=1 stride(2)=1 edges(1)=xl edges(2)=yl c Here we do the actual reading of the data, using the above parameters. c As in previous steps, status is 0 for success -1 for failure. Note c that here we use sfrcdata since climsst is integer*2 type data. Below, c when reading real*8 data, we'll use sfrdata. (in C, the same function c is used for both) write(*,*) " Reading data... " status=sfrcdata(sds_id,start,stride,edges,climsst) c Now, print out the values to an ascii file so you can be confident in c the values you have read in. In this example, we are writing out only c every 512th value, even though all values were read in. Write c out the converted pixel values: write(*,*) " Example output (every 512th clim. SST) has " write(*,*) " been written to ", OUTPUT_FILE_CLIMSST write(*,*) " " open(99,file=OUTPUT_FILE_CLIMSST) write(99,"(16(f6.3,1x))")((climsst(i,j)*scale_factor+add_offset, + i=1,dim_sizes(1),512),j=1,dim_sizes(2),512) c Now, end access to the data set we have been working with (sds_id). Since c this HDF file has more data sets in it, we can end access to this c data set and begin working with the others. So, we will now end access c to the current data set. status=sfendacc(sds_id) c Ok, if all of the above looks rather confusing, this time we will c simply select the latitude and longitude datasets, read them in, c print a few values to the screen, then end access to them without c any further commentary. Note to read these data sets, we will use c sfrdata, since they are real*8. Note that having Latitude and Longitude c included as distinct data sets in an HDF file is not really c necessary, since you can give names and values to the dimensions c of datasets HDF. c However, some people find it easier to have them as separate c datasets, so I have included them as well. So, without further c commentary, select, read, display, and close: do 30 dataset_number = 3,4 sds_id=sfselect(sd_id,dataset_number) status=sfginfo(sds_id,sd_set_name,rank,dim_sizes, + data_type,n_set_attributes) write(*,*) "-----------------------------------" write(*,*) " " write(*,*) " Dataset name = ", sd_set_name write(*,*) " dataset number = ", dataset_number write(*,*) " dim_sizes(1) = ", dim_sizes(1) write(*,*) " Example output (every 512th value):" start(1)=0 stride(1)=1 edges(1)=dim_sizes(1) write(*,*) " Reading data... " if (dataset_number .eq. 1) then status=sfrdata(sds_id,start,stride,edges,latitude) write(*," (8(F8.3,x))")(latitude(i),i=1,dim_sizes(1),512) else status=sfrdata(sds_id,start,stride,edges,longitude) write(*," (16(F8.3,x))")(longitude(i),i=1,dim_sizes(1),512) endif write(*,*) " " status=sfendacc(sds_id) 30 continue write(*,*) "-----------------------------------" c Now also read in the STDV, and send output to one text file: dataset_number=1 sds_id=sfselect(sd_id,dataset_number) status=sfginfo(sds_id,sd_set_name,rank,dim_sizes, + data_type,n_set_attributes) status = sfgcal(sds_id, scale_factor, scale_factor_err, + add_offset, add_offset_err, data_type) write(*,*) " " write(*,*) " Dataset name = ", sd_set_name write(*,*) " dataset number = ", dataset_number write(*,*) " dim_sizes(1) = ", dim_sizes(1) write(*,*) " dim_sizes(2) = ", dim_sizes(2) write(*,*) " scale_factor = ", scale_factor write(*,"(A,F4.1)") " add_offset = ", add_offset start(1)=0 start(2)=0 stride(1)=1 stride(2)=1 edges(1)=dim_sizes(1) edges(2)=dim_sizes(2) write(*,*) " Reading data... " status=sfrcdata(sds_id,start,stride,edges,climstd) write(*,*) " Example output (every 512th clim. stdv.) has " write(*,*) " been written to ", OUTPUT_FILE_CLIMSTD write(*,*) " " open(97,file=OUTPUT_FILE_CLIMSTD) write(97,"(16(f6.3,1x))")((climstd(i,j)*scale_factor+add_offset, + i=1,dim_sizes(1),512),j=1,dim_sizes(2),512) status=sfendacc(sds_id) c Now also read in the filled Clim. SST, and send output to one text file: dataset_number=5 sds_id=sfselect(sd_id,dataset_number) status=sfginfo(sds_id,sd_set_name,rank,dim_sizes, + data_type,n_set_attributes) status = sfgcal(sds_id, scale_factor, scale_factor_err, + add_offset, add_offset_err, data_type) write(*,*) " " write(*,*) " Dataset name = ", sd_set_name write(*,*) " dataset number = ", dataset_number write(*,*) " dim_sizes(1) = ", dim_sizes(1) write(*,*) " dim_sizes(2) = ", dim_sizes(2) write(*,*) " scale_factor = ", scale_factor write(*,"(A,F4.1)") " add_offset = ", add_offset start(1)=0 start(2)=0 stride(1)=1 stride(2)=1 edges(1)=dim_sizes(1) edges(2)=dim_sizes(2) write(*,*) " Reading data... " status=sfrcdata(sds_id,start,stride,edges,climfill) write(*,*) " Example output (every 512th clim. stdv.) has " write(*,*) " been written to ", OUTPUT_FILE_CLIMFILL write(*,*) " " open(90,file=OUTPUT_FILE_CLIMFILL) write(90,"(16(f6.3,1x))")((climfill(i,j)*scale_factor+add_offset, + i=1,dim_sizes(1),512),j=1,dim_sizes(2),512) status=sfendacc(sds_id) c Now terminate access to the entire HDF file, and print the status: status=sfend(sd_id) write(*,*) "-----------------------------------" write(*,*) " Done. Goodbye. " write(*,*) " " write(*,*) "====================================================" write(*,*) " " stop end