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')