我正在尝试解决一个我认为相当简单的问题，但产生了错误的结果。我从同一个图块（路径：201，行：24）的 3 个不同日期（5 月 26 日、9 月 15 日、9 月 7 日，全部在 2023 年）获取了 Landsat 8（2 级）数据。我想获得三个日期的平均地表温度值（以摄氏度为单位）。为此，我执行了下面的代码，该代码运行良好，但最终出现在
celslibrary(colorBlindness)
library(terra)
library(sf)
library(ggplot2)
library(ggspatial)
library(readxl)
library(spdep)
library(tmap)
library(tidyverse)
library(grDevices)
# https://data.cdrc.ac.uk/dataset/index-multiple-deprivation-imd
imd <- read_sf("data/english_IMD_2019/English IMD 2019/IMD_2019.shp")
# https://data.cdrc.ac.uk/system/files/Output_Area_Classification_2011/Output_Area_Classification_2011_E37000023.zip
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)
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")
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")
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 <- (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
cels <- (bt / (1 + (0.00115 * bt / 1.4388) * log(e)))
plot(
cels,
col = Blue2Orange10Steps,
xlim = c(255000, 317000),
ylim = c(5685000, 5734000),
)
如上所述，我尝试遵循这篇文章，并期望上面的代码能够正确计算所选三个日期的平均 LST（以摄氏度为单位）。
TOA辐射率和LST之间的关系并不是简单的线性关系。因此，您不应该从 3 个图像中获取平均辐射率，然后应用一次公式。相反，您应该分别计算每个图像的 LST，并仅在最后一步取平均值。
类似这样的：
coll <- c(may_26, sept_15, sept_7)
<<< calculate TOA, BT, NDVI, PV, e and cels>>>
mos <- app(cels, fun = mean, na.rm = TRUE)