Load one HDF file of ocean color data (e.g., chlorophyll, PAR, SST)

load_color_1file.m

function [this] = load_color_1file(filename,lonlim,latlim);
%LOAD_OCEANCOLOR loads hdf-format ocean color files (including PAR,
% chlorophyll and Aqua SST
%
% [THIS] = LOAD_OCEANCOLOR(FILENAME,LONLIM,LATLIM) loads the SeaWiFS level 3
% style HDF-file given by filename and optional LONLIM and LATLIM, which
% are 2 element vectors ranging from -180:180 for longitude and -90:90
% for latitutde
%
% Removes values out of range, e.g.
% For chlorophyll, removes values which are greater than 1.81001 IF file 
% doesn't have values greater than 1.82 - this is for chlorophyll values 
% which have no data at the highest value
%
% Outputs a structure with fields data, lat and lon.

% Default region is global
if nargin<3
    latlim=[-90 90];
    if nargin<2
        lonlim = [-180 180];
    end
end

%% Load the file and trim to lonlim latlim
PI=hdfinfo(filename);

% And write it into a structure
pin=[];
for k=1:58,
    nm=PI.Attributes(k).Name;nm(nm==' ')='_';
    if isstr(PI.Attributes(k).Value),
        pin=setfield(pin,nm,PI.Attributes(k).Value);
    else
        pin=setfield(pin,nm,double(PI.Attributes(k).Value));
    end
end;

% lon/lat of grid corners
lon = [pin.Westernmost_Longitude:pin.Longitude_Step:pin.Easternmost_Longitude];
lat = [pin.Northernmost_Latitude:-pin.Latitude_Step:pin.Southernmost_Latitude];

% Get the indices needed for the area of interest
[mn,ilt]=min(abs(lat-max(latlim)));
[mn,ilg]=min(abs(lon-min(lonlim)));
ltlm=fix(diff(latlim)/pin.Latitude_Step);
lglm=fix(diff(lonlim)/pin.Longitude_Step);

% load the subset of data needed for the map limits given
P=hdfread(filename,'l3m_data','Index',{[ilt ilg],[],[ltlm lglm]});
P=double(P);
P(P==255)=NaN;
P=(pin.Slope*P+pin.Intercept);   % log_10 of chla

%% Remove values out-of-range
% Convert data into log(Chla) using the equations given. Blank no-data.
if max(P(:))<1.82
    if max(P(:))>1.8
        P(find(P>=1.81001))=NaN;
    end
end
% for Aqua-modis 11mu sst daytime, remove the nans
if max(P(:))<45.1
    if max(P(:))>45
        P(find(P>45))=NaN;
    end
end
% Eliminate PAR values greater than 76.2
if max(P(:))>76
    if max(P(:))<76.3
        P(find(P>76.2))=NaN;
    end
end
% Eliminate k490 values greater than 0.794
if max(P(:))>0.79
    if max(P(:))<0.80
        P(find(P>0.794))=NaN;
    end
end

%% Reformat the data to be LONGxLATxTIME
LT = lat(ilt+[0:ltlm-1]);
LG=lon(ilg+[0:lglm-1]);

%% Output structure
this.data = P;
this.lon = LG;
this.lat = LT;