1 Correção radiométrica

Normalmente, cada detector no plano focal de um sistema de sensoriamento remoto teria a mesma resposta (média e padrão), quando estimulados pelo mesmo sinal luminoso de entrada. Assim, uma imagem de um brilho uniforme deve estar livre de qualquer variação. No entanto, isso raramente é o caso, as imagens brutas produzidas contêm variações visivelmente discerníveis e mensuráveis, chamadas ruído de padrão fixo. Essas intensidades espaciais podem ser atribuídas a diferenças no detector individual de foto-resposta, processamento de sinal eletrônico, variação em função do campo de visão do sistema óptico e da presença de contaminantes dentro da carga útil meio Ambiente (ANDERSON et al., 2011), daí a necessidade de correção radiométrica dos produtos (imagens de uma cena) de um sensor.

De modo que a correção radiométrica é o procedimento de conversão da intensidade detectada pelo sensor em unidades físicas de energia significativas para análise no campo de sensoriamento remoto.

1.1 Sensores e produto do satélite Landsat 8

O satélite Landsat 8 possui dois sensores montados na plataforma de observação da Terra: o OLI (Operational Land Imager) e o TIRS (Termal Infrared Sensor). O OLI realiza imageamento nas regiões do espectro eletromagnético: visível, infravermelho próximo (NIR) e infravermelho de ondas curtas (SWIR), dotado de 9 bandas multiespectrais de 30 metros de resolução espacial cada, com exceção da banda pancromática (banda 8), que possui resolução espacial de 15 metros. Algumas características das bandas do sensor OLI podem ser visualizadas na tabela e gráfico logo abaixo:

Tabela 1 - Bandas multiespectrais do sensor OLI
Banda Nome Comprimento.de.onda.central_nm Largura.de.banda_nm Resolucao.espacial_m
Band 1 Coastal/Aerossol 443.0 16 30
Band 2 Blue 482.0 60 30
Band 3 Green 561.5 57 30
Band 4 Red 654.5 37 30
Band 5 NIR 865.0 28 30
Band 6 SWIR 1 1608.5 85 30
Band 7 SWIR 2 2200.5 187 30
Band 8 Pan 589.5 173 15
Band 9 Cirrus 1373.5 21 30

Fonte: adaptado de U.S. Geological Survey (2015)

Já o sensor TIRS, realiza imageamento em duas bandas multiespectrais, na região do espectro infravermelho termal. Ambas possuem resolução espacial de 100 metros, com especificações (comprimento de onda central, largura de banda, e nome da banda) dadas pela tabela abaixo:

Tabela 2 - Bandas multiespectrais do sensor TIRS
Banda.T Nome.T Comprimento.de.onda.central_nm.T Largura.de.banda_nm.T Resolucao.espacial_m.T
Band 10 Thermal Infrared: TIR-1 10895 590 100
Band 11 Thermal Infrared: TIR-2 12005 1010 100

Fonte: adaptado de U.S. Geological Survey (2015)

Os produtos Landsat 8 são distribuídos oficialmente pela USGS (United States Geological Survey) através do portal Earth Explorer (pelo link https://earthexplorer.usgs.gov). Os produtos Landsat 8 são distribuídos como produto de nível 1. Esses produtos são compostos de imagens Landsat 8 nas 11 bandas detalhadas acima, ortoretificadas, georregistradas, com resolução radiométrica de 16 bits, ou seja, com 65.536 níveis de cinza e portanto, com pixels escalados no formato de número digital (DN). Além dessas 11 imagens de uma cena (Órbita/Ponto), acompanham o produto dois arquivos de texto, contendo os metadados (documentação) das imagens: são os arquivos ANG e MTL. A figura abaixo mostra um produto de nível 1 dos sensores OLI e TIRS, contendo as 11 bandas de 16 bits, além dos arquivos ANG e MTL:

knitr::include_graphics("Figura.PNG")

1.2 Objetivo da rotina

Conhecida a necessidade de conversões radiométricas dos produtos gerados pelo sensor orbital, bem como tendo conhecimento desses produtos, a rotina apresentada neste tutorial tem como objetivo realizar a conversão radiométrica de imagens Landsat 8 aplicando funções da biblioteca landsat8 para obter radiância (em \(W/(m^2 sr^{-1} \mu m^{-1})\)) das bandas 1 a 11, reflectância (\(admensional\)) no topo da atmosfera (TOA) das bandas 1 a 9, e temperatura de brilho (em \(K\)) das bandas 10 e 11.

2 Bibliotecas usadas na rotina

O objetivo desta rotina é realizar a correção radiométrica de imagens Landsat 8 usando funções da biblioteca landsat8. Contudo, funções de outros pacotes foram necessárias, como: dplyr, para tornar o código mais legível; a função glue da biblioteca glue, para criar nomes de variáveis unindo variáveis e caracteres; para importar dados raster como FormalClass do tipo “SpatialGridDataFrame” (exigido pelas funções de landsat8) foi necessária a função readGDAL da biblioteca rgdal; e por fim, para plotar resultados e exportar imagens raster, foram usadas funções das bibliotecas sp e raster, respectivamente. De modo que todas as bibliotecas são listadas no chunk abaixo.

#install.packages("landsat8")
#install.packages("sp")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("dplyr")
#install.packages("glue")
library(landsat8)
library(sp)
library(raster)
library(rgdal)
library(dplyr)
library(glue)

Caso o usuário não possua instalado os pacotes listados necessários, os comandos comentados #install.packages(“landsat8”) podem ser usados para instalar (apenas descomente-os!).

3 Rotina para correção radiométrica

3.1 Criando variáveis

A primeira etapa, e é onde o usuário interage com o código fonte, é a criação de três variáveis: caminho do produto, sufixo do produto, e número específico da banda. O caminho deve levar até a pasta onde estão as imagens e arquivos de texto ANG e MTL, pois, por mais que a pasta de produto nível 1 do Landsat 8 esteja dentro do projeto R Studio do usuário, ainda assim é necessário indicá-la. Feito isso, usa-se a função setwd para configurar a variável caminho como caminho do produto. A imagem abaixo ilustra como encontrar o conteúdo da variável caminho:

knitr::include_graphics("caminho.png")

Copie o texto selecionado em azul na imagem, e cole na criação da variável caminho (note que as barras “" devem ser invertidas para”/“), em seguida configure a variável como diretório (usando a função setwd) onde o R irá buscar os dado, como demostrado no código abaixo:

#caminho <- "C:/Users/ERLI/Downloads/LC08_L1TP_215073_20190602_20190602_01_T1/"
#setwd(caminho)

Porquê o chunk acima está comentado?

Essa primeira etapa é válida apenas quando se está trabalhando com um script .R. O knitr não permite uso do comando setwd() para definir diretório (para mais detalhes consulte a seção Note pelo comando ?knitr::knit), portanto neste tutorial os dados do produto Landsat 8 foram importados para a pasta onde se encontra o arquivo .Rmd, como ilustra a figura abaixo:

knitr::include_graphics("produtos.png")

A segunda variável a ser criada é o sufixo do produto, isto é, um trecho em comum dos nomes de arquivos no produto Landsat 8 Level 1. Veja na figura abaixo, do arquivo usado para exemplificar esta rotina, como cada arquivo possui um trecho do nome em comum:

knitr::include_graphics("produto.png")

Selecione esse nome, copie e cole na variável produto, como demonstrado no código abaixo. E por último, é necessário uma variável que faça variar a banda nas etapas posteriores (loops). Para isso é criada uma variável band, para variar o número da banda, como sequência de 1 a 11:

produto <- "LC08_L1TP_215073_20190602_20190602_01_T1_"
band <- c(1:11)

Note que nesse tópico, denominado “criando variáveis”, não são criadas variáveis para as imagens. Por quê? A rotina desenvolvida foi pensada para ser processada em hardwares com RAM em torno de 4 GB, de maneira que, trazer para o Global Environment as imagens que variam de tamanho de 150 a 480 MB, é inviável e impossibilita o processamento de todas as imagens. Portanto a rotina foi programada para pegar a imagem no disco rígido (no caminho configurado), processar e escrevê-la no disco (no mesmo caminho configurado).

3.1.1 Importanto metadados do produto

O procedimento aqui aplicado para correção atmosférica é o descrito por U.S. Geological Survey (2015) no Landsat 8 User Handbook. O procedimento realiza correções baseado em fatores de conversão de imagem, e estes são exigidos pelas funções do pacote landsat8 para realizar os cálculos.

São eles: Ml, Al, Mp, Ap, K1, K2 e sunelev. E o que são? Ml trata-se do RADIANCE_MULT_BAND_b, onde b é o número da banda: este é um fator de reescalonamento multiplicativo (offset) para radiância. De modo análogo, Mp é p fator offset para reflectância: encontrado nos metadados como REFLECTANCE_MULT_BAND_b.

Al e Ap são fatores aditivos de reescalonamento (gain), encontrados no metadados como RADIANCE_ADD_BAND_b e REFLECTANCE_ADD_BAND_b, respectivamente, onde b é o número da banda. K1 e K2 são constantes de conversão térmica específica da banda (para as bandas 10 e 11), dados no arquivo MTL como K1_CONSTANT_BAND_b e K2_CONSTANT_BAND_b, respectivamente, onde b é o número da banda.

A variável sunelev trata-se do ângulo de elevação solar no momento de aquisição das imagens de uma cena (órbita/ponto), é dado em grau decimal e é obtido a partir também dos metadados como: SUN_ELEVATION. É usado para calcular a reflectância TOA (Top Of Atmosphere).

Como esses arquivos são obtidos nos metadados (arquivo MTL), a importação deste é realizada usando o comando read.table (do pacote utils, nativo do R). É necessário fixar o número de linhas em 224, pelo seguinte motivo: a importação deve ser feita como data.frame (mesmo número de linhas em todas as colunas), e por a linha 225 possuir apenas uma variável (e não três, como as demais), a importação não seria realizada com essa estrutura.

MTL <- read.table(glue::glue("{produto}MTL.txt"), nrows = 224)
rmarkdown::paged_table(MTL, options = list(rows.print = 10, cols.print = 3))

Temos os metadados, podemos processar imagens?

Ainda não! Em operações de conversão usaremos lógica nos loops (afinal, são muitas imagens, e não vamos fazer uma a uma!), portanto vamos organizar os dados de forma lógica para essas etapas.

Feita a importação dos metadados, precisamos das variáveis descritas acima para realizar as conversões de imagens. Essas variáveis, precisam estar em data frames, pois, são específicas de cada imagem, com exceção do ângulo de elevação solar.

A estrutura de obtenção desses dados é escrita com funções do pacote dplyr, desenvolvido para ser uma gramática de manipulação de dados. Usaremos o operador ‘piper’ (%>%) para colocar funções e dados como argumentos. Então, que funções iremos usar? As seguintes:

select(): função da biblioteca dplyr para organização de colunas, essa função seleciona colunas pelo nome ou funções auxiliares;

slice(): função da biblioteca dplyr para selecionar linhas pela posição desta;

%>%: operador da biblioteca dplyr que coloca o objeto do lado esquerdo como argumento de uma função do lado direito.

Observe no chunk abaixo a obtenção de um data frame para o ângulo de elevação solar:

sun.elevation <- MTL %>% select(3) %>% slice(77) %>% type.convert()

Crie a variável sun.elevation, em seguida declare que essa variável será preenchida com dados do data frame MTL.

No data frame MTL, procure na coluna V1 a variável SUN_ELEVATION, e identifique o número (posição) da sua linha: qual a posição (número) dessa linha? Observe também na coluna V3, a classe de dados: veja que a importação trouxe essa coluna como tipo factor, certo? Guarde ambas as informações.

Note que depois da igualdade (<-), há a seguinte estrutura: MTL %>% select(3) %>% slice(77) %>% type.convert(); o que significa? Significa que, no data frame MTL, estamos ordenando que R selecione a coluna 3 (a coluna que contém os valores que nos interessa), e que nessa coluna 3 ele obtenha o item da linha 77 (a linha que você identificou a variável SUN_ELEVATION), em seguida que converta a classe desse item para numeric (precisaremos dela como número, e você observou que a classe no data frame MTL é como fator).

Como ficaria sem o %>% ?

Sem o operador piper o código ficaria um tanto menos didático (veja este comentado abaixo) e mais difícil de compreender (pense em funções maiores):

#sun.elevation <- type.convert(slice(select(MTL, 3), 77))

Obtenha as demais variáveis!

De maneira análoga à criação do data frame sun.elevation, crie 2 data frames: um para valores de radiância (gain e offset), e outro para reflectância TOA (gain e offset, também). Dessa vez como estamos juntando mais de uma coluna, precisamos englobar ambas as seleções em um data.frame, portanto use a função data.frame().

Você já sabe que a coluna da base de dados MTL que contém os valores é a coluna 3, mas verifique quais linhas são das variáveis: RADIANCE_MULT_BAND_b, REFLECTANCE_MULT_BAND_b, RADIANCE_ADD_BAND_b e REFLECTANCE_ADD_BAND_b, lembre-se do type.convert(). Conseguiu escrever o comando? Observe com os chunks abaixo se coincidem, ou se extraem eficientemente os valores desejados:

radiancia <- data.frame(MTL %>% select(3) %>% slice(177:187) %>% type.convert(),
                        MTL %>% select(3) %>% slice(166:176) %>% type.convert())
reflectancia <- data.frame(MTL %>% select(3) %>% slice(197:205) %>% type.convert(),
                           MTL %>% select(3) %>% slice(188:196) %>% type.convert())

Por fim, obtenha K1 e K2 para as bandas 10 e 11. Aqui, faremos a adição de mais colunas, com o objetivo de facilitar uma etapa posterior: um loop para calcular a temperatura de brilho de ambas as bandas. Vale o já dito para os data frames anteriores: verifique em MTL qual a posição da linha que guarda esses valores!

Crie um data frame para K1, e outro para K2. Em ambos, chame a função data.frame(), e nela crie 11 variáveis: uma para cada banda (isso mesmo, de B1 a B11). Em cada uma, iguale a NA, exceto em B10 e B11, nestas extraia os valores correspondentes a K1 (e no outro data frame a K2) de ambas, usando a estrutura já conhecida. Ficou em dúvida? Consulte o chunk abaixo:

K1_10e11 <- data.frame(B1 = NA, B2 = NA, B3 = NA, B4 = NA, B5 = NA, B6 = NA,
                       B7 = NA, B8 = NA, B9 = NA,
                       B10 = MTL %>% select(3) %>% slice(208) %>% type.convert(),
                       B11 = MTL %>% select(3) %>% slice(210) %>% type.convert())
K2_10e11 <- data.frame(B1 = NA, B2 = NA, B3 = NA, B4 = NA, B5 = NA, B6 = NA,
                       B7 = NA, B8 = NA, B9 = NA,
                       B10 = MTL %>% select(3) %>% slice(209) %>% type.convert(),
                       B11 = MTL %>% select(3) %>% slice(211) %>% type.convert())

Por quê acabamos de criar 11 colunas (variáveis) em cada data.frame? Resposta: Lembra que logo acima foi dito que posteriormente usaremos funções lógicas? Então, elas vão facilitar o nosso trabalho em um loop posterior.

Observe o que você criou

Tabela 3 - Fatores offset e gain (1ª e 2ª variável) para radiância.
V3 V3.1
-61.03303 0.0122070
-62.49857 0.0125000
-57.59189 0.0115180
-48.56473 0.0097129
-29.71919 0.0059438
-7.39089 0.0014782
-2.49113 0.0004982
-54.96191 0.0109920
-11.61493 0.0023230
0.10000 0.0003342
0.10000 0.0003342
***
Tabela 4 - Fatores offset e gain (1ª e 2ª variável) para reflectância.
V3 V3.1
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
-0.1 2e-05
***
**Tabela 5 - Constante de conversão térmica K1 para bandas 10 e 11.
B1 B2 B3 B4 B5 B6 B7 B8 B9 V3 V3.1
NA NA NA NA NA NA NA NA NA 774.8853 480.8883
***
**Tabela 6 - Constante de conversão térmica K2 para bandas 10 e 11.
B1 B2 B3 B4 B5 B6 B7 B8 B9 V3 V3.1
NA NA NA NA NA NA NA NA NA 1321.079 1201.144
***
**Tabela 7 - Ângulo de elevação solar para reflectância.
V3
40.01575

Tudo pronto?

Será que temos tudo que precisamos? Vamos conferir então se todas os conjuntos de dados foram criados como esperado: primeiro, confira se a variável type.convert() funcionou em todos os casos, use a função str() da seguinte forma:

str(sun.elevation)
## 'data.frame':    1 obs. of  1 variable:
##  $ V3: num 40
str(radiancia)
## 'data.frame':    11 obs. of  2 variables:
##  $ V3  : num  -61 -62.5 -57.6 -48.6 -29.7 ...
##  $ V3.1: num  0.01221 0.0125 0.01152 0.00971 0.00594 ...
str(reflectancia)
## 'data.frame':    9 obs. of  2 variables:
##  $ V3  : num  -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1
##  $ V3.1: num  2e-05 2e-05 2e-05 2e-05 2e-05 2e-05 2e-05 2e-05 2e-05
str(K1_10e11)
## 'data.frame':    1 obs. of  11 variables:
##  $ B1  : logi NA
##  $ B2  : logi NA
##  $ B3  : logi NA
##  $ B4  : logi NA
##  $ B5  : logi NA
##  $ B6  : logi NA
##  $ B7  : logi NA
##  $ B8  : logi NA
##  $ B9  : logi NA
##  $ V3  : num 775
##  $ V3.1: num 481
str(K2_10e11)
## 'data.frame':    1 obs. of  11 variables:
##  $ B1  : logi NA
##  $ B2  : logi NA
##  $ B3  : logi NA
##  $ B4  : logi NA
##  $ B5  : logi NA
##  $ B6  : logi NA
##  $ B7  : logi NA
##  $ B8  : logi NA
##  $ B9  : logi NA
##  $ V3  : num 1321
##  $ V3.1: num 1201

Se o tipo de dados dos valores for num, então tudo resolvido. Podemos processar as imagens!

3.2 Calculando radiância das bandas

O primeiro processamento que faremos é o cálculo da radiância de todas as bandas espectrais dos sensores OLI e TIRS. Para isso, usaremos a função radconv() do pacote landsat8. Mas antes de fazer essa operação, veremos qual a estrutura desse comando.

O cálculo da radiância é baseado em fatores de conversão de imagem, como dito anteriormente, a relação entre esses fatores é linear, como proposto e utilizado por U.S. Geological Survey (2015), e a equação usada para realizar essa conversão é:

\(L_\lambda = M_L . Q_{cal} + A_L\)

Em que, \(L_\lambda\): radiância espectral (em \(W/(m^2 sr^{-1} \mu m^{-1})\)); \(M_L\): fator de escala multiplicativo para radiância (RADIANCE_MULT_BAND_b); \(A_L\) é o fator de escala aditivo para radiância (RADIANCE_ADD_BAND_b); e \(Q_{cal}\) é o valor do pixel da imagem de nivel 1 (em número digital, DN).

De modo que esses são as variáveis que a função radconv() requer para processar cada imagem desejada. A estrutura dessa função é: radconv(x, Ml, Al), onde x representa o \(Q_{cal}\) da equação acima, mas na função são todos os valores \(Q_{cal}\) de uma imagem.

Agora vamos processar imagem

A radiância é calculada para todas as 11 bandas multiespectrais, e vai de contra as regras de um bom programador, fazer esse procedimento manualmente, por isso, usaremos um loop. Na criação das variáveis foi dito que essa rotina foi pensada para máquinas com memória RAM pequena, e que não seriam criadas variáveis para imagens. Pois bem, aqui faremos isso: a imagem será lida no disco rígido, processada e a imagem de saída será escrita no disco rígido.

O loop (usando a função for (variable in vector)) será escrito criando a variável i, fazendo esta variar no vetor band nos itens 1 a 11 deste. Lembre-se que criamos este vetor logo acima para controlar quais bandas vamos processar.

A função radconv() (bem como as demais funções do pacote landsat8) exige que a imagem (argumento x) esteja em classe de dados do tipo SpatialGridDataFrame, que é uma classe de dados para atributos espaciais que contém uma localização espacial em uma grade regular, ou seja, um destes são imagens georegistradas. A função que lê um arquivo raster como essa classe de dados é a função readGDAL(), do pacote sp.

Usando a função readGDAL() para ler imagens, usemos a função glue() para colar como caractere (sem espaços): o sufixo da banda (a variável produto que criamos), o ‘B’ do nome do produto, a variável ‘i’ do loop for bem como a extensão “.tif” que as bandas possuem no arquivo de origem.

Lida a imagem, podemos aplicar a função readconv(): a função readGDAL(glue(“{produto}B{i}.tif”)) é o primeiro argumento da função readconv(), pois usaremos o %>%! Os demais argumentos são do data frame radiancia, que criamos anteriormente. Como nesse data frame a primeira coluna refere-se ao fator RADIANCE_MULT_BAND_b e a segunda ao RADIANCE_ADD_BAND_b, faça-variar no loop usando a variável ‘i’ para mudar de linha, e pegar o valor correspondente a cada imagem.

Não entendeu até aqui? Observe o chunk abaixo (você já sabe como funciona o %>%!): observe que a variável i da função for varia em band[5:5], para ser didático e mais rápido de executar este código, mas para todas as bandas o correto é configurar: band[1:11].

for (i in band[5:5]) {
  readGDAL(glue("{produto}B{i}.tif")) %>% 
    radconv(radiancia[i,2], radiancia[i,1]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{i}_L.tif"), drivername = "GTiff")
}
## LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

Os últimos argumentos são: a função raster(), usada para transformar o processamento das funções readGDAL() e radconv() de SpatialGridDataFrame em classe RasterLayer, pois, essa classe de variável é exigida pela função writeRaster(); e a função writeRaster(), que fará a exportação do arquivo processado. De modo que os argumentos desta são: filename, para nome do arquivo; e drivername, para formato do arquivo.

O nome do arquivo, foi declarado no chunk usando novamente a função glue, adicionando o sufixo “_L" ao final do nome do arquivo, indicando que trata-se da radiância. Já o formato, será o mesmo de entrada: GTiff.

Observe (no eixo de valores) das plotagens abaixo, a primeira contendo valores em escala de DN (para imagens Landsat 8 de 0 a 65.536) e a segunda em já em \(W/(m^2 sr^{-1} \mu m^{-1})\):

sp::plot(readGDAL("LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif"))
## LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

sp::plot(readGDAL("LC08_L1TP_215073_20190602_20190602_01_T1_B5_L.tif"))
## LC08_L1TP_215073_20190602_20190602_01_T1_B5_L.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

3.3 Calculando reflectância TOA das bandas B1 a B9

O segundo processamento é o cálculo da reflectância no topo da atmosfera (TOA) das bandas espectrais do sensor OLI. Para isso, usaremos a função reflconvS() do pacote landsat8. Vejamos então qual a estrutura desse comando.

O cálculo da reflectância é dividido em duas etapas, a primeira para obter um fator de reflectância espectral no topo da atmosfera (\(\rho\lambda '\)) sem correção de ângulo solar: que também é baseado em fatores de conversão de imagem com relação linear entre fatores (U.S. GEOLOGICAL SURVEY, 2015), e a equação usada para realizar essa conversão é:

\(\rho\lambda ' = M_\rho . Q_{cal} + A_\rho\)

Em que, \(\rho\lambda '\): reflectância espectral no topo da atmosfera sem correção de ângulo solar (adimensional); \(M_\rho\): fator de escala multiplicativo para reflectância (REFLECTANCE_MULT_BAND_b); \(A_\rho\) é o fator de escala aditivo para reflectância (REFLECTANCE_ADD_BAND_b); e \(Q_{cal}\) é o valor do pixel da imagem de nivel 1 (em número digital, DN).

A segunda etapa é a correção do ângulo solar a partir da equação:

\(\rho_\lambda = \rho\lambda ' / sen \theta\)

Em que \(\rho_\lambda\) é a reflectância no topo da atmosfera (adimensional) e \(\theta\) é o ângulo de elevação solar (SUN_ELEVATION), que é único para o momento de captura das imagens.

De modo que esses são as variáveis que a função reflconvS() requer para processar cada imagem desejada. A estrutura dessa função é: reflconvS(x, Mp, Ap, sunelev), onde x representa o \(Q_{cal}\) da equação acima, mas na função são todos os valores \(Q_{cal}\) de uma imagem.

Agora vamos processar imagem

Aqui também iremos usar um loop para processar mais de uma imagem, e a estrutura do loop é semelhante ao cômputo da radiância:

  • Crie um loop com a função for, defina j para variar no vetor band[1:9] (a reflectância é calculada para as bandas do sensor OLI);

  • A estrutura muda apenas na sintaxe da função reflconvS(): lembre-se do data frame reflectancia, com os fatores Mp e Ap, bem como do data frame sun.elevation que contém uma única variável com uma observação, o ângulo de elevação solar:

  • Fixe sun.elevation[1,1]

  • Para o arquivo de saída, mude o sufixo que identifica a imagem já processada para “pTOA”.

Observe o código chunk abaixo para verificar a estrutura! Veja que o vetor band varia de 5 a 5, para fins didáticos, sabendo você que o correto é band[1:9]:

for (j in band[5:5]) {
  readGDAL(glue("{produto}B{j}.tif")) %>% 
    reflconvS(reflectancia[j,2], reflectancia[j,1], sunelev = sun.elevation[1,1]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{j}_pTOA.tif"), drivername = "GTiff")
}
## LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

Observe (no eixo de valores) das plotagens abaixo, a primeira contendo valores em escala de DN (para imagens Landsat 8 de 0 a 65.536), e a segunda em número decimal, o fator de reflectância:

sp::plot(readGDAL("LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif"))
## LC08_L1TP_215073_20190602_20190602_01_T1_B5.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

sp::plot(readGDAL("LC08_L1TP_215073_20190602_20190602_01_T1_B5_pTOA.tif"))
## LC08_L1TP_215073_20190602_20190602_01_T1_B5_pTOA.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

3.4 Calculando temperatura de brilho das bandas B10 e B11

A última etapa de processamento de imagens é a obtenção da temperatura de brilho no topo da atmosfera das bandas 10 e 11 do sensor TIRS. Para isso, usaremos a função tempconv() do pacote landsat8.

O cálculo da temperatura de brilho no tpo da atmosfera é realizado usando uma inversão da equação de Plank descrita pelo U.S. Geological Survey (2015) como:

\(T = K2/ln(K1/L_\lambda + 1)\)

Em que \(T\) é a temperatura de brilho no topo da atmosfera (em Kelvin); \(K1\) é a constante de conversão térmica da banda (K1_CONSTANT_BAND_b); \(K2\) é a constante de conversão térmica da banda (K2_CONSTANT_BAND_b), \(L_\lambda\) é a radiância espectral (em \(W m^{-2} sr^{-1} \mu m^{-1}\)).

E a estrutura da função que calcula a temperatura de brilho é a seguinte: tempcpmv(x, Ml, Al, K1, K2). A função exige Ml e Al pois realiza novamente o cômputo da radiância na sua estrutura. Para K1 e K2 criamos dois data frames para armazenar estes valores.

Agora vamos processar imagem

Aqui também iremos usar um loop para processar mais de uma imagem, e a estrutura do loop é semelhante aos anteriores:

  • Crie um loop com a função for, defina k para variar no vetor band[10:11] (a temperatura de brilho TOA é calculada para as bandas 10 e 11 do sensor TIRS);

  • A estrutura muda apenas na sintaxe da função tempconv(): lembre-se do data frame radiancia, com os fatores Ml e Al, bem como dos data frames K1_10e11 e K2_10e11:

Criamos ambos os data frames (K1_10e11 e K2_10e11) com 11 variáveis para facilitar a atuação da variável k da função for(), de modo que o loop usará os valores 10 e 11 para pegar K1 e K2 dessas duas bandas nestes data frames!

  • Para o arquivo de saída, mude o sufixo que identifica a imagem já processada para “T”.

Observe o código chunk abaixo para verificar a estrutura! Como nesta etapa apenas duas bandas são processadas, o vetor band varia de 10 a 11 (band[10:11]):

for (k in band[10:11]) {
  readGDAL(glue("{produto}B{k}.tif")) %>% 
    tempconv(radiancia[k,2], radiancia[k,1],
             K1 = K1_10e11[1,k],K2=K2_10e11[1,k]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{k}_T.tif"), drivername = "GTiff")
}
## LC08_L1TP_215073_20190602_20190602_01_T1_B10.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns
## LC08_L1TP_215073_20190602_20190602_01_T1_B11.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

Observe no eixo de valores da imagem abaixo que a temperatura no topo da atmosfera predomina em torno de 300 Kelvin (26,85 ºC), valores aceitáveis para uma região tropical de uma imagem do mês de junho:

sp::plot(readGDAL("LC08_L1TP_215073_20190602_20190602_01_T1_B10_T.tif"))
## LC08_L1TP_215073_20190602_20190602_01_T1_B10_T.tif has GDAL driver GTiff 
## and has 7751 rows and 7641 columns

4 Referências

ANDERSON, Cody et al. Radiometric correction of RapidEye imagery using the on-orbit side-slither method. Image And Signal Processing For Remote Sensing Xvii, [s.l.], 6 out. 2011. SPIE. http://dx.doi.org/10.1117/12.898411

U.S. Geological Survey. Landsat 8 (L8): Data users handbook. Sioux Falls: USGS, 2015. 97p.

[Landsat 8 (L8): Data users handbook] (https://www.usgs.gov/land-resources/nli/landsat/landsat-8-data-users-handbook)

À César o que é de César… aos criadores, sua autoria:

citation('landsat8')
## 
## To cite package 'landsat8' in publications use:
## 
##   Alexandre dos Santos (2017). landsat8: Landsat 8 Imagery
##   Rescaled to Reflectance, Radiance and/or Temperature. R package
##   version 0.1-10. https://CRAN.R-project.org/package=landsat8
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {landsat8: Landsat 8 Imagery Rescaled to Reflectance, Radiance and/or
## Temperature},
##     author = {Alexandre dos Santos},
##     year = {2017},
##     note = {R package version 0.1-10},
##     url = {https://CRAN.R-project.org/package=landsat8},
##   }
## 
## ATTENTION: This citation information has been auto-generated from
## the package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.
citation('sp')
## 
## To cite package sp in publications use:
## 
##   Pebesma, E.J., R.S. Bivand, 2005. Classes and methods for
##   spatial data in R. R News 5 (2),
##   https://cran.r-project.org/doc/Rnews/.
## 
##   Roger S. Bivand, Edzer Pebesma, Virgilio Gomez-Rubio, 2013.
##   Applied spatial data analysis with R, Second edition. Springer,
##   NY. http://www.asdar-book.org/
## 
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
citation('raster')
## 
## To cite package 'raster' in publications use:
## 
##   Robert J. Hijmans (2019). raster: Geographic Data Analysis and
##   Modeling. R package version 2.9-5.
##   https://CRAN.R-project.org/package=raster
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {raster: Geographic Data Analysis and Modeling},
##     author = {Robert J. Hijmans},
##     year = {2019},
##     note = {R package version 2.9-5},
##     url = {https://CRAN.R-project.org/package=raster},
##   }
citation('rgdal')
## 
## To cite package 'rgdal' in publications use:
## 
##   Roger Bivand, Tim Keitt and Barry Rowlingson (2019). rgdal:
##   Bindings for the 'Geospatial' Data Abstraction Library. R
##   package version 1.4-3. https://CRAN.R-project.org/package=rgdal
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {rgdal: Bindings for the 'Geospatial' Data Abstraction Library},
##     author = {Roger Bivand and Tim Keitt and Barry Rowlingson},
##     year = {2019},
##     note = {R package version 1.4-3},
##     url = {https://CRAN.R-project.org/package=rgdal},
##   }
citation('dplyr')
## 
## To cite package 'dplyr' in publications use:
## 
##   Hadley Wickham, Romain François, Lionel Henry and Kirill Müller
##   (2019). dplyr: A Grammar of Data Manipulation. R package version
##   0.8.0.1. https://CRAN.R-project.org/package=dplyr
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {dplyr: A Grammar of Data Manipulation},
##     author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller},
##     year = {2019},
##     note = {R package version 0.8.0.1},
##     url = {https://CRAN.R-project.org/package=dplyr},
##   }
citation('glue')
## 
## To cite package 'glue' in publications use:
## 
##   Jim Hester (2018). glue: Interpreted String Literals. R package
##   version 1.3.0. https://CRAN.R-project.org/package=glue
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {glue: Interpreted String Literals},
##     author = {Jim Hester},
##     year = {2018},
##     note = {R package version 1.3.0},
##     url = {https://CRAN.R-project.org/package=glue},
##   }

5 Bonus: R Script para processamento dos dados

Como as etapas anteriores tratam-se de exemplos didáticos, para que o usuário venha a fazer processamento completo de produtos de nível 1 do Landsat 8, segue no chunk abaixo um código inteiro da rotina!

#install.packages("landsat8")
#install.packages("sp")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("dplyr")
#install.packages("glue")
library(landsat8)
library(sp)
library(raster)
library(rgdal)
library(dplyr)
library(glue)
#Definindo pasta das imagens landsat 8 e nome do produto
caminho <- "C:/Users/ERLI/Downloads/LC08_L1TP_215073_20190602_20190602_01_T1/"
setwd(caminho)
produto <- "LC08_L1TP_215073_20190602_20190602_01_T1_"
band <- c(1:11)
#Importando metadados do produto
MTL <- read.table(glue("{produto}MTL.txt"), nrows = 224)
sun.elevation <- MTL %>% select(3) %>% slice(77) %>% type.convert()
radiancia <- data.frame(MTL %>% select(3) %>% slice(177:187) %>% type.convert(),
                        MTL %>% select(3) %>% slice(166:176) %>% type.convert())
reflectancia <- data.frame(MTL %>% select(3) %>% slice(197:205) %>% type.convert(),
                           MTL %>% select(3) %>% slice(188:196) %>% type.convert())
K1_10e11 <- data.frame(B1 = NA, B2 = NA, B3 = NA, B4 = NA, B5 = NA, B6 = NA,
                       B7 = NA, B8 = NA, B9 = NA,
                       B10 = MTL %>% select(3) %>% slice(208) %>% type.convert(),
                       B11 = MTL %>% select(3) %>% slice(210) %>% type.convert())
K2_10e11 <- data.frame(B1 = NA, B2 = NA, B3 = NA, B4 = NA, B5 = NA, B6 = NA,
                       B7 = NA, B8 = NA, B9 = NA,
                       B10 = MTL %>% select(3) %>% slice(209) %>% type.convert(),
                       B11 = MTL %>% select(3) %>% slice(211) %>% type.convert())
#Calculando radiância no topo da atmosfera (TOA)
for (i in band[9:11]) {
  readGDAL(glue("{produto}B{i}.tif")) %>% 
    radconv(radiancia[i,2], radiancia[i,1]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{i}_L.tif"), drivername = "GTiff")
}
#Calculando reflectância no topo da atmosfera (TOA)
for (j in band[1:9]) {
  readGDAL(glue("{produto}B{j}.tif")) %>% 
    reflconvS(reflectancia[j,2], reflectancia[j,1], sunelev = sun.elevation[1,1]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{j}_pTOA.tif"), drivername = "GTiff")
}
#Calculando temperatura de brilho (Kelvin) para as bandas 10 e 11
for (k in band[10:10]) {
  readGDAL(glue("{produto}B{k}.tif")) %>% 
    tempconv(radiancia[k,2], radiancia[k,1], K1 = K1_10e11[1,k],K2=K2_10e11[1,k]) %>% 
    raster() %>% 
    writeRaster(filename = glue("{produto}B{k}_T.tif"), drivername = "GTiff")
}