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)