% Example program to read 4 km seawifs GAC data in HDF SDS format
% The chl array is stored as type 'Mapped - Composited', see below
% The 'dimsizes' attribute gives the array dimensions, which are 972 x 833
% 
% P. Brickley peter@grampus.umeoce.maine.edu 8/23/01
% This was run on my dec alpha machine w/ Matlab release 11. 
% Change pathnames where appropriate.

filename='/usr/users/peter/comp-mapCAL2000153.hdf'
SD_id = hdfsd('start',filename,'read')
[ndatasets,nglobal_attr,status] = hdfsd('fileinfo',SD_id)

% for the 4km daily scenes, name1='Mapped - Composited'
% for the 4km 8-day composites, name1='Mapped-Composite'   - an old typo 

name1='Mapped - Composited'			; the chl array is stored by this name
id_u = hdfsd('nametoindex',SD_id,name1);
sds_id1 = hdfsd('select',SD_id, id_u)
[name,rank,dimsizes,data_type,nattrs,status] = hdfsd('getinfo',sds_id1)


% read the chl data array
[chl_0,status] = hdfsd('readdata',sds_id1,[0,0],[1,1],dimsizes);

% read some attributes of the data
att1='units'
attr_id1 = hdfsd('findattr',sds_id1, att1)
[units,status] = hdfsd('readattr',sds_id1,attr_id1)

att2='Slope'
attr_id2 = hdfsd('findattr',sds_id1, att2)
[slope,status] = hdfsd('readattr',sds_id1,attr_id2)

att3='Intercept'
attr_id3 = hdfsd('findattr',sds_id1, att3)
[intercept,status] = hdfsd('readattr',sds_id1,attr_id3)

att4='Limit'
attr_id4 = hdfsd('findattr',sds_id1, att4)
[limit,status] = hdfsd('readattr',sds_id1,attr_id4)

att5='Latitude Center'
attr_id5 = hdfsd('findattr',sds_id1, att5)
[center_lat,status] = hdfsd('readattr',sds_id1,attr_id5)

att6='Longitude Center'
attr_id6 = hdfsd('findattr',sds_id1, att6)
[center_lon,status] = hdfsd('readattr',sds_id1,attr_id6)

% close access to sds array and to HDF file
status = hdfsd('endaccess',sds_id1)
status = hdfsd('end',SD_id)

% Now scale and plot the data

chl_0=double(chl_0);  % convert int16 array into a standard Matlab double format
chl_0=chl_0';		% reverse row column format (hdf convention)

%scale to geophysical chl units

chl=chl_0*0.001 + 32.0;		%scales image to units of mg m^-3 with a range of 0 - 25

imagesc(chl);		% Create a scaled image. You probably want to log10 scale this 
%instead...
