% This example Matlab m-file demonstrates how to read in Pathfinder
% Version 5 data from NODC.  http://pathfinder.nodc.noaa.gov
%
% Created by: Kenneth S. Casey, NOAA National Oceanographic Data Center
%             Kenneth.Casey@noaa.gov
%
% Date: August 17, 2005
%
% Note:  this example code is rather simple in that it does not use
% the full HDF API found in Matlab.  To take advantage of its capabilities,
% including doing things like reading in subsets, striding the data,
% extracting more metadata from the file, etc., look at the help pages 
% on "hdfsd".

% Define the overall quality flag file name:
fname1 = '200112.m04m3pfv50-qual.hdf';

% Now read the quality data using hdfread.  You could also use full HDF API
% with matlab's "hdfsd" functions, but this is simpler. 
qual=hdfread(fname1, 'qual', 'index', {[1.0 1.0],[1.0 1.0], [4096 8192]});

% Now get the information about the hdf file and its contents. Again,
% you could use the "hdfsd" function, but this is simpler:
fileinfo = hdfinfo(fname1);  
lat=fileinfo.SDS.Dims(1).Scale; % Get the latitudes
lon=fileinfo.SDS.Dims(2).Scale; % Get the longitudes

% Now display it:
figure(1)
imagesc(lon,lat,qual); % Display the data
set(gca,'ydir','normal'); % Simply puts the display with positive lats up
title('Overall Quality Levels (0 worst to 7 best)')
caxis([0 7])
colorbar

% NOTE: The overall quality flags do not need any conversion from pixel
% values.  They are stored simply as values from 0 to 7 (in other words,
% the scale is 1 and the offset is 0).  However, if you were reading in one
% of the Pathfinder SST files, you would need to apply a scale and offset
% after reading in the data. For example:
fname2 = '200112.s04m3pfv50-sst-16b.hdf';  

% Read the data in as 16-bit unsigned integers:
sst=hdfread(fname2, 'sst', 'index', {[1.0 1.0],[1.0 1.0], [4096 8192]});

% Get the scale and offset:
fileinfo2 = hdfinfo(fname2);  
scale_factor = double(fileinfo2.SDS.Attributes(11).Value);  % = 0.075;
add_offset= double(fileinfo2.SDS.Attributes(12).Value);     % = -3.0;

% Now apply them to turn SST pixel values into deg C:
sst = double(sst)*scale_factor + add_offset;

% Now, let's subset these data to a particular region,
% apply the quality flag to the data and display it:

lonmin = -100;  % Set lon and lat limts.  Note that if you
                % want to subset across the dateline
                % to create a map of the Pacific, for example,
                % you have to be a little more clever and subset the map
                % into two pieces, one from, for example, 100E to 180E,
                % then another from -180W to 100W.  Then piece
                % the two together.
lonmax = 0;
latmin = 0;
latmax = 45;
goodlons = lon>=lonmin & lon<=lonmax;  % Find matching longitudes
goodlats = lat>=latmin & lat<=latmax; % Find matching latitudes
subset_qual = qual(goodlats,goodlons);  % Subset the quality data
subset_sst = sst(goodlats,goodlons);  % Subset the SST data
cloudypixels = subset_qual<4; % Finds all pixels with quality level < 4
masked_sst = subset_sst;
masked_sst(cloudypixels) = -3;  % Sets the cloud pixels to -3 deg C 
                                % (colder than physically possible)
sublons = lon(goodlons);  % Subset the longitudes
sublats = lat(goodlats); % Subset the latitudes
figure(2)
imagesc(sublons,sublats,masked_sst);  % Display it
set(gca,'ydir','normal'); % Simply puts the display with positive lats up
title('SST in DegC, Quality Levels 4-7 Only')
colorbar; % Add a color scale

% If you like to use the Mapping Toolbox,  do something like 
% this.  Here, we will plot the global SST field with cloudy pixels
% unmasked. Note that this step can take a while depending on the
% memory, CPU, and graphics capability of your computer:
figure(3)
ml = [8192/360 90 -180];  % Sets up the map legend
axesm('robinson');  % or whatever projection you want
sst = flipud(sst); % Puts the array into a more workable orientation
lat2 = fliplr(lat);   % Makes sure you flip the lats around too, in case
                      % you wish to use these later along with "sst"
m=meshm(sst,ml); % Display it
ll = coast;  % Load coarse resolution global coastlines
p=plotm(ll(:,1),ll(:,2),'w');  % Display coastlines on map
caxis([-3 35])
colorbar('horiz')
title('SST in Deg C, All Quality Levels')