多景Landsat影像平均地表温度(LST)计算异常求助

使用2023年5月26日、9月7日、9月15日同一分幅（PATH:201，ROW:24）的Landsat 8（Level 2）数据，计算该区域跨三个日期的平均地表温度（LST，单位为摄氏度）。代码可正常运行，但最终生成的cels栅格数值范围为-11至0℃，与所选日期的实际温度逻辑不符。参照《how to calculate land surface temperature with landsat 8 images》教程实现计算，但无法定位错误，恳请专业人士提供帮助！

代码实现 #

library(colorBlindness)
library(terra)
library(sf)
library(ggplot2)
library(ggspatial)
library(readxl)
library(spdep)
library(tmap)
library(tidyverse)
library(grDevices)

# 读取IMD数据
imd <- read_sf("data/english_IMD_2019/English IMD 2019/IMD_2019.shp")

# 读取伦敦边界数据
london <- read_sf("data/Output_Area_Classification_2011/Local_Enterprise_Partnerships/E37000023/shapefiles/E37000023.shp")
london_imd <- st_intersection(imd, london)
london_poly <- vect(london)

# 读取LST相关波段（ST_B10）
may_26 <- rast("data/LC08_L2SP_201024_20230526_20230603_02_T1/LC08_L2SP_201024_20230526_20230603_02_T1_ST_B10.TIF")
sept_15 <- rast("data/LC08_L2SP_201024_20230915_20230925_02_T1/LC08_L2SP_201024_20230915_20230925_02_T1_ST_B10.TIF")
sept_7 <- rast("data/LC09_L2SP_201024_20230907_20230913_02_T1/LC09_L2SP_201024_20230907_20230913_02_T1_ST_B10.TIF")

# 读取SR_B4波段
may_26_b4 <- rast("data/LC08_L2SP_201024_20230526_20230603_02_T1/LC08_L2SP_201024_20230526_20230603_02_T1_SR_B4.TIF")
sept_15_b4 <- rast("data/LC08_L2SP_201024_20230915_20230925_02_T1/LC08_L2SP_201024_20230915_20230925_02_T1_SR_B4.TIF")
sept_7_b4 <- rast("data/LC09_L2SP_201024_20230907_20230913_02_T1/LC09_L2SP_201024_20230907_20230913_02_T1_SR_B4.TIF")

# 读取SR_B5波段
may_26_b5 <- rast("data/LC08_L2SP_201024_20230526_20230603_02_T1/LC08_L2SP_201024_20230526_20230603_02_T1_SR_B5.TIF")
sept_15_b5 <- rast("data/LC08_L2SP_201024_20230915_20230925_02_T1/LC08_L2SP_201024_20230915_20230925_02_T1_SR_B5.TIF")
sept_7_b5 <- rast("data/LC09_L2SP_201024_20230907_20230913_02_T1/LC09_L2SP_201024_20230907_20230913_02_T1_SR_B5.TIF")

# 裁剪范围
new_ext <- ext(203985, 444615, 5606985, 5851815)

may_26 <- crop(may_26, new_ext)
sept_15 <- crop(sept_15, new_ext)
sept_7 <- crop(sept_7, new_ext)

may_26_b4 <- crop(may_26_b4, new_ext)
sept_15_b4 <- crop(sept_15_b4, new_ext)
sept_7_b4 <- crop(sept_7_b4, new_ext)

may_26_b5 <- crop(may_26_b5, new_ext)
sept_15_b5 <- crop(sept_15_b5, new_ext)
sept_7_b5 <- crop(sept_7_b5, new_ext)

# 投影转换
london_poly <- project(london_poly, crs(may_26))

# 掩膜处理
may_26 <- mask(may_26, london_poly)
sept_15 <- mask(sept_15, london_poly)
sept_7 <- mask(sept_7, london_poly)

may_26_b4 <- mask(may_26_b4, london_poly)
sept_15_b4 <- mask(sept_15_b4, london_poly)
sept_7_b4 <- mask(sept_7_b4, london_poly)

may_26_b5 <- mask(may_26_b5, london_poly)
sept_15_b5 <- mask(sept_15_b5, london_poly)
sept_7_b5 <- mask(sept_7_b5, london_poly)

# 计算多日期平均值
coll <- c(may_26, sept_15, sept_7)
mos <- app(coll, fun = mean, na.rm = TRUE)

coll_b4 <- c(may_26_b4, sept_15_b4, sept_7_b4)
mos_b4 <- app(coll_b4, fun = mean, na.rm = TRUE)

coll_b5 <- c(may_26_b5, sept_15_b5, sept_7_b5)
mos_b5 <- app(coll_b5, fun = mean, na.rm = TRUE)

# 辐射定标参数
RADIANCE_MULT_BAND_10 <- 3.3420E-04
RADIANCE_ADD_BAND_10 <- 0.10000

TOA <- RADIANCE_MULT_BAND_10 * mos + RADIANCE_ADD_BAND_10

# 亮度温度计算参数
K1_CONSTANT_BAND_10 <- 774.8853
K2_CONSTANT_BAND_10 <- 1321.0789

bt <- (K2_CONSTANT_BAND_10 / (log((K1_CONSTANT_BAND_10 / TOA)) + 1)) - 273.15

# NDVI计算
ndvi <- (mos_b5 - mos_b4) / (mos_b5 + mos_b4)

min_ndvi <- global(ndvi, fun = "min", na.rm = TRUE)[[1]]
max_ndvi <- global(ndvi, fun = "max", na.rm = TRUE)[[1]]

# 植被覆盖度计算
pv <- (ndvi - min_ndvi) / (max_ndvi - min_ndvi)

# 发射率计算
e <- 0.004 * pv + 0.986

# LST计算
cels <- (bt / (1 + (0.00115 * bt / 1.4388) * log(e)))

# 绘图
plot(
    cels,
    col = Blue2Orange10Steps,
    xlim = c(255000, 317000),
    ylim = c(5685000, 5734000),
)

如前所述，参照教程实现代码，预期能正确计算所选三个日期的平均LST（摄氏度），但结果出现异常。

内容的提问来源于stack exchange，提问作者trizzo