使用R从Landsat 8计算地表温度（LST）

贝加尔湖

嘿！
I was recently asked by a user how to 计算ulate the 地表温度（LST） 来自Landsat 8影像，并决定撰写有关此主题的文章。我希望你会喜欢它。

让’从基础开始– 什么是LST？


陆地表面温度（LST）是地面的辐射表皮温度。这取决于反照率，植被覆盖度和土壤湿度。在大多数情况下，LST是植被和裸露的土壤温度的混合物。由于两者都会对由于云层覆盖和气溶胶负载变化以及光照的日变化而引起的太阳辐射的变化做出快速响应，因此LST也会显示快速变化。反过来，LST影响地面和植被之间的能量分配，并确定地面空气温度。

定义者 哥白尼全球土地服务

有多种计算LST的方法。我将向您展示官方USGS的方法 网页 使用  乐队10 和 11 来自 陆地卫星 8卫星的热红外传感器（TIRS）。 建议使用波段10进行定量分析，因为波段11的杂散光比波段10的污染多得多。

步骤1：原始图片

让’首先看一下从其中一个下载的原始图像 陆地卫星镜子。我的图片来自2015年7月25日，显示了俄罗斯西伯利亚雄伟的贝加尔湖及其周围的植被。（RGB 432）：

贝加尔湖

大自然对贝加尔湖非常慷慨。一条针叶林项链装饰着海岸山脉，创造出独特的景观美感，庇护动物和鸟类，并为周围环境增添活力。 （Bir.irk.ru)

贝加尔湖是世界上最大的淡水湖，约占世界的20％’未冻结的表面淡水。它的特点是环境非常异质，周围是广阔的森林地区，灌木丛和草原以及非常干燥的草原地区。人们可能会认为’西伯利亚的寒冷，但贝加尔湖的年平均气温实际上是11.8°C。 7月是最热的月份，平均气温为23°C。 （让’请牢记这些值以便进行LST分析。）

步骤2：将数字（DN）转换为大气（TOA）反射率的最高值

For our analysis 和 计算ulation of the LST we need 乐队10和乐队11 只要。下载的原始Landsat 8图像中的值称为传感器的数字编号。首先，必须使用以下公式将这些转换为TOA反射率值：

DN到TOA的反射率

from http://landsat.usgs.gov/Landsat8_Using_Product.php

要解决这个方程，您需要 RADIANCE_MULT_BAND_x 和 RADIANCE_ADD_BAND_x 从每个Landsat 8下载中包含的图元文件中获取。这些值（可能）因图像而异，因此请确保在每个场景中都查找它们！只需在文本编辑器中打开图元文件，并查找Band 10和Band 11的相应值：
元文件

现在我们可以如下转换R中的值：

#Values from 元文件
RADIANCE_MULT_BAND_10 <- 3.3420E-04
RADIANCE_MULT_BAND_11 <- 3.3420E-04

RADIANCE_ADD_BAND_10 <- 0.10000
RADIANCE_ADD_BAND_11 <- 0.10000

#Load 栅格 package 和 load band 10 & 11 into R (navigate to your image directory first)
library(raster)
band_10 <- 栅格("band10.tif") #change image name accordingly
band_11 <- 栅格("band11.tif") #change image name accordingly

＃从DN计算TOA：
toa_band10 <-calc（band_10，fun = 功能（x）{RADIANCE_MULT_BAND_10 * x + RADIANCE_ADD_BAND_10}）
toa_band11 <- 计算(band_11, fun=function(x){RADIANCE_MULT_BAND_11 * x + RADIANCE_ADD_BAND_11})

第三步：  转换为卫星亮度温度

In the next step we convert our reflectance values into 卫星Brightness温度，在这种情况下为Kelvin中的LST. This is done with the following equation 和 also needs some 数据 来自 metafile:

对LST的反射

http://landsat.usgs.gov/Landsat8_Using_Product.php

一旦在元文件中找到K1_CONSTANT_BAND_10，K2_CONSTANT_BAND_10，K1_CONSTANT_BAND_11和K2_CONSTANT_BAND_11的相应值，就可以使用R来计算LST，如下所示：

#Values from 元文件
K1_CONSTANT_BAND_10 <- 774.8853
K1_CONSTANT_BAND_11 <- 480.8883
K2_CONSTANT_BAND_10 <- 1321.0789
K2_CONSTANT_BAND_11 <- 1201.1442

#Calculate LST in Kelvin for 乐队10和乐队11
temp10_kelvin <- 计算(toa_band10, fun=function(x){K2_CONSTANT_BAND_10/log(K1_CONSTANT_BAND_10/x + 1)})
temp11_kelvin <- 计算(toa_band11, fun=function(x){K2_CONSTANT_BAND_11/log(K1_CONSTANT_BAND_11/x + 1)})

#Convert Kelvin to Celsius for Band 10 和 11
temp10_celsius <- 计算(temp10_kelvin, fun=function(x){x - 273.15})
temp11_celsius <- 计算(temp11_kelvin, fun=function(x){x - 273.15})

#Export 栅格 images
writeRaster(temp10_celsius, "temp10_c.tif")
writeRaster(temp11_celsius, "temp11_c.tif")

步骤4：在QGIS中显示最终结果

这就是结果的样子。下图描绘了自2015年7月25日起，西伯利亚贝加尔湖周围地区10带的LST值：

摄氏温度

您会看到贝加尔湖深处的云层（白色）和蓝色的冷蓝色水体，平均水面温度约为5-6°C。周围的灌木丛和森林的LST约为18-22°C。图片的南部干燥的草原和沙地的温度很高（因为图片是在炎热的夏季的中午左右拍摄的）。

我希望本教程对某些人有用。如果您对代码或转换有任何疑问，请随时在评论部分提出。有关更多详细信息，您可以访问以下网站：

 

干杯马丁

关于作者

马丁出生于捷克共和国，就读于维也纳自然资源与生命科学大学。他目前在GeoVille（一家位于奥地利的地球观测公司）工作，专门从事土地监测。他的主要兴趣是：开源应用程序，例如R，（地理空间）统计和数据管理，Web映射和可视化。他喜欢旅行，寻宝，摄影和运动。

33条留言

您可以在这篇文章中发表评论。

  • 你好
    我想使用R中的上述代码来计算Landsat8数据的LST，但是遇到以下错误。
    如果您能帮助我，我将不胜感激。
    ＃从DN计算TOA：
    toa_band10<-calc（band_10，fun = 功能（x）{RADIANCE_MULT_BAND_10 * x + RADIANCE_ADD_BAND_10}）
    错误：
    找不到函数的继承方法'calc' for signature '"character", "function"'

    马里赫 5年前 回复

    • 我将需要查看更多您的代码。似乎您的栅格不是’不能正确读入R，或者您的变量不是整数而是字符。您可以发布以下变量的输出吗？
      toa_band10
      RADIANCE_MULT_BAND_10
      RADIANCE_ADD_BAND_10

      马丁 5年前 回复

    • USgs已经在其分析就绪数据（ard）中发布了LST tiff，并且可以下载已处理的数据

      阿斯兰·瓦罗夸（Aslan Varoqua） 2年前 回复

  • 可以使用此代码将Landsat TM5的热带转换为LST产品，或者仅将其用于OLI数据。
    谢谢

    马里赫 5年前 回复

  • The 计算ulations done correctly. How do I get the image in QGIS ?
    使用什么颜色组合？

    谢谢

    罗纳尔多

    罗纳尔多 在4年前 回复

  • Tif fo的喀喀湖在秘鲁。
    以下链接生成了tif（temp10_c.tif）。

    //drive.google.com/file/d/0B8q_146s7i16VUo1clpfQTY1a2M/view?usp=sharing

    虚假色彩：

    //drive.google.com/file/d/0B8q_146s7i16bnJFNV9iaVVheUE/view?usp=sharing

    你能帮我弄个类似这个博客的tif吗

    谢谢

    罗纳尔多 在4年前 回复

  • 嗨，马丁。您能告诉我如何生成云掩码吗？非常感谢你。

    吉米 在4年前 回复

  • 在第3步中，您说，“卫星Brightness温度，在这种情况下为Kelvin中的LST”. Aren’这些不同的东西吗？看起来您正在计算传感器温度，但是我认为这与发射表面温度和大气效应不同，它与陆地表面温度（LST）不同。

    曼科夫 在4年前 回复

    • 确切地说，这里计算的值是传感器亮度温度下的值。为了计算陆地表面温度，将使用其他方法，例如分割窗口算法，温度发射率分离等。

      雷马 在4年前 回复

  • 亲爱的先生，
    感谢您对我们的帮助。
    您能帮我用R找出LAI叶面积指数吗？
    和 if you could help me to get the albedo 计算ulation 使用 R?

    问候

    艾哈迈德·沙（Ahmad Shah） 在4年前 回复

  • 感谢您向我们展示您的代码，它非常有帮助！

    赖德mond 在4年前 回复

  • 你好
    我从Landsat 8获得了城市LST，但是房屋的温度很低，但事实并非如此。
    你能帮我吗？

    新生 3年前 回复

    • 我不’我不知道你在哪里，但是我是第一次从Landsat 8波段10绘制热图，我的假设是草原会凉爽，市区会很热。所以当它反过来出现时，我感到非常惊讶。然而，经过深思熟虑，卫星在日出后约4小时的上午10点越过我的站点。在这一点上，只有任何面向太阳的建筑物的侧面都受到了直接的变暖。另一半在阴影下。因此，城市地区的任何像素都将平均达到多少阴影，以及被太阳直接加热了多少。另一方面，草原基本上提供了一个相当平坦但暴露的边界层，随着早晨的进行，该边界层越来越暴露于直射的阳光下。其他因素（例如土壤颜色）也可能会影响这些草地的价值。不幸的是，我当时’t可以检验这些解释，但它们对输出有意义，并且是从不信任假设的另一个很好的指示，直到他们’重新测试。我怀疑Landsat是在午后还是下午结束时，市区的温度会比这些草原高得多，但不幸的是Landsat的数据没有’让我们检验这个假设。

      戴夫·弗尼斯（Dave Furniss） 3年前 回复

      • 非常有趣的见解！

        马丁 3年前 回复

  • 嗨，马丁，
    我想引用您的解释，并希望您的名字被MLA引用。
    很高兴与您合作。
    谢谢

    杰克 3年前 回复

    • 嘿，杰克，
      多么荣幸。您可以使用以下参考：
      Siklar，M.（2016年）。使用R.从Landsat 8计算陆地表面温度（LST）。[Blog] 地理信息系统博客.com。可在： http://www.citystoragesrq.com/wp-admin/post.php?post=872 [2017年11月17日访问]。

      让我知道您的作品何时出版。我确定’s a good read.

      干杯马丁

      马丁 3年前 回复

  • 这是在亮度温度下…。要计算陆地表面温度，需要添加发射率..对吧？

    尤努斯 3年前 回复

  • 似乎是您计算了亮度温度，土地表面温度将使用发射率和其他常数进行计算：
    LST =亮度温度/ 1 + Wevelenghth（10微米）*（亮度温度/ P）* Ln（e）
    P = h * c / s = 14380
    h =木板常数
    c =光速
    s =玻尔兹曼常数
    e =发射率
    e = 0.004 * Pv + 0.986
    Pv =植被比例=（NDVI-NDVImin / NDVImax-NDVImin）^ 2
    LST 计算ulation by NDVI Threshold Method (Sobrino et al, 2014)

    拉米因 3年前 回复

    • 波长数取决于您使用的波段。频段10为10.6，频段11为11.5

      拉明 3年前 回复

  • 你好马丁
    我们正在研发LANDSAT 7&8，我们要下载10年的数据。掩盖的区域是LANDSAT的16瓦。我需要一个代码来帮助我以较高的下载速度下载批量图像。非常感谢您的帮助

    莫欣·瓦卡斯（Mohsin Waqas） 3年前 回复

  • 非常感谢马丁

    如果您想知道后台发生了什么，上面的代码绝对是很棒的。也肯定会拓宽一个’s coding skills.

    但是，如果要在多个LANDSAT文件上执行上述任务，并希望以最少的错别字更快地完成该任务，则可能需要查看以下选项。
    http://terradue.github.io/rLandsat8/section3/bt.html
    这种方法摆脱了背景像素。但是，结果以开尔文为单位’我想添加另一行代码，将数据转换为celcius并导出.tif文件

    彼得·卡巴诺 3年前 回复

    • 嘿，彼得！
      感谢您分享链接。我没’不知道有一个rLandsat软件包。我相信这对许多读者来说非常有用。
      最好的祝福

      马丁

      马丁 3年前 回复

  • 您好马丁，非常感谢您提供的信息，它’s great.
    您能说明如何使用LANDSAT 5吗？
    最好的祝福！

    西蒙 2年前 回复

  • 嗨，马丁，
    我希望你一切都好。
    您能否告诉我如何计算城乡温差。如果有的话，也请分享任何教程。

    谢谢

    最好的祝福！
    纳亚卜

    纳亚卜 1年前 回复

  • 非常感谢！

    天尾 1年前 回复

  • 非常感谢马丁，
    但是有什么可以帮助我使用R脚本使用LST和NDVI / EVI计算温度植被干燥指数（TVDI）或植被温度条件指数（VTCI）

    谢谢

    阿卜杜拉拉西德·哈桑 10个月前 回复

  • 您能否告诉我，LST是否取决于月份？

    阿努帕姆·达斯（Anupam Das） 7个月前 回复

  • 嗨，我计算Landsat 8的LST，但最后使用导出查询时，出现警告消息
    警告信息：
    在.gd_SetProject（object，…）：未针对PROJ更新>= 6
    您能帮我解决这个问题吗？

    穆斯塔法526 7个月前 回复

发表回复

验证码

*