lst_based_on_sc_and_isc.py

# -*- coding: utf-8 -*-
"""LST based on SC and ISC.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1FQv8F05UnisdC0UIfUBEKuym-bD7hQxk
"""

import os
import numpy as np
import imageio
import imageio as im
import gdal
import rasterio 
from glob import glob
import matplotlib.pyplot as plt

"""### IMPORTANDO VARIÁVEIS DO METADADO PARA CALIBRAÇÃO RADIOMÉTRICA E CONVERSÃO TÉRMICA (TbTOA) ###"""

## 1 CONSTANTES PARA CALIBRAÇÃO RADIOMÉTRICA (transformando de ND para Radiância TOA em Wm-2sr-1um-1)

## 1.1 Extraindo fator de escala multiplicativo da banda a partir do metadata da cena (MULT)


MULT = 'MULT_BAND_10'
def encontrar(MULT):
    with open('caminho_metadata.txt', 'r') as MTL:
        linhas = MTL.readlines()
        for i in range(len(linhas)):
            if 'MULT_BAND_10' in linhas[i]:
                a = linhas[i].split('=')[-1].strip('\n')
                
    return a
ML = encontrar(MULT)

################################################################################################################################

## 1.2 Extraindo fator de escala aditivo da banda a partir do metadata da cena (ADD)


ADD = 'ADD_BAND_10'
def encontrar(ADD):
    with open('caminho_metadata.txt', 'r') as MTL:
        linhas = MTL.readlines()
        for i in range(len(linhas)):
            if 'ADD_BAND_10' in linhas[i]:
                b = linhas[i].split('=')[-1].strip('\n')
                
    return b
AL = encontrar(ADD)

## 2 CONSTANTES DE CONVERSÃO TERMAL PARA CALCULDO DA TEMPERATURA DE BRILHO (transformando de RAD TOA para TB TOA em K)

## 2.1 Extraindo constante de conversão termal da banda a partir do metadata da cena (K1)


K1 = 'K1_CONSTANT_BAND_10'
def encontrar(K1):
    with open('caminho_metadata.txt', 'r') as MTL:
        linhas = MTL.readlines()
        for i in range(len(linhas)):
            if 'K1_CONSTANT_BAND_10' in linhas[i]:
                b = linhas[i].split('=')[-1].strip('\n')
                
    return b
K1 = encontrar(K1)

################################################################################################################################

## 2.2 Extraindo constante de conversão termal da banda a partir do metadata da cena (K2)

K2 = 'K2_CONSTANT_BAND_10'
def encontrar(K2):
    with open('caminho_metadata.txt', 'r') as MTL:
        linhas = MTL.readlines()
        for i in range(len(linhas)):
            if 'K2_CONSTANT_BAND_10' in linhas[i]:
                b = linhas[i].split('=')[-1].strip('\n')
                
    return b
K2 = encontrar(K2)

"""### IMPORTANDO BANDA THERMAL EM ND ###"""

banda10 = glob('caminho_da_banda.TIF')  #inserir o diretório no qual a banda se encontra
with rasterio.open(banda10[0]) as src:
    b10= src.read(1)
b_10 = b10.astype(float)
kwargs = src.meta.copy()
kwargs.update({'dtype': 'float64'})

"""### EXECUTANDO A CALIBRAÇÃO RADIOMÉTRICA ###"""

ml = float(ML) #fator de escala multiplicativo da banda
al = float(AL) #fator de escala aditivo da banda

#função para calibração radiométrica

def radcal(b_10):
    for i in b_10:
        radiance = (ml*b_10) + al
    return radiance

Lsens = radcal(b_10).astype(float) #obtendo Radiância TOA em Wm-2sr-1um-1

with rasterio.open('Radiance_TOA_BX.tiff', 'w', **kwargs) as dst: #exportando imagem final com calibração radiométrica 
    dst.write(Lsens.astype(rasterio.float64),1)

"""### EXECUTANDO CONVERSÃO DE RADIANCIA TOA PARA TEMPERATURA DE BRILHO TOA ###"""

k1 = float(K1) #constante de conversão termal
k2 = float(K2) #constante de conversão termal

#função para conversão de rad TOA para tbTOA

def TbTOA(Lsens):
    for i in Lsens:
        TbTOA = k2/(np.log((k1/Lsens)+1))
    return TbTOA

TbTOA = TbTOA(Lsens).astype(float) #obtendo TbTOA em K

with rasterio.open('Tb_TOA_BX.tiff', 'w', **kwargs) as dst: #exportando imagem final convertida para TbTOA 
    dst.write(TbTOA.astype(rasterio.float64),1)

"""### DEFINIÇÃO DAS VARIÁVEIS DA SUPERFÍCIE E ATMOSFÉRICAS ###"""

#Definição de variáveis conhecidas para método SC e ISC#

#Cálculo Vapor d'água método de Leckner (1978)

#importando dados do perfil atmosférico 
a = open('vertical_profile.txt', 'r') 
VP = a.readlines()

#pressão da primeira camada/layer/primeira altitude do perfil vertical atmosférico
pressure = float(VP[5][1:7]) 

#temperatura próxima à superfície; temperatura da primeira camada do perfil vertical (K)
To = float(VP[5][17:21]) + 273.15 

#umidade relativa da primeira camada do perfil vertical, em fração
RH = float(VP[5][33:35])/100 

###################################################################

#APLICAÇÃO CÁLCULO VAPOR D'ÁGUA POR LECKNER#

#Cálculo da pressão parcial (Ps) de vapor d'água
Ps = np.exp(26.23-(5416/To)) 


#cálculo do vapor d'água em gcm-2
w = (0.493*RH*Ps)/To  


#Definição da temperatura atmosférica média (Ta) para método ISC, equação conforme Qin et al. (2001)
Ta = 16.011 + (0.9262*To) #para midlatsummer

###################################################################

#Variável da superfície 

emis = 0.9733 #emissividade no comprimento de onda efetivo

#Comprimento de onda efetivo da banda

lamb = 10.904 #comprimento de onda efetivo para banda 10 sensor TIR LANDSAT 8 (Jiménez-Muñoz, 2014)

#constantes da equação de Planck
c1= 1.19104*(10**8) #Wµm4m-2sr-1 
c2= 14387.7 #µmK

##################################################################

#Constantes derivadas de Planck para método SC e ISC #

#######GAMA######
gama = (((c2*Lsens)/(TbTOA**2))*((((lamb**4)/c1)*Lsens)+(1/lamb)))**-1

#######SIGMA######
sigma = -(gama*Lsens) + TbTOA

"""### RECUPERAÇÃO TEMPERATURA DA SUPERFÍCIE TERRESTRE (TST) PELA MÉTODO SINGLE CHANNEL GENERALIZADO (SCG) ###"""

#Função para recuperar TST pelo método SCG

def tempSC(Lsens):
    # DEFINIÇÃO VARÍAVEIS ATMOSFÉRICAS DENTRO DO MÉTODO #
    psi1 = 0.04019*w*w + 0.02916*w + 1.01523
    psi2 = -0.38333*w*w - 1.50294*w + 0.20324
    psi3 = 0.00918*w*w + 1.36072*w - 0.27514
    
    for i in Lsens:
        tempSC = (gama*((((psi1*Lsens) + psi2)/emis) + psi3)) + sigma
    return tempSC

tempSC = tempSC(Lsens).astype(float) #obtendo TST pelo SCG em Kelvin

with rasterio.open('TST_SCG_BX.tiff', 'w', **kwargs) as dst: #exportando imagem final de TST pelo método SCG
    dst.write(tempSC.astype(rasterio.float64),1)

"""### RECUPERAÇÃO TEMPERATURA DA SUPERFÍCIE TERRESTRE (TST) PELA MÉTODO IMPROVED SINGLE CHANNEL (ISC) ###"""

#Função para recuperar TST pelo método ISC
def tempISC(Lsens):
    
    # DEFINIÇÃO VARÍAVEIS ATMOSFÉRICAS DENTRO DO MÉTODO #
    psi1 = -7.2121979375*w*w+0.0000492124*Ta*Ta-2.4523205637*w-0.0262745276*Ta-0.0000496173*Ta*Ta*w+0.0231691781*Ta*w+0.0466282124*Ta*w*w-0.0000748260*Ta*Ta*w*w+4.4729730361
    psi2 = 89.6156888857*w*w-0.0003760208*Ta*Ta+106.5509303783*w+0.2157797227*Ta+0.0014080695*Ta*Ta*w-0.7844419527*Ta*w-0.5731956714*Ta*w*w+0.0009118768*Ta*Ta*w*w-30.3702785256
    psi3 = -14.6595491055*w*w-0.0001047275*Ta*Ta-79.9583806096*w+0.0418090158*Ta-0.0009095018*Ta*Ta*w+0.5453487543*Ta*w+0.0911362208*Ta*w*w-0.0001417749*Ta*Ta*w*w-3.7618398628
         
    for i in Lsens:
        tempISC = (gama*((((psi1*Lsens) + psi2)/emis) + psi3)) + sigma
    return tempISC

tempISC = tempISC(Lsens).astype(float) #obtendo TST pelo ISC em Kelvin

with rasterio.open('TST_ISC_BX.tiff', 'w', **kwargs) as dst: #exportando imagem final de TST pelo método ISC
    dst.write(tempISC.astype(rasterio.float64),1)