How can I find water surface temperature from landsat-9 bands?

I would like to find water surface temperature by using Lansat-9 bands.

Thanks in advance,

Riyadh

[Riyadh]

% this one hasn't work!
B10 = double(imread('LC08_L2SP_169037_20200821_20200905_02_T1_ST_B10.tif'));
% Retrieve metadata parameters
% (Extract ML, AL, K1, and K2 from the metadata file)
% Calculate Top of Atmospheric (TOA) spectral radiance
ML=3.3420E-04;
AL=0.1;
TOA = ML .* B10 + AL;
% Convert radiance to brightness temperature (in Kelvin)
K2=1321.0789;
K1=774.8853;
T_B = K2 ./ (log(K1 ./ (TOA + 1)));
% Adjust for water surface emissivity (assume ~0.98)
emissivity = 0.98;
% Calculate Land Surface Temperature (LST)
rho = 1.4388e-2; % Stefan-Boltzmann constant
LST = T_B ./ (1 + (T_B ./ rho) .* log(emissivity));

[Umar]

Hi Riyadh,

Your aim is to calculate the Land Surface Temperature (LST) from a thermal band image. The code reads a thermal band image, extracts metadata parameters, calculates the Top of Atmospheric (TOA) spectral radiance, converts radiance to brightness temperature, and finally computes the Land Surface Temperature (LST). Sounds reasonable.

The reason your above code is not working is due to incorrect adjustment for water surface emissivity in the calculation of Land Surface Temperature (LST). The formula used to adjust for emissivity is not applied properly, leading to inaccurate temperature calculations.

Hope this will help resolve your problem.

Answer 1

[Umar]

To calculate water surface temperature using Landsat-9 bands in Matlab, you can utilize the following functions:

multibandread: This function can be used to read Landsat-9 bands data.

data = multibandread('your_image.tif', [rows, cols, bands], 'datatype', 'single', 'ieee-le', 'bsq');

band2tir: This function converts Landsat-9 band data to Top of Atmosphere (TOA) brightness temperature.

tir = band2tir(data, 'Landsat8', '10', 'B10');

radiance2bt: This function converts radiance to brightness temperature.

bt = radiance2bt(tir, 'Landsat8', 'B10');

To analyze water surface temperature using Landsat-9 bands in Matlab, you can follow these steps:

Import Landsat-9 Data: Load the Landsat-9 bands into Matlab using the multibandread function.

For more information on multibandread function, please refer to

https://www.mathworks.com/help/matlab/ref/multibandread.html

Convert Digital Numbers to Radiance: Convert the digital numbers to radiance values using the provided conversion coefficients.

Calculate Brightness Temperature: Utilize the Planck's law equation to calculate brightness temperature from radiance.

Apply Atmospheric Correction: Implement atmospheric correction algorithms to enhance the accuracy of temperature estimation.

Derive Water Surface Temperature: Use specific band combinations and algorithms to derive water surface temperature from Landsat-9 data.

By combining these functions and processing the Landsat-9 bands data accordingly, you can derive the water surface temperature effectively.

[Walter Roberson]

Hmmm, I would not expect multibandread to work on .tif files. .tif files contain tagged chunks of data, but the order and size of the chunks is not fixed; multibandread() expects files to be binary blobs of data.

[Umar]

Hi Walter,

Perhaps you can use function like imread for handling specific file formats like .tif. For more information, please refer to https://www.mathworks.com/help/matlab/ref/imread.html Please correct me if I am wrong.

[Walter Roberson]

imread() could potentially work -- provided that the data of interest happened to be recorded in the first 3 bit planes. Otherwise you would need to use geotiffread or readgeoraster (mapping toolbox) or use the Tiff interface.

[Image Analyst]

Actually I think @Riyadh has already figured out how to extract the bands. Look at this filename : 'LC08_L2SP_169037_20200821_20200905_02_T1_ST_B10.tif' That looks Band 10 that has been extracted into a single gray scale TIFF image. So, once that's been done, imread() should work fine.

[Umar]

I agree

Answer 2

[Image Analyst]

Use multibandread to read bands 10 and 11. I'm not sure of the difference between those two bands (presumably you do), but their website says those are the two thermal bands:

https://landsat.gsfc.nasa.gov/satellites/landsat-9/landsat-9-bands/

help multibandread

 MULTIBANDREAD Read band interleaved data from a binary file
    X = MULTIBANDREAD(FILENAME,SIZE,PRECISION,
                       OFFSET,INTERLEAVE,BYTEORDER)
    reads band-sequential (BSQ), band-interleaved-by-line (BIL), or
    band-interleaved-by-pixel (BIP) data from a binary file, FILENAME.  X is
    a 2-D array if only one band is read, otherwise it is 3-D. X is returned
    as an array of data type double by default.  Use the PRECISION argument
    to map the data to a different data type.
 
    X = MULTIBANDREAD(FILENAME,SIZE,PRECISION,OFFSET,INTERLEAVE,
                     BYTEORDER,SUBSET,SUBSET,SUBSET)
    reads a subset of the data in the file. Up to 3 SUBSET parameters may be
    used to subset independently along the Row, Column, and Band dimensions.
 
    In addition to BSQ, BIL, and BIP files, multiband imagery may be stored 
    using the TIFF file format.  In that case, the data should be imported
    with IMREAD.
 
    Parameters:
 
      FILENAME: A character vector or string scalar containing the name of
      the file to be read.
 
      DIMS: A 3 element vector of integers consisting of
      [HEIGHT, WIDTH, N]. HEIGHT is the total number of rows, WIDTH is
      the total number of elements in each row, and N is the total number
      of bands. This will be the dimensions of the data if it read in its
      entirety.
 
      PRECISION: A character vector or string scalar to specify the format
      of the data to be read. For example, 'uint8', 'double', 'integer*4'.
      By default X is returned as an array of class double. Use the
      PRECISION parameter to format the data to a different class.  For
      example, a precision of 'uint8=>uint8' (or '*uint8') will return the
      data as a UINT8 array. 'uint8=>single' will read each 8 bit pixel and
      store it in MATLAB in single precision.  MULTIBANDREAD will attempt
      to use the efficient MEMMAPFILE function if the precision value
      corresponds to a native MATLAB type.  See the help for FREAD for a
      more complete description of PRECISION.
 
      OFFSET: The zero-based location of the first data element in the file.
      This value represents number of bytes from the beginning of the file
      to where the data begins.
 
      INTERLEAVE: The format in which the data is stored.  This can be
      either 'bsq','bil', or 'bip' for Band-Sequential,
      Band-Interleaved-by-Line or Band-Interleaved-by-Pixel respectively.
 
      BYTEORDER: The byte ordering (machine format) in which the data is
      stored. This can be 'ieee-le' for little-endian or 'ieee-be' for
      big-endian.  All other machine formats described in the help for FOPEN
      are also valid values for BYTEORDER.
 
      SUBSET: (optional) A cell array containing either {DIM,INDEX} or
      {DIM,METHOD,INDEX}. DIM is one of three values: 'Column', 'Row', or
      'Band' specifying which dimension to subset along.  METHOD is 'Direct'
      or 'Range'. If METHOD is omitted, then the default is 'Direct'. If
      using 'Direct' subsetting, INDEX is a vector specifying the indices to
      read along the Band dimension.  If METHOD is 'Range', INDEX is a 2 or
      3 element vector of [START, INCREMENT, STOP] specifying the range and
      step size to read along the dimension. If INDEX is 2 elements, then
      INCREMENT is assumed to be one.
 
    Examples:
 
    % Setup initial parameters for a dataset.
    rows=3; cols=3; bands=5;
    filename = tempname;
 
    % Define the dataset.
    fid = fopen(filename, 'w', 'ieee-le');
    fwrite(fid, 1:rows*cols*bands, 'double');
    fclose(fid);
 
    % Read the every other band of the data using the Band-Sequential format.
    im1 = multibandread(filename, [rows cols bands], ...
              'double', 0, 'bsq', 'ieee-le', ...
              {'Band', 'Range', [1 2 bands]} )
 
    % Read the first two rows and columns of data using
    % Band-Interleaved-by-Pixel format.
    im2 = multibandread(filename, [rows cols bands], ...
              'double', 0, 'bip', 'ieee-le', ...
              {'Row', 'Range', [1 2]}, ...
              {'Column', 'Range', [1 2]} )
 
    % Read the data using Band-Interleaved-by-Line format.
    im3 = multibandread(filename, [rows cols bands], ...
              'double', 0, 'bil', 'ieee-le')
 
    % Delete the file that we created.
         delete(filename);
 
    % The FITS file 'tst0012.fits' contains int16 BIL data starting at
    % byte 74880.
    im4 = multibandread( 'tst0012.fits', [31 73 5], ...
              'int16', 74880, 'bil', 'ieee-be', ...
              {'Band', 'Range', [1 3]} );
    im5 = double(im4)/max(max(max(im4)));
    imagesc(im5);
 
    See also FREAD, FWRITE, IMREAD, MEMMAPFILE, MULTIBANDWRITE.

    Documentation for multibandread
       doc multibandread