clear
clc
ImgNameQC=input('输入LSTQC图像文件名:','s');%写QC的图像文件名
%img=imread('fire2geo.tif');%提取图像
[QC,geo]=geotiffread(ImgNameQC);%提取图像
info=geotiffinfo(ImgNameQC);
OutImgName=strcat(extractBefore(ImgNameQC,"."),"LST.tif");
ImgNamelst=input('输入LST图像文件名:','s');%写lst的图像文件名
[lstdata,geo]=geotiffread(ImgNamelst);%提取图像
%%制作QC文件掩膜
[x, y] = size(QC);%获取矩阵的行列
QC_dec = dec2bin(QC,8);%十进制转二进制%按列转换
QC_78 = QC_dec(:,7:8);%LST的Mandatory QA是0bit和1bit,也就是二进制八位数的最右两位,即第七和第八位
QC78 = str2num(QC_78);%字符串转数值，判断是否符合要求
data0 = reshape(QC78,x,y);
data0(data0==10)=-999;%将质量差的像元设为Nodata
data0(data0==11)=-999;
data0(data0==00)=1;%将可用的像元设为1
data0(data0==01)=1;
data1 = double(lstdata);%为了保存小数点，double可以保存15位小数%%lstdata是之前读进去的lst数据，此处未展示完整代码
qc = double(data0);
data2 = data1*0.02-273.15;%LST的系数为0.02，单位为开尔文，根据需要换算成摄氏度
data3 = data2.*qc;%利用质量掩膜来提取质量好的像元
geotiffwrite(OutImgName, data3, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
MOD11A2 11个时间LST取平均：
( ((b1 eq 272876.85)*0+(b1 lt 272876.85)*b1)+((b2 eq 272876.85)*0+(b2 lt 272876.85)*b2)+((b3 eq 272876.85)*0+(b3 lt 272876.85)*b3)+((b4 eq 272876.85)*0+(b4 lt 272876.85)*b4)+((b5 eq 272876.85)*0+(b5 lt 272876.85)*b5)+((b6 eq 272876.85)*0+(b6 lt 272876.85)*b6)+((b7 eq 272876.85)*0+(b7 lt 272876.85)*b7)+((b8 eq 272876.85)*0+(b8 lt 272876.85)*b8)+((b9 eq 272876.85)*0+(b9 lt 272876.85)*b9)+((b10 eq 272876.85)*0+(b10 lt 272876.85)*b10)+((b11 eq 272876.85)*0+(b11 lt 272876.85)*b11) )/( ((b1 lt 272876.85)+(b2 lt 272876.85)+(b3 lt 272876.85)+(b4 lt 272876.85)+(b5 lt 272876.85)+(b6 lt 272876.85)+(b7 lt 272876.85)+(b8 lt 272876.85)+(b9 lt 272876.85)+(b10 lt 272876.85)+(b11 lt 272876.85))>0 )