%
% This example code illustrates how to access and visualize LAADS MYD02HKM v6
% HDF-EOS2 Swath file in MATLAB. 
%
%  If you have any questions, suggestions, comments on this example, please 
% use the HDF-EOS Forum (http://hdfeos.org/forums). 
%
%  If you would like to see an  example of any other NASA HDF/HDF-EOS data 
% product that is not listed in the HDF-EOS Comprehensive Examples page 
% (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org or 
% post it at the HDF-EOS Forum (http://hdfeos.org/forums).
%
% Usage: save this script and run (without .m at the end)
% 
% $matlab -nosplash -nodesktop -r MYD02HKM_A2015125_0115_006_2015125153005_hdf
%
% Tested under: MATLAB R2017a
% Last updated: 2018-02-19

import matlab.io.hdfeos.*
import matlab.io.hdf4.*

% Open file.
FILE_NAME='MYD02HKM.A2015125.0115.006.2015125153005.hdf';
file_id = sw.open(FILE_NAME, 'rdonly');

% Open swath.
SWATH_NAME='MODIS_SWATH_Type_L1B';
swath_id = sw.attach(file_id, SWATH_NAME);

% Read the dataset.
DATAFIELD_NAME='EV_500_RefSB';
data_raw = sw.readField(swath_id, DATAFIELD_NAME, [], [], []);

% Detach Swath object.
sw.detach(swath_id);

% Close file.
sw.close(file_id);


% Read attributes from the data field.
SD_id = sd.start(FILE_NAME, 'rdonly');
sds_index = sd.nameToIndex(SD_id, DATAFIELD_NAME);
sds_id = sd.select(SD_id, sds_index);

% Read _FillValue from data field.
fillvalue_index = sd.findAttr(sds_id, '_FillValue');
fillvalue = sd.readAttr(sds_id, fillvalue_index);

% Get the long name from data field.
long_name_index = sd.findAttr(sds_id, 'long_name');
long_name = sd.readAttr(sds_id, long_name_index);

% Read units from the data field.
units_index = sd.findAttr(sds_id, 'reflectance_units');
units = sd.readAttr(sds_id, units_index);

% Read scale_factor from the data field.
scale_index = sd.findAttr(sds_id, 'reflectance_scales');
scale = sd.readAttr(sds_id, scale_index);
scale = double(scale(1));

% Read add_offset from the data field.
offset_index = sd.findAttr(sds_id, 'reflectance_offsets');
offset = sd.readAttr(sds_id, offset_index);
offset = double(offset(1));

% Read valid_range from the data field.
range_index = sd.findAttr(sds_id, 'valid_range');
range = sd.readAttr(sds_id, range_index);

% Terminate access to the corresponding data set.
sd.endAccess(sds_id);

% Close the file.
sd.close(SD_id);

% Subset 3-D data and make 2D.
lev=0;
data=squeeze(data_raw(:,:,lev+1));


% Read lat/lon info from the outputs of  eo2dump file.
lat1D = ...
    load('lat_MYD02HKM.A2015125.0115.006.2015125153005.output');

lon1D = ...
    load('lon_MYD02HKM.A2015125.0115.006.2015125153005.output');
[xdimsize, ydimsize] = size(data);
lat = reshape(lat1D, xdimsize, ydimsize);
lon = reshape(lon1D, xdimsize, ydimsize);

% Convert the data to double type for plot.
data=double(data);
lon=double(lon);
lat=double(lat);

% Replace the filled value with NaN.
data(data==fillvalue) = NaN;
data(data > double(range(2))) = NaN;
data(data < double(range(1))) = NaN;

% Multiply scale and add offset, the equation is scale *(data-offset).
data = scale*(data-offset);

% Set the map parameters.
lon_c = lon(xdimsize/2, ydimsize/2);
lat_c = lat(xdimsize/2, ydimsize/2);
latlim=ceil(max(max(lat))) - floor(min(min(lat)));

% create the graphics figure -- 'visible'->'off' = off-screen rendering
f=figure('Name', FILE_NAME, ...
	 'Renderer', 'zbuffer', ...
	 'Position', [0,0,800,600], ...
	 'visible', 'off');


% FlatLimit will give us a zoom-in effect in Ortho projection.
axesm ('ortho', 'Frame', 'on', 'Grid', 'on', ...
       'FLatLimit', [-Inf, latlim], ...
       'origin', [lat_c, lon_c]);
mlabel('equator');
plabel(0); 
plabel('fontweight','bold');

% Load the coastlines data file.
coast = load('coast.mat');

% Plot coastlines in color black ('k').
plotm(coast.lat,coast.long,'k');

% Plot data. 
lat = lat(:)';
lon = lon(:)';
data = data(:)';

% Use every 5th point to save memory.
% scatterm(lat, lon, 1, data);
scatterm(lat(1:5:end), lon(1:5:end), 1, data(1:5:end));


% Put colormap.
colormap('Jet');
h=colorbar();
set (get(h, 'title'), 'string', units);

% Set the title using long_name.
field_info = hdfinfo(FILE_NAME, 'eos');
dname = field_info.Swath.DataFields(3).Dims(1).Name;
title({FILE_NAME; ...
      ['Reflectance derived from ',long_name ]; ...
      ['at ', strrep(dname,'_','\_'),'=',int2str(lev)]}, ...
      'FontSize',10,'FontWeight','bold');
saveas(f, [FILE_NAME '.m.png']);
exit;