: 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;
dist_lim = .1; %initial radius
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
%reads file coordinates
Lat = double(hdfread(full_name,'Latitude')); %latitudes
Lon = double(hdfread(full_name,'Longitude')); %longitudes
%reads file variables
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);
% load Variables\lons
% load Variables\lats
% 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);
% load lons
% load lats
% 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);
% % load lons
% % load lats
% % 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);
% load lons
% load lats
% 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
No comments:
Post a Comment