function [imagedata] = modis_level1b_read(filename, band, keywords, area)

%FUNCTION [IMAGEDATA] = MODIS_LEVEL1B_READ(FILENAME, BAND, KEYWORDS, AREA)
%
% NAME: 
%    MODIS_LEVEL1B_READ
%
% PURPOSE:
%    Read a single band from a MODIS Level 1B HDF product file at
%    1000, 500, or 250 m resolution.
%
%    The output image is available in the following units 
%    (Radiance is default):
%    RAW DATA:    Raw data values as they appear in the HDF file
%    CORRECTED:   Corrected counts
%    RADIANCE:    Radiance (Watts per square meter per steradian per micron)
%    REFLECTANCE: Reflectance (Bands 1-19 and 26; 
%                 *without* solar zenith correction)
%    TEMPERATURE: Brightness Temperature (Bands 20-25 and 27-36, Kelvin)
%
%    This procedure uses only HDF calls (it does *not* use HDF-EOS),
%    and only reads from SDS and Vdata arrays 
%    (it does *not* read ECS metadata).
%
% CATEGORY:
%    MODIS utilities.
%
% CALLING SEQUENCE:
%    MODIS_LEVEL1B_READ(FILENAME, BAND)
%
% INPUTS:
%    FILENAME     (string)  Name of MODIS Level 1B HDF File
%    BAND         (integer) Band number to be read
%                 (1-36 for 1000 m, 1-7 for 500 m, 1-2 for 250m)
%                 (Use 13, 13.5, 14, 14.5 for 13lo, 13hi, 14lo,
%                 14hi)
% OPTIONAL INPUTS:
%    KEYWORDS =   cellstring of keywords such as:
%                 {'reflectance', 'temperature', 'raw', 'corrected',
%                  'latitude', 'longitude'}.
%          KEYWORD PARAMETERS:
%          RAW         If set, image data are returned as raw HDF values
%                      (default is to return image data as radiance).
%          CORRECTED   If set, image data are returned as corrected counts
%                      (default is to return image data as radiance).
%          REFLECTANCE If set, image data for bands 1-19 and 26 only are
%                      returned as reflectance *without* solar
%                      zenith angle correction
%                      (default is to return image data as radiance).
%          TEMPERATURE If set, image data for bands 20-25 and
%                      27-36 only are returned as brightness temperature
%                      (default is to return image data as
%                      radiance).
%          LATITUDE    On return, an array containing the reduced 
%                      resolution latitude data for the entire
%                      granule (degrees, every 5th line and pixel).
%          LONGITUDE   On return, an array containing the reduced 
%                      resolution longitude data for the entire
%                      granule (degrees, every 5th line and pixel)
%    AREA        (optional)  a four element vector specifying the area
%                 to be read [XO, YO, NX, NY].
%
% OUTPUTS:
%    IMAGEDATA    A struct with the format: 
%             .MIMAGE     a two dimensional array of image data in
%                         the requested units
%             .UNITS      a string describing the image units
%             .PARAMETER  a string describing the image (e.g. 'Radiance')
%	      .BAND       band number integer 
%             .FILETYPE   a string describing the file type (e.g. 
%                         'MOD021KM')
%             .DATESTRING a string describing the image date
%                         (e.g. '2000092')
%             .TIMESTRING a string describing the image time
%                         (e.g. '1710' for 17:10 UTC) 
%             .FILENAME   the HDF filename (e.g. 
%                         'MOD021KM.A2000092.1710.002.2000094041335.hdf'
%
% OPTIONAL OUTPUTS:
%    IMAGEDATA.LATITUDE   (see above)
%    IMAGEDATA.LONGITUDE  (see above)
%
% RESTRICTIONS:
%    Only tested on Matlab version 5.3.0 (R11) - performance under
%    other versions of Matlab is unknown
%
%    Requires the following functions by Shaima.Nasiri@ssec.wisc.edu
%    MODIS_BRIGHT
%    BRITE_M  (within modis_bright.m file)
% 
% EXAMPLES:
%
% Read band 1 in radiance units from a 1000 m resolution file
% filename = 'MOD021KM.A2000062.1020.002.2000066023928.hdf'
% imagedata = modis_level1b_read(filename, 01)
%
% Read band 31 in temperature units from a 1000 m resolution file
% filename = 'MOD021KM.A2000062.1020.002.2000066023928.hdf'
% imagedata = modis_level1b_read(filename, 31, {'temperature'})
%
% read a band 1 subset in reflectance units from a 1000 m resolution
% filename = 'MOD021KM.A2000062.1020.002.2000066023928.hdf'
% imagedata = modis_level1b_read(filename, 01, {'reflectance'},...
%                                [20 20 400 400])
%
% read a band 2 subset in reflectance units from a 250 m resolution
% filename = 'MOD02QKM.A2000161.1420.002.2000163022951.hdf'
% imagedata = modis_level1b_read(filename, 01, {'reflectance'},...
%                                [200 200 600 600])
%  MODIFIED from Liam Gumley's modis_level1b_read.pro IDL routine
%
%  


% Check Inputs
 nargs = nargin;
 if (nargs < 2) 
    error(['Usage:  MODIS_LEVEL1B_READ(FILENAME,BAND)']);
 elseif (nargs == 2)
    keywords = {' '};
 end
 if ( (nargs <= 4) & (~iscellstr(keywords)) )
       error(['Third argument must be a cellstring of keywords']);
 end
 if (nargs == 4)
    if (~isnumeric(area) | (length(area) ~= 4))
       error(['AREA argument  must be 4 element vector of the form' ...
	      ' [XO YO NX NY]']);
    end
 end
 keywords = lower(keywords);
 
    
%- Check for valid MODIS L1B HDF file, and get file type (1km, 500 m, or 250 m)
 if (~exist(filename,'file'))
    error(['Filename : ' filename ' was not found']);
 end

 len_filename = length(filename);
 if (~strcmp( filename(len_filename-3:len_filename), '.hdf'))
    error(['Filename: ' filename ' is not an HDF file']);
 end

% open hdf file for reading
 sd_id = hdfsd( 'start', filename, 'read' );
 if (sd_id < 0)
    error(['MODIS L1B HDF file is NOT found: ' filename]);
 end
 
 
% find out what sort of datafile we are looking at
 index_1km = hdfsd('nametoindex',sd_id,'EV_1KM_Emissive');
 index_500 = hdfsd('nametoindex',sd_id,'EV_500_RefSB');
 index_250 = hdfsd('nametoindex',sd_id,'EV_250_RefSB');
 
 
 status = hdfsd('end',sd_id);
 if (status < 0)
    warning(['HDF file' filename ' was not closed.']);
 end

 
 if (index_1km > 0)
    filetype = 'MOD021KM';
 elseif (index_500 > 0)
    filetype = 'MOD02HKM';
 elseif (index_250 > 0)
    filetype = 'MOD02QKM';
 else
    error(['Filename ' filename ' is not a Level 1B HDF file']);
 end
% disp(['     Filetype is ' filetype]);

%- check band number and keywords which depend on band number
% check band number
 switch filetype
  case 'MOD021KM'  
   if ( (band < 1) | (band > 36) )
     error([ 'BAND range is 1-36 for this MODIS type => ' filetype ]);
   end
  case 'MOD02HKM'  
   if ( (band < 1) | (band > 7) )
     error([ 'BAND range is 1-7 for this MODIS type => ' filetype ]);
   end
  case 'MOD02QKM'  
   if ( (band < 1) | (band > 2) )
     error([ 'BAND range is 1-2 for this MODIS type => ' filetype ]);
   end
  otherwise
   error(['invalid filetype => ' filetype])
 end

% check for invalid request for reflectance units
 if ( ((band >= 20) & (band ~= 26)) & ismember('reflectance',keywords) )
    error('Reflectance units valid for bands 1-19, 26 only');
 end
% check for invalid request for temperature units
 if ( ((band <= 19) | (band == 26)) & ismember('temperature',keywords) )
    error('Temperature units valid for bands 20-25, 27-36 only');
 end
 

% Set variable name for image data 
 switch filetype
  case 'MOD021KM'
     if ( (band >= 1) & (band <= 2) )
        sds_name = 'EV_250_Aggr1km_RefSB';
     elseif ( (band >= 3) & (band <= 7) )
        sds_name = 'EV_500_Aggr1km_RefSB';
     elseif ( (band >= 8) & (band <= 19) )
        sds_name = 'EV_1KM_RefSB';
     elseif ( (band == 26) )
        sds_name = 'EV_Band26';
     elseif ( (band >= 20) & (band <= 36) )
        sds_name = 'EV_1KM_Emissive';
     end
  case 'MOD02HKM'      
     if ( (band >= 1) & (band <= 2) )
        sds_name = 'EV_250_Aggr500_RefSB';
     elseif ( (band >= 3) & (band <= 7) )
        sds_name = 'EV_500_RefSB';
     end
  case 'MOD02QKM'      
     sds_name = 'EV_250_RefSB';
 end %switch

 % Set attribute names for images data
 fillvalue_name = '_FillValue';
 if ismember('reflectance', keywords)   
    scale_name  = 'reflectance_scales';
    offset_name = 'reflectance_offsets';
    units_name  = 'reflectance_units';
    parameter   = 'Reflectance';
 elseif ismember('corrected', keywords)  
    scale_name  = 'corrected_counts_scales';
    offset_name = 'corrected_counts_offsets';
    units_name  = 'corrected_counts_units';
    parameter   = 'Corrected Counts';
 else
    scale_name  = 'radiance_scales';
    offset_name = 'radiance_offsets';
    units_name  = 'radiance_units';
    parameter   = 'Radiance Counts';
 end
 
% Open the file again and open the sds dataset
 sd_id = hdfsd( 'start', filename, 'read' );
 if (sd_id < 0)
    error(['HDF file' filename ' was not opened.']);
 end
 sds_index = hdfsd('nametoindex',sd_id, sds_name);
 sds_id = hdfsd('select',sd_id, sds_index) ;
 if (sds_id < 0)
    error(['HDF file' filename ' dataset ' sds_name ' was not opened.'])
 end
  
  

% Check band order and get band index
 if (band ~= 26) 
 % Get the actual band order
   band_name = 'band_names' ;
   attr_index = hdfsd('findattr',sds_id, band_name);
   [band_data,status] = hdfsd('readattr',sds_id,attr_index) ;
   if (status < 0)
      status = hdfsd('endaccess',sds_id);
      status = hdfsd('end',sd_id);
      error(['Attribute ' band_data ' not found.']);
   end
 % Set expected band order
   if ( (band >= 1 ) &  (band <= 2) )
      band_order = '1,2';
   elseif ( (band >= 3) & (band <= 7) )
      band_order = '3,4,5,6,7';
   elseif ( (band >= 8) & (band <= 19) )
      band_order = '8,9,10,11,12,13lo,13hi,14lo,14hi,15,16,17,18,19,26';
   else  %( (band >= 8) & (band <= 19) )
      band_order = '20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36';
   end		    
 
 % check actual band order against expected band order
   if ~(strcmp(band_data,band_order))
      status = hdfsd('endaccess',sds_id);
      status = hdfsd('end',sd_id);
      error(['Unexpected band order in image array']);
   end
   
 % get band index
   if ( (band >= 1 ) &  (band <= 2) )
      band_index = band - 1;
   elseif ( (band >= 3 ) &  (band <= 7) )
      band_index = band - 3;
   elseif ( (band >= 8 ) &  (band <= 12) )
      band_index = band - 8;
   elseif ( (band >= 13 ) &  (band <= 14) )
      band_index = 2*band - 21;
   elseif ( (band >= 15 ) &  (band <= 19) )
      band_index = band - 6;
   elseif ( (band >= 20 ) &  (band <= 25) )
      band_index = band - 20;
   elseif ( (band >= 27 ) &  (band <= 36) )
      band_index = band - 21;
   end %if
 
 else  %(band = 26)
    band_index = 0;
 end

% Get Info about the image array
 [varinfo.name, varinfo.rank, varinfo.dims, varinfo.data_type, ...
  varinfo.nattrs, varinfo.status] = hdfsd('getinfo',sds_id);
 if isempty(varinfo.name)
    status = hdfsd('endaccess',sds_id);
    status = hdfsd('end',sd_id);
    error(['Image array was not found: ' sds_name ]);
 end
 if (band ~= 26)
    nbands = varinfo.dims(1);
    npixels_across = varinfo.dims(2);
    npixels_along = varinfo.dims(3);
 else %(band == 26)
    npixels_across = varinfo.dims(1);
    npixels_along = varinfo.dims(2);
 end
 
 
% set start and count values
 if (band ~= 26)
    start = [band_index ; 0 ; 0 ];
    count = [1 ; 1 ; 1];
    edge = [1 ; npixels_across ;  npixels_along];
 else %(band == 26)
    start = [0 ; 0];
    count = [1 ; 1];
    edge = [npixels_across ; npixels_along];
 end

% if Area argument is supplied
 if ( (nargs == 4) & (band == 26) )
    start(1) = min( [max([area(1) 0]) (npixels_across - 1)] );
    start(2) = min( [max([area(2) 0]) (npixels_along - 1)] );
    edge(1) = min( [max([area(3) 0]) (npixels_across - start(1))] );
    edge(2) = min( [max([area(4) 0]) (npixels_along - start(2))] );
 elseif ( (nargs == 4) & (band ~= 26) )
    start(2) = min( [max([area(1) 0]) (npixels_across - 1)] );
    start(3) = min( [max([area(2) 0]) (npixels_along - 1)] );
    edge(2) = min( [max([area(3) 0]) (npixels_across - start(2))] );
    edge(3) = min( [max([area(4) 0]) (npixels_along - start(3))] );
 end

 
 
% read the image array
  [data,status] = hdfsd('readdata',sds_id,start,count,edge);
  %disp(['status:  ' num2str(status)]);
  image_temp = data;
 
 
% read scale attribute    
 attr_index = hdfsd('findattr',sds_id, scale_name);
 [scale, status] = hdfsd('readattr',sds_id,attr_index) ;
 if (status < 0)
    warning(['Attribute ' scale_data ' not found.']);
 end

% read offset attribute    
 attr_index = hdfsd('findattr',sds_id, offset_name);
 [offset, status] = hdfsd('readattr',sds_id,attr_index) ;
 if (status < 0)
    warning(['Attribute ' offset_data ' not found.']);
 end

% read units attribute    
 attr_index = hdfsd('findattr',sds_id, units_name);
 [units, status] = hdfsd('readattr',sds_id,attr_index) ;
 if (status < 0)
    warning(['Attribute ' units_data ' not found.']);
 end

% read units attribute    
 attr_index = hdfsd('findattr',sds_id, units_name);
 [units, status] = hdfsd('readattr',sds_id,attr_index) ;
 if (status < 0)
    warning(['Attribute ' units_data ' not found.']);
 end
  
% read fill value
attr_index = hdfsd('findattr',sds_id, fillvalue_name);
[fillvalue, status] = hdfsd('readattr',sds_id,attr_index);
if (status < 0)
  warning(['Attribute ' fillValue_data ' not found.']);
end
%disp( ['fillvalue = ' num2str(double(fillvalue))] );

% Close the sds dataset 
 status = hdfsd('endaccess',sds_id) ;
 if (status < 0)
    warning(['HDF file' filename ' dataset ' sds_name ' was not closed.']);
 end
 
% Latitude and longitude stuff
 if (ismember('latitude',keywords))
   index_lat = hdfsd('nametoindex',sd_id,'Latitude'); 
   sds_id = hdfsd('select',sd_id, index_lat);
   [name,rank,dimsizes,data_type,nattrs,status] = ...
       hdfsd('getinfo', sds_id);
   [data,status] = hdfsd('readdata',sds_id,[0 ;0],[1;1],dimsizes) ;
   status = hdfsd('endaccess',sds_id);
   if (status < 0)
      warning(['HDF file' filename ' dataset Latitude was not closed.']);
   end

   if index_1km > 0
     resDiff = 1;
   elseif index_500 > 0
     resDiff = 2;
   else
     resDiff = 1;
   end

   % find fill value index
   bad_lat_ind = find(data == -999);
   data(bad_lat_ind) = NaN;

   % interpolate if needed
   tmpData = zeros(size(data));
   if resDiff > 1
     dims = resDiff * size(data);
     tmpData = interpArray(data, dims);
     clear dims;
   else
     tmpData = data;
   end

   % get data in the selected sub-section
   if (nargs == 4)
     latitude = tmpData(area(2):area(2)+area(4)-1, area(1):area(1)+area(3)-1);
   else
     latitude = tmpData;
   end

   bad_lat_ind = find(isnan(latitude) == 1);
   latitude(bad_lat_ind) = fillvalue;

   clear tmpData; clear data; clear bad_lat_ind;
 end

 if (ismember('longitude',keywords))
   index_lat = hdfsd('nametoindex',sd_id,'Longitude'); 
   sds_id = hdfsd('select',sd_id, index_lat);
   [name,rank,dimsizes,data_type,nattrs,status] = ...
       hdfsd('getinfo', sds_id);
   [data,status] = hdfsd('readdata',sds_id,[0 ;0],[1;1],dimsizes) ;
   status = hdfsd('endaccess',sds_id);
   if (status < 0)
      warning(['HDF file' filename ' dataset Longitude was not closed.']);
   end

   if index_1km > 0
     resDiff = 5;
   elseif index_500 > 0
     resDiff = 2;
   else
     resDiff = 1;
   end

   % find fill value index
   bad_lon_ind = find(data == -999);
   data(bad_lon_ind) = NaN;

   % interpolate if needed
   tmpData = zeros(size(data));
   if resDiff > 1
     dims = resDiff * size(data);
     tmpData = interpArray(data, dims);
     clear dims;
   else
     tmpData = data;
   end

   % get data in the selected sub-section
   if (nargs == 4)
     longitude = tmpData(area(2):area(2)+area(4)-1, area(1):area(1)+area(3)-1);
   else
     longitude = tmpData;
   end

   bad_lon_ind = find(isnan(longitude) == 1);
   latitude(bad_lon_ind) = fillvalue;

   clear tmpData; clear data; clear bad_lon_ind;
 end

% Close the file 
 status = hdfsd('end',sd_id);
 if (status < 0)
    warning(['HDF file' filename ' was not closed.']);
 else
 %   disp(['HDF file ' filename ' successfully closed']);
 end

% adjust the image data
% Convert to requested units
 ind_fillvalue = find(image_temp < 0 | image_temp > 32767);
 if ismember('raw',keywords)
 %leave as HDF values
    units = 'Unsigned 16-bit integers';
    parameter = 'Raw HDF Values';
    mimage = image_temp;
 else
 % convert image to unscaled valued   
    mimage = double(scale(band_index+1)) * ...
	    ( double(image_temp) - double(offset(band_index+1)) );

    % Convert Radiance to Brightness Temperature for IR bands
    if ismember('temperature',keywords)
       mimage = modis_bright(mimage, band, 1);
       units = 'Kelvin';
       parameter = 'Brightness Temperature';
       ind = find(mimage < 0);
       mimage(ind) = fillvalue;
    end
 end
 mimage(ind_fillvalue) = fillvalue;
 
% Get data range (min/max of valid image values)
% Not yet coded

% get information from filename
% assuming filename of format:
% 'MOD021KM.A2000092.1710.002.2000094041335.hdf'
 dot_ind = findstr('.', filename);
 platform= filename(dot_ind(1)+1);


 datestring = filename(dot_ind(1)+2 : dot_ind(2)-1);
 timestring = filename(dot_ind(2)+1 : dot_ind(3)-1);

 imagedata = struct('mimage',[], 'units',[], 'parameter',[], ...
		    'band',[], 'filetype',[], ...
		    'datestring',[], 'timestring',[], 'filename',[], ...
                    'fillvalue', []);
% have to flip and rotate image -- NOT
%imagedata.mimage = flipud(rot90(mimage));
 imagedata.mimage = mimage;
 imagedata.units = units;
 imagedata.parameter = parameter;
 imagedata.band = band;
 imagedata.filetype = filetype;
% imagedata.platform = platform;  % not including at this time
 imagedata.datestring = datestring;
 imagedata.timestring = timestring;
 imagedata.filename = filename;
 imagedata.fillvalue = fillvalue;
 if ismember('longitude',keywords)
%   imagedata.longitude = flipud(rot90(longitude));
    imagedata.longitude = longitude;
 end
 if ismember('latitude',keywords)
%   imagedata.latitude = flipud(rot90(latitude));
    imagedata.latitude = latitude;
 end
