%  This example code illustrates how to access and visualize OMI Swath file
% in MATLAB with Google Earth. 
%
%  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).
%
% Acknowledgement
%
%  The KML conversion part of this example code is provided by
%  Kelly Luekemeyer, Mapping Toolbox Development, MathWorks, Inc.

% Usage: save this script and run (without .m at the end)
%
% $matlab -nosplash -nodesktop -r OMI_L2_OMNO2_CloudFraction_kml
%
% Tested under: MATLAB R2012a
% Last updated: 2013-11-08

clear

% We assume that the KML Toolbox is decompressed under /tmp directory.
addpath('/tmp/kmltoolbox v2.71/');

% Open the HDF5 File.
FILE_NAME = 'OMI-Aura_L2-OMNO2_2008m0720t2016-o21357_v003-2008m0721t101450.he5';
file_id = H5F.open (FILE_NAME, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');

% Open the dataset.
DATAFIELD_NAME = '/HDFEOS/SWATHS/ColumnAmountNO2/Data Fields/CloudFraction';
data_id = H5D.open (file_id, DATAFIELD_NAME);

Lat_NAME='HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Latitude';
lat_id=H5D.open(file_id, Lat_NAME);

Lon_NAME='HDFEOS/SWATHS/ColumnAmountNO2/Geolocation Fields/Longitude';
lon_id=H5D.open(file_id, Lon_NAME);

% Read the dataset.
data=H5D.read (data_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');

lat=H5D.read(lat_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');

lon=H5D.read(lon_id,'H5T_NATIVE_DOUBLE', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT');

% Read the units.
ATTRIBUTE = 'Units';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
units = H5A.read(attr_id, 'H5ML_DEFAULT');

% Read the offset.
ATTRIBUTE = 'Offset';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
offset = H5A.read(attr_id, 'H5ML_DEFAULT');

% Read the scale.
ATTRIBUTE = 'ScaleFactor';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
scale = H5A.read(attr_id, 'H5ML_DEFAULT');

% Read the fill value.
ATTRIBUTE = '_FillValue';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
fillvalue=H5A.read (attr_id, 'H5T_NATIVE_DOUBLE');

% Read the missing value.
ATTRIBUTE = 'MissingValue';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
missingvalue=H5A.read (attr_id, 'H5T_NATIVE_DOUBLE');

% Read title attribute.
ATTRIBUTE = 'Title';
attr_id = H5A.open_name (data_id, ATTRIBUTE);
long_name=H5A.read (attr_id, 'H5ML_DEFAULT');

% Close and release resources.
H5A.close (attr_id)
H5D.close (data_id);
H5F.close (file_id);

% Replace the fill value with NaN.
data(data==fillvalue) = NaN;
% Replace the missing value with NaN.
data(data==missingvalue) = NaN;

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

% Create an overlay image in PNG.
f = figure('Name', FILE_NAME, ...
    'Renderer', 'zbuffer', ...
    'Position', [0,0,1440,720], ...
    'Visible', 'on', ...
    'PaperPositionMode', 'auto', ...
    'Colormap', jet(2048));
ax = axesm('MapProjection', 'eqdcylin', 'Frame', 'on', 'Grid', 'on');
set(ax, ...
    'PlotBoxAspectRatio', [2 1 1], ...
    'Box', 'off', ...
    'Position', [0 0 1 1], ...
    'Units', 'normalized', ...
    'DataAspectRatioMode', 'auto');
tightmap;

surfacem(lat,lon,data);

coast = load('coast.mat');
plotm(coast.lat,coast.long,'k')

saveas(f, [FILE_NAME, '.m.png']);

% Make the white background transparent.
transparent = 255;
A = imread([FILE_NAME, '.m.png']);
alpha = ones(size(A,1), size(A,2));
index = A(:,:,1) == transparent ...
        & A(:,:,2) == transparent ...
        & A(:,:,3) == transparent;
alpha(index) = 0;
imwrite(A, [FILE_NAME, '.m.png'], 'Alpha', alpha);

% Save it as KMZ.
k  = kml([FILE_NAME, '.m.kml']);
k.overlay(180.0, -180.0, -90.0, 90.0, 'file',[FILE_NAME '.m.png']);
k.run;
exit;