% Read height fields and coordinate arrays
hgt = double(ncread('hgt.2010.nc','hgt'));
lat = double(ncread('hgt.2010.nc','lat'));
lon = double(ncread('hgt.2010.nc','lon'));
time = double(ncread('hgt.2010.nc','time'));
level = double(ncread('hgt.2010.nc','level'));

% Set up map
figure; worldmap([30, 90], [-90, 45]);

% Get 600 mb height surface on April 16
hgt600 = double(hgt(:,:,5,106));

% Read ash plume image for April 16 2010 and map pixels onto 
% latitudes and longitudes
ash = imread('16042010_00z.png');
[rows, columns, channels] = size(ash);
R = georefcells([-90,90],[-180,180], [rows, columns], ...
    'ColumnsStartFrom','north');

% Plot ash plume
geoshow(ash,R);

% Plot coastlines
load coastlines;
plotm(coastlat, coastlon, 'k');

% Plot height surface
[C, h] = contourm(lat,lon,hgt600',3000:50:5000,'k');
clabelm(C, h, 'manual');

% Set title and save
title('16 April 2010');
exportgraphics(gca, '1604.png');

% Repeat for April 17
hgt600 = double(hgt(:,:,5,107));
ash = imread('17042010_00z.png');
[rows, columns, channels] = size(ash);
R = georefcells([-90,90],[-180,180], [rows, columns], ...
    'ColumnsStartFrom','north');
geoshow(ash,R);
plotm(coastlat, coastlon, 'k');
[C, h] = contourm(lat,lon,hgt600',3000:50:5000,'k');
clabelm(C, h, 'manual');
title('17 April 2010');
exportgraphics(gca, '1704.png');

% Repeat for April 18
hgt600 = double(hgt(:,:,5,108));
ash = imread('18042010_00z.png');
[rows, columns, channels] = size(ash);
R = georefcells([-90,90],[-180,180], [rows, columns], ...
    'ColumnsStartFrom','north');
geoshow(ash,R);
plotm(coastlat, coastlon, 'k');
[C, h] = contourm(lat,lon,hgt600',3000:50:5000,'k');
clabelm(C, h, 'manual');
title('18 April 2010');
exportgraphics(gca, '1804.png');

% Repeat for April 19
hgt600 = double(hgt(:,:,5,109));
ash = imread('19042010_00z.png');
[rows, columns, channels] = size(ash);
R = georefcells([-90,90],[-180,180], [rows, columns], ...
    'ColumnsStartFrom','north');
geoshow(ash,R);
plotm(coastlat, coastlon, 'k');
[C, h] = contourm(lat,lon,hgt600',3000:50:5000,'k');
clabelm(C, h, 'manual');
title('19 April 2010');
exportgraphics(gca, '1904.png');