In an earlier post, I talked about how to download MODIS data in NetCDF format. Now I will tell you how to process the downloaded data in MATLAB.

Note that there can be be data outages on certain days so you might see less data than expected. For example, MODIS Aqua data for the year 2008 will have only 364 days of data although 366 (2008 is a leap year) is expected. Check the day stamp in the file name to know exactly which day is missing. I will be processing the 2008 data below. Further, I have noticed that on certain days (2-3 days in a year), the filenames are extra long with some added information during data processing in Giovanni. If you see such unusual names, you have to manually rename the files to make the names consistent with other files.

Read data

  1. Copy all the files in your MATLAB workspace.
  2. Since there are so many daily files, we need to read all the files using a loop. So follow the following steps: 
file_names = dir ('*.nc'); % extract filenames and associated information
num_files = length(file_names); %num_files in this case is 364

deep_aod = zeros (360, 180, num_files); % 180 and 360 are latitude and longitude grids in each MODIS data
for m = 1:num_files;
deep_aod (:, :, m) = ncread(file_names(m).name, 'Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean');
end;

The above code will create a 3-d workspace variable named deep_aod (360*180*364) in which 360, 180, and 364 represents longitude, latitude, and days of year 2008, respectively.

Calculate average for plotting

% change the dimensions of latitude and longitude to make it more intuitive (we want to see latitude in the rows and longitudes in the columns)

deep_aod_arng = permute(deep_aod, [2 1 3]); 

%bring the time dimension to the front for easier processing
deep_aod_arng_time = permute(deep_aod_arng, [3 1 2]); 

% convert the 3-d data to 2-d data so that we can apply vector operations in MATLAB, the resulting rows represent each day and columns represent each grid points in the entire globe.
deep_aod_arng_points = reshape(deep_aod_arng_time, [364 64800]); 

Upon examining the data files we see that the days 107 and 309 days are missing. So lets insert two arrays in the respective place so that the date represents complete year and it is easier for further processing. 

%create a dummy row of NaN to be inserted

dummy = NaN (1, 64800); 

%insert the dummy row at the respective rows

deep_aod_arng_points_complete = [deep_aod_arng_points (1:106, :);dummy;deep_aod_arng_points (107:308, :);dummy;deep_aod_arng_points (309:end, :)];

Now create a timetable using this data; it will make it very easy for further calculations, for example, for calculating monthly mean.

%first create a datetime array for 2008

t1 = datetime (2008, 01, 01);
t2 = datetime (2008, 12, 31);
t = (t1:days:t2)';

% now create the timetable

deep_aod_2008_timetable = timetable(t, deep_aod_arng_points_complete);

%now lets calculate the monthly mean AOD values using this timetable

aod_monmean_2008 = retime (deep_aod_2008_timetable, 'monthly', 'mean'); %

Plot the data

% To plot the data, we need to read the latitude and longitude data from any one MODIS file:

MYD_lat = ncread('MYD08_D3.A2008001.061.2018031095350.hdf.nc', 'YDim');
MYD_lon = ncread('MYD08_D3.A2008001.061.2018031095350.hdf.nc', 'XDim');

% create a 2-d grid of latitude and longitude to be used for plotting

lat_mat = repmat(MYD_lat, 1, 360);
lon_mat = repmat(MYD_lon', 180, 1);

Now plot data with the below script using geoshow.

h1 = figure('position', [50 50 1200 600]); 
axesm mercator
p = worldmap([-90 90],[-180 180]);
%borders('countries', 'Color', 'black', 'linewidth', 0.75) 
gridm('on');
setm (p, 'mlabellocation', 40); % for reference_longitude
setm (p, 'plabellocation', 20); % for reference_latitude
geoshow('landareas.shp', 'linewidth', 1.2, 'facecolor', 'none')
h2 = geoshow (lat_mat, lon_mat, reshape(aod_monmean_2008.deep_aod_arng_points_complete(1, :), [180 360]), 'displaytype', 'surface');
caxis([0, 1]);        
title('MODIS AOD Jan 2008');
colormap jet; brighten(0.5);
h3 = colorbar;
set(h3,'fontsize', 12);

You should get a plot like below:

Note: always verify the accuracy of the plotted data from your intuition or secondary sources. For example, in this case, I know that bodele in central Africa is one of the most active dust source region so it should have high AOD which I can see in this plot.

Thank you!