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

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);
   
   
   
%     load lons
%     load lats
%     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