使用 R 进行辐射校正
Posted
技术标签:
【中文标题】使用 R 进行辐射校正【英文标题】：Radiometric corrections with R 【发布时间】：2021-11-30 21:38:04 【问题描述】：我正在尝试通过 R / RStoolbox 将 Landsat 8 文件转换为反射率
我正在使用带有以下代码的脚本：
metaData <- readMeta("LC08_L1TP_183033_20210623_20210630_02_T1_MTL.txt")
lsat <- stackMeta("LC08_L1TP_183033_20210623_20210630_02_T1_MTL.txt")
lsat_sref <- radCor(lsat, metaData, method = "dos")
我遇到了这个错误：
Error in CRS(paste0(c("+proj=", "+zone=", "+units=m +datum="), pars, collapse = " ")) :
No spaces permitted in PROJ4 argument-value pairs: +proj= +zone= +units=m +datum=
你能帮帮我吗？
【问题讨论】：
你好，看那里：gis.stackexchange.com/questions/390084/… 非常感谢曼罗。我已经看过这个页面了。但这对我没有帮助。添加 raw=T 适用于 readMeta 函数，但不适用于 staskMeta 函数。所以我不知道如何通过。你有其他想法吗？ metaData 不，不。我从未使用过这个库。但我认为，如果你在那里问：gis.stackexchange.com - 你得到一些答案的机会已经长大了。 【参考方案1】：问题在于您的第二行代码。您需要按如下方式更改它，它应该可以工作：
lsat <- stackMeta(metaData)
请在下面找到一个小代表来向您展示它是如何工作的：
REPREX
library(RStoolbox)
mtlFile <- system.file("external/landsat/LT52240631988227CUB02_MTL.txt", package="RStoolbox")
metaData <- readMeta(mtlFile)
lsat <- stackMeta(metaData)
lsat_sref <- radCor(lsat, metaData, method = "dos")
lsat_sref
#> class : RasterStack
#> dimensions : 310, 287, 88970, 7 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 619395, 628005, -419505, -410205 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=22 +datum=WGS84 +units=m +no_defs
#> names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_bt, B7_sre
#> min values : 0.01327889, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 293.37508120, 0.00000000
#> max values : 0.2028536, 0.1935375, 0.2243498, 0.4356817, 0.3393489, 299.8284592, 0.2617160
由reprex package 创建于 2021-10-12 (v2.0.1)
编辑
summary()的输出
summary(metaData)
#> Scene: LT52240631988227CUB02
#> Satellite: LANDSAT5
#> Sensor: TM
#> Date: 1988-08-14
#> Path/Row: 224/63
#> Projection: +proj=utm +zone=22 +datum=WGS84 +units=m +no_defs
#>
#> Data:
#> FILES QUANTITY CATEGORY
#> B1_dn LT52240631988227CUB02_B1.TIF dn image
#> B2_dn LT52240631988227CUB02_B2.TIF dn image
#> B3_dn LT52240631988227CUB02_B3.TIF dn image
#> B4_dn LT52240631988227CUB02_B4.TIF dn image
#> B5_dn LT52240631988227CUB02_B5.TIF dn image
#> B6_dn LT52240631988227CUB02_B6.TIF dn image
#> B7_dn LT52240631988227CUB02_B7.TIF dn image
#>
#> Available calibration parameters (gain and offset):
#> dn -> radiance (toa)
#> dn -> brightness temperature (toa)
由reprex package (v2.0.1) 于 2021 年 10 月 12 日创建
编辑 2
我想我发现了你的问题。从网站https://earthexplorer.usgs.gov/ 下载时，除了 .mtl 元数据文件之外，您通常还会获得总共 12 个辐射波段（即 B1-B11，和一个“专用评估质量波段”-即 Bqa）。所有这些文件必须在同一个文件夹中，并且您不得删除其中任何一个（其中一些用于辐射校正）。换句话说，包RSToolbox 需要这些不同的文件。
如果你这样做，一切都会正常。
这是我下载的 landsat 8 图像的示例（不是代表），所有代码都有效：
library(RStoolbox)
metaData <- readMeta("DATA_09/Landsat8/Landsat8_RAW/LC08_L1TP_016030_20160927_20170220_01_T1_MTL.txt")
# NB : Please, note the presence of all the bands (B1-B11 + Band qa):
summary(metaData)
#> Scene: LC80160302016271LGN01
#> Satellite: LANDSAT8
#> Sensor: OLI_TIRS
#> Date: 2016-09-27
#> Path/Row: 16/30
#> Projection: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
#> Data:
#> FILES QUANTITY CATEGORY
#> B1_dn LC08_L1TP_016030_20160927_20170220_01_T1_B1.TIF dn image
#> B2_dn LC08_L1TP_016030_20160927_20170220_01_T1_B2.TIF dn image
#> B3_dn LC08_L1TP_016030_20160927_20170220_01_T1_B3.TIF dn image
#> B4_dn LC08_L1TP_016030_20160927_20170220_01_T1_B4.TIF dn image
#> B5_dn LC08_L1TP_016030_20160927_20170220_01_T1_B5.TIF dn image
#> B6_dn LC08_L1TP_016030_20160927_20170220_01_T1_B6.TIF dn image
#> B7_dn LC08_L1TP_016030_20160927_20170220_01_T1_B7.TIF dn image
#> B9_dn LC08_L1TP_016030_20160927_20170220_01_T1_B9.TIF dn image
#> B10_dn LC08_L1TP_016030_20160927_20170220_01_T1_B10.TIF dn image
#> B11_dn LC08_L1TP_016030_20160927_20170220_01_T1_B11.TIF dn image
#> B8_dn LC08_L1TP_016030_20160927_20170220_01_T1_B8.TIF dn pan
#> QA_dn LC08_L1TP_016030_20160927_20170220_01_T1_BQA.TIF dn qa
#>
#> Available calibration parameters (gain and offset):
#> dn -> radiance (toa)
#> dn -> reflectance (toa)
#> dn -> brightness temperature (toa)
lsat <- stackMeta(metaData)
lsat
#> class : RasterStack
#> dimensions : 504, 541, 272664, 10 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 318195, 334425, 4737735, 4752855 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
#> names : B1_dn, B2_dn, B3_dn, B4_dn, B5_dn, B6_dn, B7_dn, B9_dn, B10_dn, #> B11_dn
#> min values : 8492, 7617, 6777, 6026, 5290, 5185, 5142, 4995, 22328, #> 20666
#> max values : 19081, 19935, 21595, 22716, 27590, 38013, 39165, 5105, 30111, #> 27137
lsat_sref <- radCor(lsat, metaData, method = "dos")
lsat_sref
#> class : RasterStack
#> dimensions : 504, 541, 272664, 10 (nrow, ncol, ncell, nlayers)
#> resolution : 30, 30 (x, y)
#> extent : 318195, 334425, 4737735, 4752855 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
#> names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre, B9_sre, B10_bt, B11_bt
#> min values : 1.078127e-01, 7.822588e-02, 5.389152e-02, 3.082883e-02, 8.681202e-03, 5.291948e-03, 3.968224e-03, 0.000000e+00, 2.847589e+02, 2.830717e+02
#> max values : 4.347280e-01, 4.464235e-01, 5.032560e-01, 5.323340e-01, 6.762622e-01, 9.449113e-01, 9.547490e-01, 3.183901e-03, 3.039071e+02, 3.018956e+02
【讨论】：
嗨@Olivier Theureaux，如果我的回答解决了您的问题，请将其标记为已接受，以便其他用户更容易找到正确的答案。如果您不这样做，请写出问题所在。 感谢爱。不幸的是，它没有用。我有一个新错误 # Erreur dans file.exists(file) : argument 'file' wrong 好的。那是另一个问题。您收到的错误消息意味着文件“LC08_L1TP_183033_20210623_20210630_02_T1_MTL.txt”不可访问。有两种可能性：文件名在您的代码中包含错误，或者该文件不在您的工作目录中。在后一种情况下，您必须使用命令setwd() 设置工作目录（在括号中指明包含要读取的文件的文件夹的完整路径） C:/documents/mydirectory 放置斜杠（即“/”）而不是反斜杠（即“\”）。让我知道它是否有效。
元数据
这很好以上是关于使用 R 进行辐射校正的主要内容，如果未能解决你的问题，请参考以下文章
gdal学习笔记利用python 的gdal，以及相关库进行遥感图像处理（影像裁剪，辐射定标，大气校正，异常值去除）——以基于landsat8数据提取NDVI为例