iteration limit reached in robustfit matlab

Warning: Iteration limit reached.
> In stats\private\statrobustfit at 76
  In robustfit at 106

So, I tracked down the robustfit function. In MATLAB 2013a, it is located at:
C:\MATLAB\R2013a\toolbox\stats\stats\private
as statrobustfit.
Around line 70 was the iteration limit set to 50. Changed that to 250. So far my data does not seem to get any improvements by this; but it is here if it could be useful.

iter = 0;
iterlim = 250; % modified from 50
wxrank = xrank;    % rank of weighted version of x
while((iter==0) || any(abs(b-b0) > D*max(abs(b),abs(b0))))
   iter = iter+1;
   if (iter>iterlim)
      warning(message('stats:statrobustfit:IterationLimit'));

creating csv with header text in matlab

The csvwrite function in matlab does not work well with the text. There is a table feature in the new matlab version. However, here is a quick way to create csv with text header in matlab

csvData = [Elev TPW SkinT SurfacePr dayyes' ProfType Uyear Umo Uhr latitude longitude];

header=['Elev,TPW,SkinT,SurfacePr,dayyes,ProfType, year, mo,HH, latitude,longitude'];
outid = fopen('SeeBor.csv', 'w+');
fprintf(outid, '%s', header);
fclose(outid);
dlmwrite ('SeeBor.csv',csvData,'roffset',1,'-append')

Latitude Longitude and Data on Map @matlab

I wanted to have a nice map of the global dataset given the lat/lon and the dataset.
This code is being put here so that it could be used later.

boxmap = 0
% close all

latlim = [-90, 90];
lonlim = [-180, 180];
load coast
figure

if boxmap
    axesm('MapProjection','lambcyln','grid','on', ...
        'MapLatLimit',latlim,'MapLonLimit',lonlim,...
        'MLineLocation',15, 'PLineLocation',15)
else
    ax = worldmap(latlim, lonlim);
 
end
geoshow(lat, long,'Color', 'black' )
set(gcf,'Color','white')

for jj = 1:length(cimms)
   
    [SunRiseSet,Day,Dec,Alt,Azm,Rad] = suncycle( latitude(jj) , longitude(jj) , [Uyear(jj) Umo(jj) Uday(jj)] , 2880 );
 
    if SunRiseSet(2)<SunRiseSet(1)
        SunRiseSet(2) = SunRiseSet(2)+24;
    end
 
    if Uhr(jj)>SunRiseSet(1) & Uhr(jj)<SunRiseSet(2)
        plotm(latitude(jj),longitude(jj), jj,'o', 'MarkerEdgeColor','r','MarkerSize',8);
     
        dayyes(jj) = 1;
    else
        plotm(latitude(jj),longitude(jj), jj,'o', 'MarkerEdgeColor','b','MarkerSize',8);
     
        dayyes(jj) = 0;
     
    end
 
    hold on
 
    if rem(jj,1000)==0
        textm(latitude(jj),longitude(jj),  sprintf('%.1d', jj), 'color', 'g')
     
        pause(0.5)
     
    end
 
 
end

Byobu shortcuts for screen/ssh

Byobu is a lifesaver for multiple jobs to be done on the command line via ssh client.
Currently I am using it in putty.
You can create, close and detach... and come right back to the place you left off!

Here are some of the handy commands

Start with byobu

Ctrl-a, c: create a new window
F8: rename
Ctrl-a C-a - Go back to the previous window

Ctrl-a <0-9> - Switch to screen # 0-9.
If you want to go to screen # 10,
Ctrl-a ' - Enter a screen number to switch to two digit number

Ctrl-a " - View a list of the current screens, which will allow you to select one from the list

Ctrl-a d - Detach the whole screen session and fork to the background. Very useful for remote sessions you want to leave open. The command "screen -r" will resume your screen session.

Ctrl-a, k: kill the screen. Usually detach is good; but sometimes you need to kill the screen/process.

If you hit F9, you can set options to start with byobu screen.

Ctrl-a, n or p does next and p

Map data matlab plot

Plotting Map in Matlab with the colorbar and NaNs set as transparent white.

Just call it with the lat,lon and Z (same size), and do caxis of your desire at the end. Voila!
Save it.

function hhf = mapme(Latitude, Longitude, Z)
hhf = figure
load coast
latlim=[floor(min(min(Latitude))),ceil(max(max(Latitude)))];
lonlim=[floor(min(min(Longitude))),ceil(max(max(Longitude)))];
ax = worldmap(latlim, lonlim);
pcolorm(Latitude, Longitude, Z);
geoshow(lat, long,'Color', 'black' )
colormap;
set(gcf,'Color','white')
% map2 = colormap; map2( 64, : ) = 1; colormap(map2);
colorbar

Map/ Figure with white color bar assignment to the lowest value

Sometimes we need the figures or plots showing the maps. The oceans/missing values are indicated in the data as NaN or -9999. This code will plot the data with the lowest value being shown as white color.

A = rand(55);

figure;
imagesc(A)
colormap(jet(33))
map2 = colormap; map2( 1, : ) = 1; colormap(map2)
cmap = colormap; % cmap nicely puts colormap into 3 col data
colorbar

Open Chrome or Firefox from command using matlab

This simple script in Windows 7/8 worked fine.

Sometime I need to open multiple urlwant string
urlwant = ["nabinkm.com"]
Voila!

   system(['start chrome ',urlwant])

Turn off the spotlight and turn it on

Turn off the spotlight and turn it on using these commands in terminal.

sudo mdutil -i off

sudo mdutil -E /

Monitor Temperature in Mac using bash script command line

The code is simple thanks to the temperature monitor app
Just do ./runTapp.sh &
in the terminal, and you are good to go!

#!/bin/sh
cd /Applications/TemperatureMonitor.app/Contents/MacOS/
./tempmonitor -th > ~/Desktop/CPUtempData.csv
while true
do
    ./tempmonitor -tv >> ~/Desktop/CPUtempData.csv
    /bin/sleep 60
done

Reference:
http://www.bresink.com/osx/216202/Docs-en/commandline.html

Map Subsection USA States using matlab

I had the regional PM2.5 map, but I wanted to have it by states. So, came up with this code which can generate state level maps. If you use it, please let me know. Current example uses the solar radiation from GCM.

Highlight one, highlight all!

There is also R version available, at:

% Create the northeast mask
% Created by NM
% August 4, 2014
%
% http://www.mathworks.com/help/map/_f7-10852.html#f7-9156
% Use shaperead to get the patch data for the boundary:
close all
% http://www.mathworks.com/help/map/understanding-raster-geodata.html#f20-12849
pcs = {'New York', 'New Jersey','Massachusetts', 'Vermont','Connecticut' ,...
     'North Carolina','Ohio','Maine', 'New Hampshire', 'Rhode Island', ...
    'Maryland', 'Delaware', 'Virginia','West Virginia','District of Columbia' };
% pcs = {     'Pennsylvania'  };

cUS = shaperead('usastatelo.shp',...
    'UseGeoCoords', true,...
    'Selector',{@(name)any(strcmpi(name,pcs),2), 'Name'});

inLat = [cUS.Lat];
inLon = [cUS.Lon];
%%
% nys = shaperead('usastatehi.shp',...
%     'UseGeoCoords', true,...
%     'Selector', {@(name)strcmpi('New York',name), 'Name'});
% inLat = nys.Lat;
% inLon = nys.Lon;

% Set the grid density to be 40 cells per degree, and use vec2mtx to rasterize the boundary and generate a referencing vector for it:
gridDensity = 100;
[inGrid, inRefVec] = vec2mtx(inLat, inLon, gridDensity);
whos

% set the limits
figure(1)
axesm eqdcyl
meshm(inGrid, inRefVec)
colormap jet(4)

[latlim, lonlim] = limitm(inGrid, inRefVec);
setm(gca, 'Flatlimit', latlim, 'FlonLimit', lonlim)
tightmap

% To fill (recode) the interior of Indiana, you need a seed point
% (which must be identified by row and column) and
% a seed value (to be allocated to all cells within the polygon).
% Select the middle row and column of the grid and choose an index value of 3
% to identify the territory when calling encodem to generate a new grid:
inPt = round([size(inGrid)/2, 3]);
inGrid3 = encodem(inGrid, inPt,1);
% redraw the map using the filled grid:
inGrid3= 1-inGrid3; % need for multiple states

meshm(inGrid3, inRefVec)


%% Now try the Data
figure(2);
udGrid = flipud(inGrid3);
imagesc(udGrid)

[m,n] = size(udGrid);


% latlim=[36 46];
%     lonlim=[-82 -69];
     [Mlat, Mlon] = ndgrid(linspace(min(inLat),max(inLat),m), linspace(min(inLon),max(inLon),n));   %Changed from 35 46, -85 -69
%     newCF = griddata(latmatf,lonmatf,flipud(data),Dlat, Dlon);
 
    figure(3)
    ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
%     geoshow(ax, states, 'FaceColor', [1 1 1]);
    surfacem(Mlat,Mlon,(inGrid3));
 
 
    %% Now try the data
  load('C:\Users\CREST\Copy\CCNY-14\downscale-now\Regression-2006-1.mat')

    %     newCF = griddata(latmatf,lonmatf,flipud(data),Dlat, Dlon);
figure(4)
  ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Dlat,Dlon,isimipAve);
 
    %%
        newCF = griddata(Dlat,Dlon,isimipAve,Mlat, Mlon);
     
        inGrid3c = inGrid3;
        inGrid3c(inGrid3c<=0)= NaN;
        inGrid3c(inGrid3c>0) = 1;
        plotCF = inGrid3c.*newCF;
       
        figure(5)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,plotCF);
 
 
%% For Penn because When added Penn, it "inverts". So adding two with nanmean
pcs = {   'New York', 'Pennsylvania',  'Massachusetts', 'Vermont','Connecticut' ,...
     'North Carolina','Ohio','Maine', 'New Hampshire', 'Rhode Island', ...
    'Maryland',  'Virginia','West Virginia','District of Columbia' };

cUS = shaperead('usastatelo.shp',...
    'UseGeoCoords', true,...
    'Selector',{@(name)any(strcmpi(name,pcs),2), 'Name'});

inLat = [cUS.Lat];
inLon = [cUS.Lon];

% nys = shaperead('usastatehi.shp',...
%     'UseGeoCoords', true,...
%     'Selector', {@(name)strcmpi('New York',name), 'Name'});
% inLat = nys.Lat;
% inLon = nys.Lon;

% Set the grid density to be 40 cells per degree, and use vec2mtx to rasterize the boundary and generate a referencing vector for it:
gridDensity = 100;
[inGrid, inRefVec] = vec2mtx(inLat, inLon, gridDensity);
whos

% set the limits
% figure(1)
% axesm eqdcyl
% meshm(inGrid, inRefVec)
% colormap jet(4)

% [latlim, lonlim] = limitm(inGrid, inRefVec);
% setm(gca, 'Flatlimit', latlim, 'FlonLimit', lonlim)
% tightmap

% To fill (recode) the interior of Indiana, you need a seed point
% (which must be identified by row and column) and
% a seed value (to be allocated to all cells within the polygon).
% Select the middle row and column of the grid and choose an index value of 3
% to identify the territory when calling encodem to generate a new grid:
inPt = round([size(inGrid)/2, 3]);
inGrid3 = encodem(inGrid, inPt,1);
% redraw the map using the filled grid:
inGrid3=   inGrid3; % need for multiple states

meshm(inGrid3, inRefVec)
%

  inGrid3c = inGrid3;
        inGrid3c(inGrid3c<=0)= NaN;
        inGrid3c(inGrid3c>0) = 1;
        plotCF2 = inGrid3c.*newCF;
  figure(6)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,plotCF2);
 
    plotCF3d(:,:,1) = plotCF;
    plotCF3d(:,:,2) = plotCF2;
 
    iplotCF = nanmean(plotCF3d,3);
 
      figure(7)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,iplotCF);  

Selected Land class data plot

I wanted to make the land class plot

for four selected datasets.

load('Land_Class_by_Percentage_3646N_8269W.mat')
[lonmat,latmat]=meshgrid(lonv_landclass,latv_landclass);
iwant = find(lonmat <-82.05 | lonmat >-69 | latmat >46.05 | latmat <36);
iLand_Class(iwant) = [];

iA = reshape(iLand_Class, 201,261);
lonv_landclass(lonv_landclass<-82.05 | lonv_landclass>-69) = [];
latv_landclass(latv_landclass>46.05 | latv_landclass<36.03) = [];
[lonmat,latmat]=meshgrid(lonv_landclass,latv_landclass);

%%
latlim=[36 46];
lonlim=[-82 -69];
% figure to see
figure
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(latmat,lonmat,iA);
title('Land Classification')
colormap;
map2 = colormap; map2( 1, : ) = 1; colormap(map2);
colorbar

iwant2 = find(iA ~=5 & iA ~=13 & iA ~= 14& iA ~= 15);
iA2 =  (iA);
iA2(iwant2) = NaN;
iA2 =(reshape(iA2, size(iA)));

iA2(iA2==5) = 1;
iA2(iA2==13) = 2;
iA2(iA2==14) = 3;
iA2(iA2==15) = 4;
% figure;
% imagesc(iA2)
% get current colormap
map = colormap;
%  adjust for number of colors you want
rows = uint16(linspace(1, size(map,1), 4)) ;
map = map(rows, :);

%  apply   colormap
colormap(map);

labels = {'broadleaf forest','croplands','urban', 'vegetation'};
lcolorbar(labels,'fontweight','bold');

%%
latlim=[36 46];
lonlim=[-82 -69];
% figure to see
figure
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(latmat,lonmat,iA2);
title('Land Classification')
%# get current colormap
map = colormap;

%# adjust for number of colors you want
rows = uint16(linspace(1, size(map,1), 4)) ;
map = map(rows, :);
%# and apply the new colormap
colormap(map);
% colorbar
% caxis([1 4])
labels = {'broadleaf forest','croplands','urban', 'vegetation'};
lcolorbar(labels,'fontweight','bold');

%%
% test to display the land classes
for jj = 1:length(classification)
   
    disp([num2str(jj), ':', classification{jj}])
end

Land Surface Temperature (K) map from MODIS

Just a plot showing LST for Jan 2006. (Working with Mark) Daytime and Nighttime LST.

fullfn = '...\MODIS\2006\MYD11C3.A2006001.005.2008107201018.hdf'
QC_Day = hdfread(fullfn,  '/MODIS_MONTHLY_0.05DEG_CMG_LST/Data Fields/QC_Day', 'Index', {[1  1],[1  1],[240  250]});
QC_Night = hdfread(fullfn, 'MODIS_MONTHLY_0.05DEG_CMG_LST', 'Fields', 'QC_Night');
 scale = 0.02;
% LST_Day_CMG = (hdfread(fullfn,  '/MODIS_MONTHLY_0.05DEG_CMG_LST/Data Fields/LST_Day_CMG', 'Index', {[860  1950],[1  1],[240  250]}));
LST_Day_CMG = (hdfread(fullfn,  '/MODIS_MONTHLY_0.05DEG_CMG_LST/Data Fields/LST_Day_CMG', 'Index', {[1  1],[1  1],[3600/2 7200/2]}));
LST_Day_CMG(LST_Day_CMG>65535) = NaN;
LST_Day_CMG(LST_Day_CMG<7500) = NaN;
LST_Day_CMG = scale *double(LST_Day_CMG);
LST_Day_CMG(LST_Day_CMG==0) = NaN;
% Big Endian
Bit0=bitget(QC_Day,1);% 1
Bit1=bitget(QC_Day,2);% 1
Bit2=bitget(QC_Day,3);% 1

Bit4=bitget(QC_Day,5);% 1
Bit5=bitget(QC_Day,6);% 1
Bit6=bitget(QC_Day,7);% 1

% find high quality flag
% In our case, we want Bit 0 and Bit 1 to be 0, this means good quality
ind=find(Bit0==1&Bit1==1);%&Bit2==1&Bit4==1&Bit5==1&Bit6==1);
LST_Day_CMG(ind) = NaN;


LST_Night_CMG = (hdfread(fullfn, '/MODIS_MONTHLY_0.05DEG_CMG_LST/Data Fields/LST_Night_CMG', 'Index', {[1  1],[1  1],[3600/2 7200/2]} )) ;
LST_Night_CMG(LST_Night_CMG>65535) = NaN;
LST_Night_CMG(LST_Night_CMG<7500) = NaN;
LST_Night_CMG = scale *double(LST_Night_CMG);
% Big Endian
Bit0=bitget(QC_Day,1);% 1
Bit1=bitget(QC_Day,2);% 1
Bit2=bitget(QC_Day,3);% 1

Bit4=bitget(QC_Day,5);% 1
Bit5=bitget(QC_Day,6);% 1
Bit6=bitget(QC_Day,7);% 1

% find high quality flag
% In our case, we want Bit 0 and Bit 1 to be 0, this means good quality
ind=find(Bit0==1&Bit1==1);%&Bit2==1&Bit4==1&Bit5==1&Bit6==1);
LST_Night_CMG(ind) = NaN;

%%
Y = linspace(0,90,3600/2);
X = linspace(-180,0,7200/2);

[Daylat, Daylon] = ndgrid(Y,X);

[Dlat, Dlon] = ndgrid(35:0.05:46, -85:0.05:-69);
asdf = rot90(LST_Day_CMG,3)';
daymet  = griddata(Daylat,Daylon,asdf, Dlat,Dlon , 'natural'); % gather all tiles

asdf1 = rot90(LST_Night_CMG,3)';
nightmet =  griddata(Daylat,Daylon,asdf1, Dlat,Dlon , 'natural'); % gather all tiles

%% test for outliers
iout = reshape(daymet(~isnan(daymet)),[],1);
mout = median(iout);
yy = iqr(iout) ;
iwant = find(abs(daymet<(mout - 3*yy)));
daym = daymet;
daym(iwant) = NaN;

figure
latlim =([35 46]);
lonlim =([-80 -69]);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(Dlat, Dlon,daym);
colorbar

%%
nightmet(nightmet==0)=NaN;
iout = reshape(nightmet(~isnan(nightmet)),[],1);
mout = median(iout);
yy = iqr(iout) ;
iwant = find(abs(nightmet<(mout - 3*yy)));
nightm = nightmet;
nightmet(iwant) = NaN;

figure
latlim =([35 46]);
lonlim =([-80 -69]);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(Dlat, Dlon,nightmet);
colorbar

scatterhist with pdf on the side

Plotting scatterhist with the pdf on the side is much informative!

iwant = find(saveMO==6 | saveMO==7 | saveMO ==8 |...
    saveMO==12 | saveMO==1 | saveMO ==2);
iTE = saveTE(iwant);
iPM = savePM(iwant);
figure
h = scatterhist(iTE,iPM,'Group',saveMO(iwant),...
    'Location','SouthEast',...
    'Direction','out','Color','mgbrkc','LineStyle',{'-.','-.','-.','-','-.','-'},...
    'LineWidth',[2,2,2,2,2,2],'Marker','+ods<*','MarkerSize',[4,5,6, 4,5,6]);
xlabel('Temperature (K)')
ylabel('Station PM2.5 (\mugm^{-3})')
title('Winter and Summer Correlation (Tempr and PM2.5)')

Map with inset

So, I had to zoom into NYC area with the counts of data shown for lat/lon points
here is a figure:

NYC Zoomed in to show more details

close all
latlim = [40, 45];
lonlim = [-80,-72];
% latlim = [40, 42];
% lonlim = [-80,-73];
% find number of stations lat lon and number of data
h1 = figure(1);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'Facecolor', [1 1 0.9]);
hold on
map = colormap;
colorvector = map(1:64,:);
pmvect = linspace(0, max(saveI), 64);
for jj = 1:41
    jj
   
    %      textm(st_lat(jj),st_lon(jj), num2str(jj));
    colorb = interp1(pmvect, colorvector, saveI(jj));
%     plotm(Lat1(jj),Lon1(jj), jj,'o', 'MarkerEdgeColor','k','MarkerSize',saveI(jj)+1);
plotm(Lat1(jj),Lon1(jj), saveI(jj),'o', 'MarkerEdgeColor','k','MarkerFaceColor', colorb, 'MarkerSize', 5);
hold on
%     textm(Lat1(jj),Lon1(jj),sprintf('%.1d', jj))
   
end
caxis([0  max(saveI)])
% colorbar
title('The number of Episodic (PM2.5>35\mug/m^{-3}) events in NY')
hold off


h2 = figure(2);
latlim = [40.5, 40.9];
lonlim = [-74.3,-73.5];
ax = worldmap(latlim, lonlim);
set(ax, 'FontSize', 1)

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

colorvector = map(1:64,:);
pmvect = linspace(0, max(saveI), 64);
for jj = 1:41
    jj
    %      textm(st_lat(jj),st_lon(jj), num2str(jj));
    colorb = interp1(pmvect, colorvector, saveI(jj));
%     plotm(Lat1(jj),Lon1(jj), jj,'o', 'MarkerEdgeColor','k','MarkerSize',saveI(jj)+1);
plotm(Lat1(jj),Lon1(jj), saveI(jj),'o', 'MarkerEdgeColor','k','MarkerFaceColor', colorb, 'MarkerSize', 5);
hold on
%     textm(Lat1(jj),Lon1(jj),sprintf('%.1d', jj))
   
end
caxis([0  max(saveI)])
% set(findobj(gcf, 'type','axes'), 'Visible','off')
hold off
% trying to inset the figure
inset_size=0.4;
inset_size=inset_size*.7;
figure;
new_fig=gcf;
main_fig = findobj(h1,'Type','axes');
h_main = copyobj(main_fig,new_fig);
set(h_main,'Position',get(main_fig,'Position'))
colorbar
inset_fig = findobj(h2,'Type','axes');
h_inset = copyobj(inset_fig,new_fig);
ax=get(main_fig,'Position');
% set(h_inset,'Position', [.7*ax(1)+ax(3)-inset_size .5*ax(2)+ax(4)-inset_size inset_size inset_size])
set(h_inset,'Position', [0.3*ax(1) 1.5*ax(2) inset_size*2 inset_size])

Creating Stack Bar Plot

Creating stack plot experiment in Matlab.
First count items in the range, then make distribution stack plot by month

for jj = 1:12
iwant = find(dates(:,2)==jj);
iPBL = []
iPBL = PBL(iwant);
PBL5(jj) = length(iPBL(iPBL<500));
PBL1(jj) = length(iPBL(iPBL<1000 &  iPBL>=500));
PBL15(jj) =length(iPBL(iPBL<1500 &  iPBL>=1000));
PBL2(jj)=length(iPBL(iPBL<2000 &  iPBL>=1500));
PBL25(jj) =length(  iPBL(iPBL>=2000));
end

 % Create a stacked bar chart using the bar function
figure;
bar([1:12], [PBL5; PBL1; PBL15; PBL2; PBL25]', 'stack');

% Add title and axis labels
title('PBL distribution by month');
xlabel('Month');
ylabel('PBL range count');

% Add a legend
legend('<.5', '.5\leqc<1', '1<c\leq1.5', '1.5<c\leq2', '2\leq c');

Applying griddata to MODIS granule

Text here So, the 'v4' option takes a lot of memory, crashes your small computer.
While I can show 'natural', 'linear' (default) and 'nearest' effects on the griddata.

At the least, do not use the nearest option because it seems to have some spurious effect due to re-gridding.

full_name = '/zxcv.hdf' scale = 0.001; %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 = scale*squeeze(Corr_Opt(2,:,:)); %collects the data in the second wavelength 550. fill = -9999*scale AOD550_land(AOD550_land==fill)=NaN; %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; [ilatv, ilonv] = meshgrid(latv,lonv); AODgrid = griddata(Lat, Lon, AOD550_land, ilatv,ilonv, 'natural'); %% figure(1) latlim =([36 39]); lonlim =([-85 -74]); ax = worldmap(latlim, lonlim); states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']); geoshow(ax, states, 'FaceColor', [1 1 1]); surfacem(ilatv,ilonv,AODgrid); colorbar title([ 'natural 2006199'], 'FontSize', 20)

Shift in the latitude in Rsig data

There was a mysterious shift in the latitude, now corrected by the following back-projection. Thanks to Todd!

latnc = ilat(:,:,1);
WGS84AxisRatioSquared = 0.9933056199957391;
latitudeWGS84Radians = latnc.*pi/180;
latitudeSphereRadians = atan( tan( latitudeWGS84Radians )  * WGS84AxisRatioSquared );
newLat =    latitudeSphereRadians *180/pi;

The Latitude projection into this Lambert Conic map seems to have some shift of about 0.2 degree (12 km resolution).

Once the Backprojection is applied, the latitude appears to align nicely.

C05 and C06 side by side

Just comparing the two MODIS AOD products side by side (DT)
 

The improvements on DB from C6 (left) is amazing!

fn = '/host/data/C6modis/aqua/2006A/MYD04_L2.A2006026.1840.006.2014036014118.hdf';
% fn = 'MYD04_L2.A2006026.1840.005.2007129090710.hdf'
Longitude = double(hdfread(fn, 'Longitude'));
Latitude = double(hdfread(fn, 'Latitude'));

scale =  0.001;

% AOD =scale* double(hdfread(fn, 'Corrected_Optical_Depth_Land'));
AOD =scale* double(hdfread(fn, 'Deep_Blue_Aerosol_Optical_Depth_550_Land'));


AOD550 = (AOD); % this is AOD .550 um

fill = -9999*scale
AOD550(AOD550==fill)=NaN;


% figure;
% imagesc(AOD550)
figure
subplot(1,2,1)
latlim =([32.5 42]);
lonlim =([-85 -73]);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(Latitude,Longitude,AOD550);
colorbar
title([ '2006026 C6 DB'])

%%
 fn = '/host/data/modisDATA/AQUA_2006/MYD04_L2.A2006026.1840.005.2007129090710.hdf'
Longitude = double(hdfread(fn, 'Longitude'));
Latitude = double(hdfread(fn, 'Latitude'));

scale =  0.001;

% AOD =scale* double(hdfread(fn, 'Corrected_Optical_Depth_Land'));
AOD =scale* double(hdfread(fn, 'Deep_Blue_Aerosol_Optical_Depth_550_Land'));


AOD550db = (AOD); % this is AOD .550 um

fill = -9999*scale
AOD550db(AOD550db==fill)=NaN;


% figure;
% imagesc(AOD550)
subplot(1,2,2)
 
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(Latitude,Longitude,AOD550db);
colorbar
title(['2006026 C5 DB'])

for input arguments of type 'single' and attributes 'sparse 2d real'

For the problem of the kind:

Undefined function 'plus' for input arguments of type 'single' and attributes 'sparse 2d real'.

I forgot to convert the data into double. do data = double(data); Problem solved!

Tic Toc for time required

Wanted to make grids, and also be informed about the time required to complete:
Implemented tic/toc to count the time... to completion.

clear all
% file location
location = 'E:\scienceFTP\';
fn = 'PBL_2006.nc'
full_fn = [location, fn];
% information1 = ncinfo(full_fn)

lon= double(squeeze(ncread(full_fn, 'longitude')));
lat = double(squeeze(ncread(full_fn, 'latitude')));
PBL = double(squeeze(ncread(full_fn, 'PBL')));
% height = double(squeeze(ncread(full_fn, 'height')));

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

%vector of coordinatesI want
Dlat = Lat_min:0.1:Lat_max;
Dlon = Lon_min:0.1:Lon_max;
[ilat, ilon] = meshgrid(Dlat, Dlon);

% available one
mylat = squeeze(double(lat(:,:,1)));
mylon =  squeeze(double(lon(:,:,1)));
rsdata = double(PBL);
tic
for kk = 1:size(rsdata,3);
 
    % This is isimip original
    rs1 =  (rsdata(:,:,kk));
    % [latmat1f, lonmat1f] = ndgrid(mylat, mylon);
    NEgrid(:,:,kk)  = griddata(mylat, mylon,rs1, ilat,ilon );
 
    if kk ==10
        it = toc;
    end
 
    if rem(kk,100)==20
     
        t  = toc;
        kk
               timewant = it*size(rsdata,3)/10 - t;
        disp(['time remaining = ', num2str(timewant) ]);
    end
 
 
end

save saveMAT/PBL2006.mat NEgrid ilat ilon

% later I can do the following to see the plots
latlim =([25 46]);
lonlim =([-85 -69]);

  kk = 16
 ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(ilat,ilon,NEgrid(:,:,kk));

Northeast region map

Basic code in matlab to create northeast states.

latlim = [ 35  48];
lonlim = [-85 -65];

% want states
ilatlim = [ 37  46];
ilonlim = [-80 -69];


figure
ax = usamap(latlim,lonlim);
 axis off

% axesm('MapProjection','conic')
 % The Lambert Conformal Conic Projection is often used for maps of the conterminous United States.
% Here is the map projection usamap selected:
% getm(gca,'MapProjection')

% ans =
% lambert
% Next, use shaperead to read U.S. state polygon boundaries from the usastatehi shapefile into a geostruct named states:
states = shaperead('usastatehi',...
    'UseGeoCoords', true, 'BoundingBox', [ilonlim', ilatlim']);
% Make a symbolspec to create a political map using the polcmap function:
faceColors = makesymbolspec('Polygon',...
    {'INDEX', [1 numel(states)], ...
    'FaceColor', polcmap(numel(states))});
% Display the filled polygons with geoshow:
geoshow(ax, states, 'SymbolSpec', faceColors)
% Extract the names for states within the window from the geostruct and use textm to plot them at the label points provided by the geostruct:
for k = 1:numel(states)
  labelPointIsWithinLimits =...
    37 < states(k).LabelLat &&...
    latlim(2) > states(k).LabelLat &&...
    -82.5 < states(k).LabelLon &&...
    lonlim(2) > states(k).LabelLon;
  if labelPointIsWithinLimits
    textm(states(k).LabelLat,...
    states(k).LabelLon, states(k).Name, ...
        'HorizontalAlignment', 'center', 'FontSize', 15)
  end
end
 
tightmap

Find cases which Temp>85

data.mat has data with temp

clear all
load('data.mat')


ind = 1
row = 1;
% for jj = 1:length(temp)
Doit = 1;
kk = 1;
while Doit
    if temp(kk)>=85
        count = 1;
     
        saveI(ind) = kk;
        ind = ind+1;
        while temp(kk) >= 85
            saveS(row, count) = kk; % save these for the case.
            saveT(count) = temp(kk);
         
            saveC(kk) = count ;
            count = count+1;
         
            kk = kk+1;
         
        end
        row = row+1;
     
    end
    if kk>=length(temp)
        Doit = 0;
    end
    kk= kk+1
 
end

%% Now we find the guys which lasted more than 8 times > 85
for kk = 1:size(saveS)
wantRows(kk) = sum(saveS(kk,:)>0)
end

findRows = find(wantRows>8);

These = saveS(findRows,:);
% These are the instances in which Temp was greater than 85.

Creating DINEOF movie for AOD 2006-2008

Basic code for getting movie for the tabfilled variable

 load('box_Lat_lon.mat');
  
    for kk = 1:length(tabfilled)
      
        %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,tabfilled(:,:,kk));
        h1=gca;
        set(h1,'CLimMode','Manual','CLim',[0, 1]);
        h2=colorbar('vert','FontSize',12);
        title(['2008 Day: ', num2str(kk)],'FontSize',14);
      
      
        M(kk) = getframe(h1);
       
    end
  
    % Output the movie as an avi file
    movie2avi(M,'AOD_IDWMovie2008r1.avi', 'compression','None', 'fps',2);
  
   

Merging PDF with latex

Use pdfpages package in the packet manager, if necessary. Specially in windows 7 I had to install the pdfpackage.sty using packet manager as an admin.

\documentclass{article}

\usepackage{pdfpages}

\begin{document}

\includepdfmerge{One.pdf,-}

\includepdfmerge{Two,-}

\includepdfmerge{Three_part.pdf, 3-14} % mention range

\end{document}

AOD from MODIS example

Text here

fn = 'MOD04_L2.A2014093.1635.051.2014094021045.hdf';

Longitude = double(hdfread(fn, 'Longitude'));
Latitude = double(hdfread(fn, 'Latitude'));

scale =  0.001;

AOD =scale* double(hdfread(fn, 'Corrected_Optical_Depth_Land'));

AOD550 = squeeze(AOD(2,:,:)); % this is AOD .550 um

fill = -9999*scale
AOD550(AOD550==fill)=NaN;


% figure;
% imagesc(AOD550)
figure
latlim =([25 46]);
lonlim =([-100 -69]);

ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(Latitude,Longitude,AOD550);
colorbar
title(fn(11:17))

saveas(gcf, [fn(11:17), '.png'], 'png')

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

load NYstation_data.mat

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

Number of data available for PM

There are some stations with zero data availability.

% How many stations are available?
load('numpoints.mat')
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

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 = xlsread('/home/nabin/dineof-3.0/DainelTry/Daniel_Codes/Day_correspondence.xlsx',1);
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;
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