How to calculate land surface temperture from Landsat8 imagery in Matlab?

조회 수: 5(최근 30일)
Yoni Verhaegen -WE-1718-
Yoni Verhaegen -WE-1718- 2021년 8월 21일
I have downloaded some Landsat 8 imagery and I want to calculate the LST (Land Surface Temperature) from these images in MATLAB.
However, something does not seem to go right. For example, when I checked the calculated temperatures over a glacier surface, the displayed values are more than 44 degrees Celcius. I don't think that this is possible, as a melting ice surface is of course not able to get positive higher than 0 degrees. Can anyone help me to check where things go wrong? I have not been able to find the error yet.
The code I have used can be found under the image. The formulas I used were found on: https://giscrack.com/how-to-calculate-land-surface-temperature-with-landsat-8-images/. I have uploaded the corresponding geotiff images at: https://drive.google.com/file/d/1CXxt7hHWNbXujt6rMQPgGw103Zo6vbIc/view?usp=sharing.
Thanks!
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 K1 from the metadata
idx = find(~cellfun('isempty',strfind(text,'K2_CONSTANT_BAND_10')));
K2 = string(text(idx));
K2 = extractAfter(K2,"= ");
K2 = str2double(K2);
clearvars idx fid metadata_file
%% Calculate TOA spectral radiance
% TOA = ML * Qcal + AL
%
% 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;
%% 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 = real(1321.0789 ./ log(774.8853./TOA +1.0)) - 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(isnan(LST))=0;
%% Plot result
figure,
imagesc(LST),colorbar
caxis([0 80])

답변(0개)

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by