0

我已经下载了一些 Landsat 8 图像，我想在 MATLAB 中从这些图像中计算 LST（地表温度）。

然而，有些事情似乎不太对劲。例如，当我检查冰川表面的计算温度时，显示的值超过 54 摄氏度。我认为这是不可能的，因为融化的冰面当然不能高于 0 度。谁能帮我检查哪里出了问题？我还没有找到错误。

我使用的代码可以在图像下找到。我使用的公式位于：https ://www.usgs.gov/media/files/landsat-8-data-users-handbook 。我已在以下位置上传了相应的 geotiff 图像：https ://drive.google.com/file/d/1CXxt7hHWNbXujt6rMQPgGw103Zo6vbIc/view?usp=sharing 。

谢谢！

在此处输入图像描述

clear all
close all

B4 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_SR_B4.tif'));
B5 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_SR_B5.tif'));
B10 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_ST_B10.tif'));

%%

% Search for metafile and extract parameters

directory = pwd;                       % Full path of the directory to be searched in
filesAndFolders = dir('*MTL.txt');     % Returns all the files and folders in the directory
metadata_file = filesAndFolders.name;

fid  = fopen(metadata_file,'r');
text = textscan(fid,'%s','Delimiter','');
text = text{1};

% Band-specific multiplicative rescaling factor ML for band 10

idx = find(~cellfun('isempty',strfind(text,'RADIANCE_MULT_BAND_10')));
ML = string(text(idx));
ML = extractAfter(ML,"= ");
ML = str2double(ML);

% Band-specific additive rescaling factor AL from the metadata

idx = find(~cellfun('isempty',strfind(text,'RADIANCE_ADD_BAND_10')));
AL = string(text(idx));
AL = extractAfter(AL,"= ");
AL = str2double(AL);

% Band-specific thermal conversion constant from the metadata K1 from the metadata

idx = find(~cellfun('isempty',strfind(text,'K1_CONSTANT_BAND_10')));
K1 = string(text(idx));
K1 = extractAfter(K1,"= ");
K1 = str2double(K1);

% Band-specific thermal conversion constant from the metadata K2 from the metadata

idx = find(~cellfun('isempty',strfind(text,'K2_CONSTANT_BAND_10')));
K2 = string(text(idx));
K2 = extractAfter(K2,"= ");
K2 = str2double(K2);

% Solar elevation SE from the metadata

idx = find(~cellfun('isempty',strfind(text,'SUN_ELEVATION')));
SE = string(text(idx));
SE = extractAfter(SE,"= ");
SE = str2double(SE);

clearvars idx fid metadata_file

%% Calculate TOA spectral radiance 
% TOA = (ML * Qcal + AL) / sind(SE)
%
% where 
% ML = Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number).
% Qcal = corresponds to band 10
% AL = Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number).

TOA = ((ML * B10 + AL)./sind(SE));

%% TOA to Brightness Temperature conversion
% BT = (K2 / (ln (K1 / TOA) + 1)) - 273.15
%
% where 
% K1 = Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the thermal band number).
% K2 = Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the thermal band number).
% TOA (Top of Atmospheric) spectral radiance.

BT = (K2 ./ log((K1./TOA)+1)) - 273.15;

%% Calculate the NDVI
% NDVI = (Band 5 – Band 4) / (Band 5 + Band 4)

NDVI = (B5-B4)./(B5+B4);

%% Calculate the proportion of vegetation Pv
% Pv = ((NDVI – NDVImin) / (NDVImax – NDVImin))^2
% where 
% NDVI is the Normalized Difference Vegetation Index

PV = ((NDVI - min(NDVI(:))) ./ (max(NDVI(:)) - min(NDVI(:)))).^2;

%% Calculate Emissivity E
% E = 0.004 * Pv + 0.986
% where 
% Pv is the proportion of vegetation

E = (0.004 .* PV) + 0.986;

%% Calculate final Land Surface Temperature (LST)
% LST = (BT / (1 + (0.00115 * BT / 1.4388) * ln(E)))
% where
% BT = Brightness Temperature
% E = Emissivity

LST = (BT) ./ (1+(0.00115 .* (BT ./ 1.4388) .* log(E)));
LST(LST<-10)=NaN;

%% Plot result

LST(1:2449,:)=[];
LST(249:end,:)=[];
LST(:,1:6250)=[];
LST(:,275:end)=[];

imagesc(LST),colorbar
4

0 回答 0