SMAC Python

Here is the SMAC atmospheric correction routine in Python language. More about SMAC.

Voici la routine de correction atmosphérique SMAC en python En savoir plus sur SMAC

 

#! /usr/bin/env python
# -*- coding: iso-8859-1 -*-
#=============================================================================================
# library for atmospheric correction using SMAC method (Rahman and Dedieu, 1994)
# Contains :
#      smac_inv : inverse smac model for atmospheric correction
#                          TOA==>Surface
#      smac dir : direct smac model
#                          Surface==>TOA
#      coefs : reads smac coefficients
#      PdeZ : #      PdeZ : Atmospheric pressure (in hpa) as a function of altitude (in meters)
 
# Written by O.Hagolle CNES, from the original SMAC C routine
#=============================================================================================
 
from math import *
import numpy as np
 
#=============================================================================================
 
def PdeZ(Z) :
    """
    PdeZ : Atmospheric pressure (in hpa) as a function of altitude (in meters)
 
    """
    p = 1013.25 * pow( 1 - 0.0065 * Z / 288.15 , 5.31 )
    return(p)
#=============================================================================================
 
class coeff:
  def __init__(self,smac_filename):
    with file(smac_filename) as f:
      lines=f.readlines()
    #H20
    temp=lines[0].strip().split()
    self.ah2o=float(temp[0])
    self.nh2o=float(temp[1])
    #O3
    temp=lines[1].strip().split()
    self.ao3=float(temp[0])
    self.no3=float(temp[1])
    #O2
    temp=lines[2].strip().split()
    self.ao2=float(temp[0])
    self.no2=float(temp[1])
    self.po2=float(temp[2])
    #CO2
    temp=lines[3].strip().split()
    self.aco2=float(temp[0])
    self.nco2=float(temp[1])
    self.pco2=float(temp[2])
    #NH4
    temp=lines[4].strip().split()
    self.ach4=float(temp[0])
    self.nch4=float(temp[1])
    self.pch4=float(temp[2])
    #NO2
    temp=lines[5].strip().split()
    self.ano2=float(temp[0])
    self.nno2=float(temp[1])
    self.pno2=float(temp[2])
    #NO2
    temp=lines[6].strip().split()
    self.aco=float(temp[0])
    self.nco=float(temp[1])
    self.pco=float(temp[2])
 
    #rayleigh and aerosol scattering
    temp=lines[7].strip().split()
    self.a0s=float(temp[0])
    self.a1s=float(temp[1])
    self.a2s=float(temp[2])
    self.a3s=float(temp[3])
    temp=lines[8].strip().split()
    self.a0T=float(temp[0])
    self.a1T=float(temp[1])
    self.a2T=float(temp[2])
    self.a3T=float(temp[3])
    temp=lines[9].strip().split()
    self.taur=float(temp[0])
    self.sr=float(temp[0])
    temp=lines[10].strip().split()
    self.a0taup  = float(temp[0])
    self.a1taup  = float(temp[1])
    temp=lines[11].strip().split()
    self.wo      = float(temp[0])
    self.gc      = float(temp[1])
    temp=lines[12].strip().split()
    self.a0P     = float(temp[0])
    self.a1P     = float(temp[1])
    self.a2P     = float(temp[2])
    temp=lines[13].strip().split()
    self.a3P     = float(temp[0])
    self.a4P     = float(temp[1])
    temp=lines[14].strip().split()
    self.Rest1   = float(temp[0])
    self.Rest2   = float(temp[1])
    temp=lines[15].strip().split()
    self.Rest3   = float(temp[0])
    self.Rest4   = float(temp[1])
    temp=lines[16].strip().split()
    self.Resr1   = float(temp[0])
    self.Resr2   = float(temp[1])
    self.Resr3   = float(temp[2])
    temp=lines[17].strip().split()
    self.Resa1   = float(temp[0])
    self.Resa2   = float(temp[1])
    temp=lines[18].strip().split()
    self.Resa3   = float(temp[0])
    self.Resa4   = float(temp[1])
 
#======================================================================
def smac_inv( r_toa, tetas, phis, tetav, phiv,pressure,taup550, uo3, uh2o, coef) :
    """
    r_surf=smac_inv( r_toa, tetas, phis, tetav, phiv,pressure,taup550, uo3, uh2o, coef)
    Corrections atmosphériques
    """
    ah2o    	= coef.ah2o
    nh2o    	= coef.nh2o
    ao3     	= coef.ao3
    no3     	= coef.no3
    ao2     	= coef.ao2
    no2     	= coef.no2
    po2     	= coef.po2
    aco2 	= coef.aco2
    nco2 	= coef.nco2
    pco2 	= coef.pco2
    ach4 	= coef.ach4
    nch4 	= coef.nch4
    pch4 	= coef.pch4
    ano2 	= coef.ano2
    nno2 	= coef.nno2
    pno2 	= coef.pno2
    aco  	= coef.aco
    nco  	= coef.nco
    pco  	= coef.pco
    a0s  	= coef.a0s
    a1s  	= coef.a1s
    a2s  	= coef.a2s
    a3s  	= coef.a3s
    a0T  	= coef.a0T
    a1T  	= coef.a1T
    a2T  	= coef.a2T
    a3T     	= coef.a3T
    taur    	= coef.taur
    sr      	= coef.sr
    a0taup  	= coef.a0taup
    a1taup  	= coef.a1taup
    wo      	= coef.wo
    gc     	= coef.gc
    a0P     	= coef.a0P
    a1P     	= coef.a1P
    a2P     	= coef.a2P
    a3P     	= coef.a3P
    a4P     	= coef.a4P
    Rest1   	= coef.Rest1
    Rest2   	= coef.Rest2
    Rest3   	= coef.Rest3
    Rest4   	= coef.Rest4
    Resr1   	= coef.Resr1
    Resr2   	= coef.Resr2
    Resr3   	= coef.Resr3
    Resa1   	= coef.Resa1
    Resa2   	= coef.Resa2
    Resa3   	= coef.Resa3
    Resa4   	= coef.Resa4
 
    cdr=pi/180
    crd=180/pi
 
    #/*------:  calcul de la reflectance de surface  smac               :--------*/
 
    us = cos (tetas * cdr)
    uv = cos (tetav * cdr)
    Peq= pressure/1013.25 
 
    #/*------:  1) air mass */
    m =  1/us + 1/uv
 
    #/*------:  2) aerosol optical depth in the spectral band, taup     :--------*/
    taup = (a0taup) + (a1taup) * taup550 
 
    #/*------:  3) gaseous transmissions (downward and upward paths)    :--------*/
    to3 = 1.
    th2o= 1.
    to2 = 1.
    tco2= 1.
    tch4= 1. 
 
    uo2 =  (Peq ** (po2))
    uco2=  (Peq ** (pco2))
    uch4=  (Peq ** (pch4))
    uno2=  (Peq ** (pno2))
    uco =  (Peq ** (pco))
 
    #/*------:  4) if uh2o <= 0 and uo3 <=0 no gaseous absorption is computed  :--------*/
    to3   = exp ( (ao3)  * ( (uo3 *m)  ** (no3)  ) )
    th2o  = exp ( (ah2o) * ( (uh2o*m)  ** (nh2o) ) )
    to2   = exp ( (ao2)  * ( (uo2 *m)  ** (no2)  ) )
    tco2  = exp ( (aco2) * ( (uco2*m)  ** (nco2) ) )
    tch4  = exp ( (ach4) * ( (uch4*m)  ** (nch4) ) )
    tno2  = exp ( (ano2) * ( (uno2*m)  ** (nno2) ) )
    tco   = exp ( (aco)  * ( (uco *m)  ** (nco) ) )
    tg      = th2o * to3 * to2 * tco2 * tch4 * tco * tno2
 
    #/*------:  5) Total scattering transmission                      :--------*/
    ttetas = (a0T) + (a1T)*taup550/us + ((a2T)*Peq + (a3T))/(1.+us)  #/* downward */
    ttetav = (a0T) + (a1T)*taup550/uv + ((a2T)*Peq + (a3T))/(1.+uv)  #/* upward   */
 
    #/*------:  6) spherical albedo of the atmosphere                 :--------*/
    s = (a0s) * Peq +  (a3s) + (a1s)*taup550 + (a2s) * (taup550 ** 2) 
 
    #/*------:  7) scattering angle cosine                            :--------*/
    cksi = - ( (us*uv) + (sqrt(1. - us*us) * sqrt (1. - uv*uv)*cos((phis-phiv) * cdr) ) )
    if (cksi < -1 ) :
	   cksi=-1.0 
 
    #/*------:  8) scattering angle in degree 			 :--------*/
    ksiD = crd*acos(cksi) 
 
    #/*------:  9) rayleigh atmospheric reflectance 			 :--------*/
    ray_phase = 0.7190443 * (1. + (cksi*cksi))  + 0.0412742
    ray_ref   = ( taur*ray_phase ) / (4*us*uv)
    ray_ref   = ray_ref*pressure / 1013.25
    taurz=(taur)*Peq
 
    #/*------:  10) Residu Rayleigh 					 :--------*/
    Res_ray= Resr1 + Resr2 * taur*ray_phase / (us*uv) + Resr3 * ( (taur*ray_phase/(us*uv)) ** 2) 
 
    #/*------:  11) aerosol atmospheric reflectance			 :--------*/
    aer_phase = a0P + a1P*ksiD + a2P*ksiD*ksiD +a3P*(ksiD ** 3) + a4P * (ksiD ** 4)
 
    ak2 = (1. - wo)*(3. - wo*3*gc)
    ak  = sqrt(ak2)
    e   = -3*us*us*wo /  (4*(1. - ak2*us*us) )
    f   = -(1. - wo)*3*gc*us*us*wo / (4*(1. - ak2*us*us) )
    dp  = e / (3*us) + us*f
    d   = e + f
    b   = 2*ak / (3. - wo*3*gc)
    delta = np.exp( ak*taup )*(1. + b)*(1. + b) - np.exp(-ak*taup)*(1. - b)*(1. - b)
    ww  = wo/4.
    ss  = us / (1. - ak2*us*us)
    q1  = 2. + 3*us + (1. - wo)*3*gc*us*(1. + 2*us)
    q2  = 2. - 3*us - (1. - wo)*3*gc*us*(1. - 2*us)
    q3  = q2*np.exp( -taup/us )
    c1  =  ((ww*ss) / delta) * ( q1*np.exp(ak*taup)*(1. + b) + q3*(1. - b) )
    c2  = -((ww*ss) / delta) * (q1*np.exp(-ak*taup)*(1. - b) + q3*(1. + b) )
    cp1 =  c1*ak / ( 3. - wo*3*gc )
    cp2 = -c2*ak / ( 3. - wo*3*gc )
    z   = d - wo*3*gc*uv*dp + wo*aer_phase/4.
    x   = c1 - wo*3*gc*uv*cp1
    y   = c2 - wo*3*gc*uv*cp2
    aa1 = uv / (1. + ak*uv)
    aa2 = uv / (1. - ak*uv)
    aa3 = us*uv / (us + uv) 
 
    aer_ref = x*aa1* (1. - np.exp( -taup/aa1 ) )
    aer_ref = aer_ref + y*aa2*( 1. - np.exp( -taup / aa2 )  )
    aer_ref = aer_ref + z*aa3*( 1. - np.exp( -taup / aa3 )  )
    aer_ref = aer_ref / ( us*uv )
 
    #/*------:  12) Residu Aerosol  					:--------*/
    Res_aer= ( Resa1 + Resa2 * ( taup * m *cksi ) + Resa3 * ( (taup*m*cksi ) ** 2) )+ Resa4 * ( (taup*m*cksi) ** 3)
 
    #/*------:  13)  Terme de couplage molecule / aerosol		:--------*/
    tautot=taup+taurz
    Res_6s= ( Rest1 + Rest2 * ( tautot * m *cksi )   + Rest3 * ( (tautot*m*cksi) ** 2) ) + Rest4 * ( (tautot*m*cksi) ** 3)
 
    #/*------:  14) total atmospheric reflectance  			:--------*/
    atm_ref = ray_ref - Res_ray + aer_ref - Res_aer + Res_6s
 
    #/*------:  15) Surface reflectance  				:--------*/
 
    r_surf = r_toa - (atm_ref * tg)
    r_surf = r_surf / ( (tg * ttetas * ttetav) + (r_surf * s) ) 
 
    return r_surf
 
#=======================================================================================================
def smac_dir ( r_surf, tetas, phis, tetav, phiv,pressure,taup550, uo3, uh2o, coef):
    """
    r_toa=smac_dir ( r_surf, tetas, phis, tetav, phiv,pressure,taup550, uo3, uh2o, coef)
    Application des effets atmosphériques
    """
 
    ah2o    	= coef.ah2o
    nh2o    	= coef.nh2o
    ao3     	= coef.ao3
    no3     	= coef.no3
    ao2     	= coef.ao2
    no2     	= coef.no2
    po2     	= coef.po2
    aco2 	= coef.aco2
    nco2 	= coef.nco2
    pco2 	= coef.pco2
    ach4 	= coef.ach4
    nch4 	= coef.nch4
    pch4 	= coef.pch4
    ano2 	= coef.ano2
    nno2 	= coef.nno2
    pno2 	= coef.pno2
    aco  	= coef.aco
    nco  	= coef.nco
    pco  	= coef.pco
    a0s  	= coef.a0s
    a1s  	= coef.a1s
    a2s  	= coef.a2s
    a3s  	= coef.a3s
    a0T  	= coef.a0T
    a1T  	= coef.a1T
    a2T  	= coef.a2T
    a3T     	= coef.a3T
    taur    	= coef.taur
    sr      	= coef.sr
    a0taup  	= coef.a0taup
    a1taup  	= coef.a1taup
    wo      	= coef.wo
    gc     	= coef.gc
    a0P     	= coef.a0P
    a1P     	= coef.a1P
    a2P     	= coef.a2P
    a3P     	= coef.a3P
    a4P     	= coef.a4P
    Rest1   	= coef.Rest1
    Rest2   	= coef.Rest2
    Rest3   	= coef.Rest3
    Rest4   	= coef.Rest4
    Resr1   	= coef.Resr1
    Resr2   	= coef.Resr2
    Resr3   	= coef.Resr3
    Resa1   	= coef.Resa1
    Resa2   	= coef.Resa2
    Resa3   	= coef.Resa3
    Resa4   	= coef.Resa4
 
    cdr=pi/180
    crd=180/pi
 
    #/*------:  calcul de la reflectance de surface  smac               :--------*/
 
    us = cos (tetas * cdr)
    uv = cos (tetav * cdr)
    Peq= pressure/1013.25 
 
    #/*------:  1) air mass */
    m =  1/us + 1/uv
 
    #/*------:  2) aerosol optical depth in the spectral band, taup     :--------*/
    taup = (a0taup) + (a1taup) * taup550 
 
    #/*------:  3) gaseous transmissions (downward and upward paths)    :--------*/
    to3 = 1.
    th2o= 1.
    to2 = 1.
    tco2= 1.
    tch4= 1. 
 
    uo2 =  (Peq ** (po2))
    uco2=  (Peq ** (pco2))
    uch4=  (Peq ** (pch4))
    uno2=  (Peq ** (pno2))
    uco =  (Peq ** (pco))
 
    #/*------:  4) if uh2o <= 0 and uo3<= 0 no gaseous absorption is computed  :--------*/
    to3   = exp ( (ao3)  * ( (uo3 *m)  ** (no3)  ) )
    th2o  = exp ( (ah2o) * ( (uh2o*m)  ** (nh2o) ) )
    to2   = exp ( (ao2)  * ( (uo2 *m)  ** (no2)  ) )
    tco2  = exp ( (aco2) * ( (uco2*m)  ** (nco2) ) )
    tch4  = exp ( (ach4) * ( (uch4*m)  ** (nch4) ) )
    tno2  = exp ( (ano2) * ( (uno2*m)  ** (nno2) ) )
    tco   = exp ( (aco)  * ( (uco *m)  ** (nco) ) )
    tg      = th2o * to3 * to2 * tco2 * tch4 * tco * tno2
 
    #/*------:  5) Total scattering transmission                      :--------*/
    ttetas = (a0T) + (a1T)*taup550/us + ((a2T)*Peq + (a3T))/(1.+us)  #/* downward */
    ttetav = (a0T) + (a1T)*taup550/uv + ((a2T)*Peq + (a3T))/(1.+uv)  #/* upward   */
 
    #/*------:  6) spherical albedo of the atmosphere                 :--------*/
    s = (a0s) * Peq +  (a3s) + (a1s)*taup550 + (a2s) * (taup550 ** 2) 
 
    #/*------:  7) scattering angle cosine                            :--------*/
    cksi = - ( (us*uv) + (sqrt(1. - us*us) * sqrt (1. - uv*uv)*cos((phis-phiv-360) * cdr) ) )
    if (cksi < -1 ) :
       cksi=-1.0 
 
    #/*------:  8) scattering angle in degree            :--------*/
    ksiD = crd*acos(cksi) 
 
    #/*------:  9) rayleigh atmospheric reflectance              :--------*/
    ray_phase = 0.7190443 * (1. + (cksi*cksi))  + 0.0412742
    ray_ref   = ( taur*ray_phase ) / (4*us*uv)
    ray_ref   = ray_ref*pressure / 1013.25
    taurz=(taur)*Peq
 
    #/*------:  10) Residu Rayleigh                      :--------*/
    Res_ray= Resr1 + Resr2 * taur*ray_phase / (us*uv) + Resr3 * ( (taur*ray_phase/(us*uv)) ** 2) 
 
    #/*------:  11) aerosol atmospheric reflectance          :--------*/
    aer_phase = a0P + a1P*ksiD + a2P*ksiD*ksiD +a3P*(ksiD ** 3) + a4P * (ksiD ** 4)
 
    ak2 = (1. - wo)*(3. - wo*3*gc)
    ak  = sqrt(ak2)
    e   = -3*us*us*wo /  (4*(1. - ak2*us*us) )
    f   = -(1. - wo)*3*gc*us*us*wo / (4*(1. - ak2*us*us) )
    dp  = e / (3*us) + us*f
    d   = e + f
    b   = 2*ak / (3. - wo*3*gc)
    delta = np.exp( ak*taup )*(1. + b)*(1. + b) - np.exp(-ak*taup)*(1. - b)*(1. - b)
    ww  = wo/4.
    ss  = us / (1. - ak2*us*us)
    q1  = 2. + 3*us + (1. - wo)*3*gc*us*(1. + 2*us)
    q2  = 2. - 3*us - (1. - wo)*3*gc*us*(1. - 2*us)
    q3  = q2*np.exp( -taup/us )
    c1  =  ((ww*ss) / delta) * ( q1*np.exp(ak*taup)*(1. + b) + q3*(1. - b) )
    c2  = -((ww*ss) / delta) * (q1*np.exp(-ak*taup)*(1. - b) + q3*(1. + b) )
    cp1 =  c1*ak / ( 3. - wo*3*gc )
    cp2 = -c2*ak / ( 3. - wo*3*gc )
    z   = d - wo*3*gc*uv*dp + wo*aer_phase/4.
    x   = c1 - wo*3*gc*uv*cp1
    y   = c2 - wo*3*gc*uv*cp2
    aa1 = uv / (1. + ak*uv)
    aa2 = uv / (1. - ak*uv)
    aa3 = us*uv / (us + uv) 
 
    aer_ref = x*aa1* (1. - np.exp( -taup/aa1 ) )
    aer_ref = aer_ref + y*aa2*( 1. - np.exp( -taup / aa2 )  )
    aer_ref = aer_ref + z*aa3*( 1. - np.exp( -taup / aa3 )  )
    aer_ref = aer_ref / ( us*uv )
 
    #/*------:  12) Residu Aerosol                      :--------*/
    Res_aer= ( Resa1 + Resa2 * ( taup * m *cksi ) + Resa3 * ( (taup*m*cksi ) ** 2) )+ Resa4 * ( (taup*m*cksi) ** 3)
 
    #/*------:  13)  Terme de couplage molecule / aerosol       :--------*/
    tautot=taup+taurz
    Res_6s= ( Rest1 + Rest2 * ( tautot * m *cksi )   + Rest3 * ( (tautot*m*cksi) ** 2) ) + Rest4 * ( (tautot*m*cksi) ** 3)
 
    #/*------:  14) total atmospheric reflectance           :--------*/
    atm_ref = ray_ref - Res_ray + aer_ref - Res_aer + Res_6s
 
    #/*------:  15) TOA reflectance                 :--------*/
 
    r_toa = r_surf*tg*ttetas*ttetav/(1-r_surf*s) + (atm_ref * tg) 
 
    return r_toa
 
#=============================================================================
if __name__=="__main__" :
    #example
    theta_s = 45
    theta_v = 5
    phi_s =  200
    phi_v = -160
    r_toa=0.2
    ######################################lecture des coefs_smac
    nom_smac ='COEFS/coef_FORMOSAT2_B1_CONT.dat'
    coefs=coeff(nom_smac)
    bd=1
    r_surf = smac_inv(r_toa , theta_s, phi_s, theta_v, phi_v,1013,0.1,0.3,0.3, coefs)
    r_toa2 = smac_dir(r_surf, theta_s, phi_s, theta_v, phi_v,1013,0.1,0.3,0.3, coefs)       
 
    print r_toa, r_surf,r_toa2

5 comments

  • Michon on 20/05/2014 at 20:15 said:

    Bonjour,
    voilà je voudrais utilisée la routine SMAC pour corriger atmosphériquemment des images landsat 7 mais je n'arrive à intégrer la routine smac sous python . est-ce que quelqu'un pourrait m'aider je débute en programmation?
    merci d'avance
    Stéphanie

  • G. Bippus on 28/01/2016 at 16:17 said:

    Hello,
    thank you very much for the code!
    I found a difference in the Python code compared with the C-code of SMAC available at the CESBIO website:
    in Step 10 - "Residual Rayleigh", in the Python code the factor "taur" is used as provided in the coefficient files, while in the C-code this factor is weighted in Step 9 - "Rayleigh atmospheric reflectance" with Peq (taurz = taur * Peq).

    Could you please clarify if this is a typo in the Python code, or if this change was made consciously compared to the C-code?
    In the latter case, why?

    Thank you very much in advance,
    Gabi

    • Olivier Hagolle on 29/01/2016 at 16:32 said:

      Dear Gabriele,
      thanks for pointing this inconsistency !
      I do not have the answer yet, we have dug in all the multiple versions we have of the code, and some of them use taurz, others taur, we thus need to dig a little more to find out which one is the right one.
      As it is a regression, the answer should be in the software which computes the coefficients, but we have several versions of it...

      I tested the effect of the error, for a blue band and a pressure of 900 hpa. Difference is only 0.35%, but still, we should correct it.
      I'll come back to you soon.
      Olivier

      • G. Bippus on 21/04/2016 at 14:46 said:

        Dear Olivier,
        did you already had time to clarify the inconsistencies in the C- and Python-codes of SMAC regarding the parameters taur and taurz, respectively?

        Further, do you have any experiences with combining the atmospheric correction with a topographic correction?

        Thanks!
        Gabi

        • Olivier Hagolle on 13/05/2016 at 06:40 said:

          Dear Gabriele,
          sorry for replying so late, last weeks have been hectic.
          We did not have time to look at your question yet. Please remind us from time to time !
          Olivier

Leave a Reply to Michon Cancel reply

Your email address will not be published. Required fields are marked *

*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>