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