## Wednesday, March 26, 2014

### Matlab Plot numbers on map

Wanted to plot a map with the numbers. I am adding +1 to the value if any zero observations were made to that station. It would be nice to show columns instead of the circles.

% plot the station numbers on the map

close all
latlim = [40, 45];
lonlim = [-80,-73];

figure

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'Facecolor', [1 1 0.9]);

for jj = 1:41
jj

%      textm(st_lat(jj),st_lon(jj), num2str(jj));

plotm(st_lat(jj),st_lon(jj), jj,'o', 'MarkerEdgeColor','k','MarkerSize',dataavailable(jj)+1);

hold on
textm(st_lat(jj),st_lon(jj),  sprintf('%.1d', jj))

end

## Tuesday, March 18, 2014

### Number of data available for PM

There are some stations with zero data availability.
% How many stations are available?
clc
imagesc(Data(:,1:3));figure(gcf);
colorbar
map2 = colormap; map2( 1, : ) = 1; colormap(map2)
cmap = colormap; % cmap nicely puts colormap into 3 col data
xlabel('Year: 2006 2007 2008')
ylabel('Station Number')
title('Number of Data Available')
xlim([0.5 3.5])
% set(gca,'xtick',[])
% set(gca,'xticklabel',[])
set(gca,'XLim',[.5 3.5]);% This automatically sets the
% XLimMode to manual.
% Set XTick so that only the integer values that
% range from 0.5 - 12.5 are used.
set(gca,'XTick',[1:3])  % This automatically sets
years = ['2006'; '2007'; '2008'];
set(gca,'XTickLabel',years)

asdf

## Monday, March 17, 2014

### Script to stitch the images

Call the MODISdaily Script to stitch the images
clc
clear all
close all
cd /home/nabin/dineof-3.0/DainelTry/Daniel_Codes

%Julian date and corresponding calendar date
N_JT = N(:,1);
N_mm = N(:,2);
N_dd = N(:,3);

disp('working on aqua and terra for 2008')

%% variables
clc

%Grid limits
Lat_min = 35;
Lat_max = 46;
Lon_min = -85;
Lon_max = -69;

%vector of coordinates
latv = Lat_min:0.1:Lat_max;
lonv = Lon_min:0.1:Lon_max;

%creates box grid
for i=1:length(latv)
for j=1:length(lonv)
box_lat(i,j) = latv(i);
box_lon(i,j) = lonv(j);
end
end
clear latv lonv

saveAODboxIDW = [];
saveAODday = [];

% end
for daynum =1:366

[box_lat box_lon AODboxIDW AODday] = daynumMODIS(daynum, N, box_lat, box_lon);

saveAODboxIDW(:,:,daynum) = AODboxIDW;
saveAODday(:,:,daynum)= AODday;

end

save box_Lat_lon_IDW_AODav2008.mat box_lat box_lon saveAODboxIDW saveAODday

### Stitch the MODIS images

: called by another script
function [box_lat box_lon AODboxIDW AODday] = daynumMODIS(daynum, N, box_lat, box_lon)
%  daynum = 102

scale_factor = 0.001;
fig_num = 1;
n_file = 0;

N_JT = N(:,1);
N_mm = N(:,2);
N_dd = N(:,3);

%For plots
latlim =([35 46]);
lonlim =([-85 -69]);

Lat_min = 35;
Lat_max = 46;
Lon_min = -85;
Lon_max = -69;

% aTarget =  '/media/Seagate Backup Plus Drive/modis/AQUA_2006/';
% tTarget =  '/media/Seagate Backup Plus Drive/modis/TERRA_2006/';

aTarget =  '/host/data/modisDATA/AQUA_2008/';
tTarget =  '/host/data/modisDATA/TERRA_2008/';

Target = {aTarget, tTarget}

n_file = 0;

for satellite = 1:2
satellite
dir_path = Target{satellite};
dir_info = dir(dir_path);

iwant = [];
for jj = 3: length(dir_info)
if (strcmp(dir_info(jj).name(15:17),sprintf('%03.f',daynum)))
iwant = [iwant jj];
end

end

if isempty(iwant)
disp(['Empty: ', num2str(daynum)])
% do nothing
else

%For each file in the directory
for nf=iwant

file_name = dir_info(nf).name;                      %name of file
file_size = dir_info(nf).bytes;                     %size of file

%if correct file and not empty
if (length(file_name)==44 && file_size>1)
full_name = [dir_path file_name];               %full file name
file_yy = str2num(file_name(11:14));            %gets year of the file
file_nd = str2num(file_name(15:17));            %gets number of day of the file
file_hh = str2num(file_name(19:20));            %UTC hour
file_mn = str2num(file_name(21:22));            %UTC minutes

%finds month and day of file
ind  = find(N_JT==file_nd);
file_mm = N_mm(ind);                            %file month
file_dd = N_dd(ind);                            %file day
clear ind

Corr_Opt   = double(hdfread(full_name,'Corrected_Optical_Depth_Land'));  %Corrected Optical Thickness at 0.47, 0.55, and 0.66 �m

%gets 550nm wavelength values only
AOD550_land = squeeze(Corr_Opt(2,:,:));         %collects the data in the second wavelength 550.

%%%%%%%%%%% need to see figure

%Final image is the original average AOD
%      figure(fig_num)
%      ax = worldmap(latlim, lonlim);
%     states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
%     geoshow(ax, states, 'FaceColor', [1 1 1]);
%     surfacem(Lat,Lon,AOD550_land);
%     fig_num = fig_num+1

sprintf('MODIS AOD %i-%i-%i',file_mm,file_dd,file_yy)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Find bounding box
[indrv, indcv]  = find((Lat<=Lat_max & Lat>=Lat_min) & (Lon>=Lon_min & Lon<=Lon_max));

%initializes matrix that holds AOD values
[row, col] = size(box_lat);
box_AOD = NaN(row, col);

%if intersection between bounding box and file coordinates is found
for Q=1:length(indrv)
%transforms file coordinates (larger precision) into grid coordinates (0.1 degree precision)
cur_lat = round(Lat(indrv(Q),indcv(Q))*10)/10;  %gets latitude
cur_lon = round(Lon(indrv(Q),indcv(Q))*10)/10;  %gets longitude

%nearby points for IDW calculation
[x,y] = find((abs(cur_lat-Lat)<dist_lim) & (abs(cur_lon-Lon)<dist_lim));

%corrects negative values and counts valid AOD readings
valid = 0;
for np=1:length(x)
if(AOD550_land(x(np),y(np))<0)
AOD550_land(x(np),y(np)) = NaN;
else
valid = valid + 1;
end
end

top = 0;
bottom = 0;

%At least 25% of the points found need to be known
if valid>=(length(x)*(.25))
%per each point
for k=1:length(x)
%distance weigthed
c = Lat(x(k),y(k));
d = Lon(x(k),y(k));
point_AOD = scale_factor*AOD550_land(x(k),y(k));

%calculates distance
sum = ((cur_lat - c)^2) + ((cur_lon - d)^2);
distance = sqrt(sum);

%distance weighted average
top = top + (1/distance^2)*point_AOD;
bottom = bottom + (1/distance^2);
end

%finds corresponding coordinates in bounding box grid and
%appends corrected AOD
[cur_row, cur_col] = find(box_lat==cur_lat & box_lon==cur_lon);
box_AOD(cur_row,cur_col) = (top/bottom);
end
end

%         %checks if there are AOD values to display in image.
%         cur_max = max(max(box_AOD));
%         if(isnan(cur_max)==0)
%             fig_title = sprintf('MODIS AOD %i-%i-%i-%i:%i UTC',file_mm,file_dd,file_yy,file_hh,file_mn);
%             figure(fig_num)
%             orient landscape;
%             ax = worldmap(latlim, lonlim);
%             states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
%             geoshow(ax, states, 'FaceColor', [1 1 1]);
%             surfacem(box_lat,box_lon,box_AOD);
%             h1=gca;
%             set(h1,'CLimMode','Manual','CLim',[0, 1]);
%             h2=colorbar('vert','FontSize',12);
%             title(fig_title,'FontSize',14);
%             geoshow(lats,lons,'DisplayType','Line', 'Color', [0 0 0])
%             fig_num = fig_num + 1;
%         end
n_file = n_file+1;
file_AOD(:,:,n_file) = box_AOD;
end
end % end of file list

end
n_file
end % end of satellite number

if  isempty(iwant) % need to pass empty box to the image output
disp(['Empty: ', num2str(daynum)])

[row, col] = size(box_lat);
%                 box_AOD = NaN(row, col);

AODboxIDW = NaN(row, col);

AODday = NaN(row, col);

else

Ave_AOD = nanmean(file_AOD,3);
%     clear file_AOD

%Final image is the original average AOD
fig_title = sprintf('MODIS AOD %i-%i-%i',file_mm,file_dd,file_yy);
figure(1)
orient landscape;
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(box_lat,box_lon,Ave_AOD);
h1=gca;
set(h1,'CLimMode','Manual','CLim',[0, 1]);
h2=colorbar('vert','FontSize',12);
title(fig_title,'FontSize',14);
%     geoshow(lats,lons,'DisplayType','Line', 'Color', [0 0 0])
fig_num = fig_num + 1;

%Applies IDW with different radius to improve spatial coverage
old_AOD_image = Ave_AOD;
new_AOD_image = Ave_AOD;
dist_array =[0.3, 0.5, 1];

[IDWrow,IDWcol] = size(Ave_AOD);
%Per each distance radius in dist_array
for z=1:length(dist_array)
dist_lim = dist_array(z);   %current distance radius
for P=1:IDWrow
for L=1:IDWcol
IDW_AOD = old_AOD_image(P,L);

%if nan value
if((isnan(IDW_AOD))==1)
%gets coordinates
cur_lat = box_lat(P,L);
cur_lon = box_lon(P,L);

%gets surrounding points
[x,y] = find((abs(cur_lat-box_lat)<dist_lim) & (abs(cur_lon-box_lon)<dist_lim));

%if there are points around
if(~isempty(x))

top = 0;
bottom = 0;

%if you have at least 25% valid data
total_points = length(x);       %number of points found around
j=1;
for i=1:total_points
cur_x = x(i);
cur_y = y(i);
if (isnan(old_AOD_image(cur_x,cur_y))==0)
AOD_around(j) = old_AOD_image(cur_x,cur_y);
lat_around(j) = box_lat(cur_x,cur_y);
lon_around(j) = box_lon(cur_x,cur_y);
j = j+1;
end

end

%if at least 25% of points are valid
if ((j-1)>=(total_points*0.25))
%per each point
for k=1:length(AOD_around)

%distance weigthed
sum = ((cur_lat - lat_around(k))^2) + ((cur_lon - lon_around(k))^2);
distance = sqrt(sum);
top = top + (1/distance^2)*AOD_around(k);
bottom = bottom + (1/distance^2);
end

%finds corresponding coordinates in bounding box grid and
%appends corrected AOD
new_AOD_image(P,L) = (top/bottom);
end
clear AOD_around lat_around lon_around

end
clear x y
end
end
end

%         fig_title = sprintf('MODIS AOD %i-%i-%i (IDW %0.1f)',file_mm,file_dd,file_yy,dist_lim);
%         figure(fig_num)
%         orient landscape;
%         ax = worldmap(latlim, lonlim);
%         states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
%         geoshow(ax, states, 'FaceColor', [1 1 1]);
%         hold on
%         surfacem(box_lat,box_lon,new_AOD_image);
%         h1=gca;
%         set(h1,'CLimMode','Manual','CLim',[0, 1]);
%         h2=colorbar('vert','FontSize',12);
%         title(fig_title,'FontSize',14);
% %         geoshow(lats,lons,'DisplayType','Line', 'Color', [0 0 0])
%         fig_num = fig_num + 1;

old_AOD_image = new_AOD_image;
end

AOD_value = new_AOD_image;

fig_title = sprintf('MODIS AOD %i-%i-%i',file_mm,file_dd,file_yy)
figure(2)
orient landscape;
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
hold on
surfacem(box_lat,box_lon,new_AOD_image);
h1=gca;
set(h1,'CLimMode','Manual','CLim',[0, 1]);
h2=colorbar('vert','FontSize',12);
title(fig_title,'FontSize',14);
%     geoshow(lats,lons,'DisplayType','Line', 'Color', [0 0 0])
%     fig_num = fig_num + 1;

AODboxIDW = new_AOD_image;

AODday = Ave_AOD;
end

clear Ave_AOD new_AOD_image
return
%  close all

## Friday, March 7, 2014

### Make the movie matlab

Making movie in matlab is not hard. The following script was used to create the movie...

plotTEST = 1

if plotTEST

for kk = 1:365

%For plots
latlim =([35 46]);
lonlim =([-85 -69]);

%     fig_title = sprintf('MODIS AOD %i-%i-%i',file_mm,file_dd,file_yy);
h1 = figure(1)
%     orient landscape;
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(box_lat,box_lon,saveAODboxIDW(:,:,kk));
h1=gca;
set(h1,'CLimMode','Manual','CLim',[0, 1]);
h2=colorbar('vert','FontSize',12);
title(['2006 Day: ', num2str(kk)],'FontSize',14);

M(kk) = getframe(h1);

h2 = figure(2)
%     orient landscape;
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(box_lat,box_lon,saveAODday(:,:,kk));
h1=gca;
set(h1,'CLimMode','Manual','CLim',[0, 1]);
h2=colorbar('vert','FontSize',12);
title(['2006 Day: ', num2str(kk)],'FontSize',14);

MM(kk) = getframe(h2);

%     geoshow(lats,lons,'DisplayType','Line', 'Color', [0 0 0])
end

% Output the movie as an avi file
movie2avi(M,'AOD_IDWMovie07.avi', 'compression','None', 'fps',2);

movie2avi(MM,'AOD_Movie07.avi', 'compression','None', 'fps',2);

end

### Check the hdf file for corruption

So, some of my hdf files were corrupted

% this script checks the sanity of hdf files.
aTarget =  '/host/data/modisDATA/AQUA_2007/';
tTarget =  '/host/data/modisDATA/TERRA_2007/';

Target = {aTarget, tTarget}

for daynum = 1:365
daynum

for satellite = 1:2
satellite
dir_path = Target{satellite};
dir_info = dir(dir_path);

iwant = [];
for jj = 3: length(dir_info)

if (strcmp(dir_info(jj).name(15:17),sprintf('%03.f',daynum)))
iwant = [iwant jj];
end

end

%For each file in the directory
for nf=iwant

file_name = dir_info(nf).name;                      %name of file
file_size = dir_info(nf).bytes;                     %size of file

full_name = [dir_path file_name];

try

Corr_Opt   = double(hdfread(full_name,'Corrected_Optical_Depth_Land'));  %Corrected Optical Thickness at 0.47, 0.55, and 0.66 �m
catch
disp(['delete ', full_name])

end

end
end

end