KdJametPython/class_KdJamet.py

# -*- coding: utf-8 -*-
# class for Jamet Neural Network Kd inversion.
# Python version :
# 2016-05-11 D.Dessailly david.dessailly@univ-littoral.fr
# LOG (Laboratoire d'Oceanoligie et Geoscience)
import numpy as np
import sys
import os

Sensor = {'idSEAWIFS' : 0, 'idMODIS' : 1, 'idMERIS' : 2}
NBWL = 6

#
# Class to handle Jamet's Kd processing for MERIS/OLCI, SEAWIFS and MODIS Sensor
#
class KdJamet:
  # -- default __init__ Method +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Input :
  # idSensor : valid values -> ['idSEAWIFS', 'idMODIS', 'idMERIS']
  # --
  def __init__(self, idSensor):
    self.sensor = idSensor
    try:
      self.lutpath = os.getenv('KD_LUTS_PATH')
    except:
      print "KD_LUTS_PATH not set, add it to your environment variables\n for bash user add the following line to your .bashrc file :\n export KD_LUTS_PATH=/where/are/my/KD/LUTS\n then:\n source ~/.bashrc"
      sys.exit()
    
    if idSensor in Sensor:
      # 
      self.init_variables()
      try:
        self.read_LUTs_Kd()
      except:
        print "Failed to read %s LUTs" % (idSensor[2:])
        raise
        sys.exit()
    else:
      print "invalid Sensor ID : %s\n valid values : idSEAWIFS, idMODIS, idMERIS" % (idSensor)
      sys.exit()
  
  # -- METHOD init_variables ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Sensor dependants variables and Perceptron parameters setting 
  # --
  def init_variables(self):
    print "class KdJamet.init_variables Method"
    # =======================================================================================================
    # MERIS / OLCI  
    if self.sensor == 'idMERIS':
      self.WL = np.array([413, 443, 490, 510, 560, 665])
      # Variables definiton for network Rrs490/Rrs555 >= .85 -------------------------------------------------
      # hidden layer number 
      self.rsup_NC=2
      # input data number
      self.rsup_NE=5
      # number of neuron in the first hidden layer 
      self.rsup_NC1=7
      # number of neuron in the second hidden layer
      self.rsup_NC2=4
      # nb output
      self.rsup_NS=1
      # nb input + nb output*/
      self.rsup_NES=6
      # LUTs file names 
      self.rsup_LUT_POIDS = "KdOLCI_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_555_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_7_4_switch_110416.sn"
      self.rsup_LUT_MOY = "Moy_KdOLCI_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_ss_620.dat"
      self.rsup_LUT_ECART = "Ecart_KdOLCI_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_ss_620.dat"

      # Variables definiton for network < .85 -----------------------------------
      # nombre de couches cachees - hidden layer number
      self.rinf_NC=2
      # input data number*/
      self.rinf_NE=6
      # number of neuron in the first hidden layer
      self.rinf_NC1=5
      # number of neuron in the second hidden layer
      self.rinf_NC2=5
      # nb output
      self.rinf_NS=1
      # nb input + nb output
      self.rinf_NES=7
      #LUTs file names
      self.rinf_LUT_POIDS = "KdOLCI_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_670_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_5_5_switch_110416.sn"
      self.rinf_LUT_MOY = "Moy_KdOLCI_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_ss_620.dat"
      self.rinf_LUT_ECART = "Ecart_KdOLCI_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_ss_620.dat"
      
    # =======================================================================================================
    # MODIS
    if self.sensor == 'idMODIS':
      self.WL = np.array([412, 443, 488, 531, 547, 667])
      # Variables definiton for network Rrs490/Rrs555 >= .85 -------------------------------------------------
      # hidden layer number 
      self.rsup_NC=2
      # input data number
      self.rsup_NE=5
      # number of neuron in the first hidden layer 
      self.rsup_NC1=5
      # number of neuron in the second hidden layer
      self.rsup_NC2=4
      # nb output
      self.rsup_NS=1
      # nb input + nb output*/
      self.rsup_NES=6
      # LUTs file names 
      self.rsup_LUT_POIDS = "KdMODIS_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_555_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_5_4_publi_CSIRO.sn"
      self.rsup_LUT_MOY = "Moy_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_555_log_Kd_lambda_merge_seuil_15_ss_645_publi_sup_CSIRO.dat"
      self.rsup_LUT_ECART = "Ecart_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_555_log_Kd_lambda_merge_seuil_15_ss_645_publi_sup_CSIRO.dat"

      # Variables definiton for network < .85 -----------------------------------
      # nombre de couches cachees - hidden layer number
      self.rinf_NC=2
      # input data number*/
      self.rinf_NE=6
      # number of neuron in the first hidden layer
      self.rinf_NC1=10
      # number of neuron in the second hidden layer
      self.rinf_NC2=7
      # nb output
      self.rinf_NS=1
      # nb input + nb output
      self.rinf_NES=7
      #LUTs file names
      self.rinf_LUT_POIDS = "KdMODIS_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_670_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_10_7_publi_CSIRO.sn"
      self.rinf_LUT_MOY = "Moy_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_ss_645_publi_inf_CSIRO.dat"
      self.rinf_LUT_ECART = "Ecart_KdMODIS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_ss_645_publi_inf_CSIRO.dat"
      
    # =======================================================================================================
    # SEAWIFS
    if self.sensor == 'idSEAWIFS':
      self.WL = np.array([412, 443, 490, 510, 555, 670])
    # Variables definiton for network Rrs490/Rrs555 >= .85 -------------------------------------------------
      # hidden layer number 
      self.rsup_NC=2
      # input data number
      self.rsup_NE=5
      # number of neuron in the first hidden layer 
      self.rsup_NC1=4
      # number of neuron in the second hidden layer
      self.rsup_NC2=4
      # nb output
      self.rsup_NS=1
      # nb input + nb output*/
      self.rsup_NES=6
      # LUTs file names 
      self.rsup_LUT_POIDS = "KdSeaWiFS_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_555_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_4_4_publi_200715.sn"
      self.rsup_LUT_MOY = "Moy_KdSeaWiFS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_160715.dat"
      self.rsup_LUT_ECART = "Ecart_KdSeaWiFS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_160715.dat"

      # Variables definiton for network < .85 -----------------------------------
      # nombre de couches cachees - hidden layer number
      self.rinf_NC=2
      # input data number*/
      self.rinf_NE=6
      # number of neuron in the first hidden layer
      self.rinf_NC1=5
      # number of neuron in the second hidden layer
      self.rinf_NC2=5
      # nb output
      self.rinf_NS=1
      # nb input + nb output
      self.rinf_NES=7
      #LUTs file names
      self.rinf_LUT_POIDS = "KdSeaWiFS_poids_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_443_to_670_log_Kd_lambda_merge_seuil_15_angle_ascii_6x1_hh_5_5_publi_160715_ter.sn"
      self.rinf_LUT_MOY = "Moy_KdSeaWiFS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_160715.dat"
      self.rinf_LUT_ECART = "Ecart_KdSeaWiFS_IOCCG_NOMAD_BOUM_ss_12_COASTCOLOUR_412_to_670_log_Kd_lambda_merge_seuil_15_angle_publi_160715.dat"
      
    # LUTs arrays definition
    self.rsup_b1 = np.zeros((self.rsup_NC1), dtype=np.float64) 
    self.rsup_b2 = np.zeros((self.rsup_NC2), dtype=np.float64) 
    self.rsup_b3 = 0.
    self.rsup_w1 = np.zeros((self.rsup_NE,self.rsup_NC1), dtype=np.float64)
    self.rsup_w2 = np.zeros((self.rsup_NC1,self.rsup_NC2), dtype=np.float64)
    self.rsup_w3 = np.zeros((self.rsup_NC2), dtype=np.float64)
    self.rsup_moy = np.zeros((self.rsup_NES), dtype=np.float64)
    self.rsup_ecart = np.zeros((self.rsup_NES), dtype=np.float64)

    self.rinf_b1 = np.zeros((self.rinf_NC1), dtype=np.float64)
    self.rinf_b2 = np.zeros((self.rinf_NC2), dtype=np.float64)
    self.rinf_b3 = 0.
    self.rinf_w1 = np.zeros((self.rinf_NE,self.rinf_NC1), dtype=np.float64)
    self.rinf_w2 = np.zeros((self.rinf_NC1,self.rinf_NC2), dtype=np.float64)
    self.rinf_w3 = np.zeros((self.rinf_NC2), dtype=np.float64)
    self.rinf_moy = np.zeros((self.rinf_NES), dtype=np.float64)
    self.rinf_ecart = np.zeros((self.rinf_NES), dtype=np.float64)
  
  # -- METHOD read_LUTs_Kd ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Perceptron LUTs reading
  # LUTS are expected to be located in the directory defined by environment variable $KD_LUTS_PATH (cf: __init__)
  # --
  def read_LUTs_Kd(self):
    print "class KdJamet.read_LUTs_Kd Method"
    # ----------------------------------------------------------------------------------
    #  < .85 weights settings
    #print "  KdJamet.read_LUTs_Kd : read weights file for Rrs%d/Rrs%d < .85: %s/%s" % (self.WL[2], self.WL[4],self.lutpath, self.rinf_LUT_POIDS)
    weights = np.loadtxt( "%s/%s" % (self.lutpath, self.rinf_LUT_POIDS), skiprows=1)
    base = 0
    self.rinf_b1 = weights[base:base+self.rinf_NC1,2]
    base = base+self.rinf_NC1
    self.rinf_b2 = weights[base:base+self.rinf_NC2,2]
    base = base+self.rinf_NC2
    self.rinf_b3 = weights[base:base+self.rinf_NS,2]
    base = base+self.rinf_NS
    for i in range(self.rinf_NE):
      self.rinf_w1[i,:] = weights[base:base+self.rinf_NC1,2]
      base = base+self.rinf_NC1
    for i in range(self.rinf_NC1):
      self.rinf_w2[i,:] = weights[base:base+self.rinf_NC2,2]
      base = base+self.rinf_NC2
    self.rinf_w3 = weights[base:base+self.rinf_NC2,2]
    
    usedData = np.array([1,2,3,4,5,6,8])
    if self.sensor == 'idMODIS':
      usedData = np.array([1,2,3,4,6,7,9])
    #print "  KdJamet.read_LUTs_Kd : read mean file : %s/%s" % (self.lutpath, self.rinf_LUT_MOY)
    moy = np.loadtxt( "%s/%s" % (self.lutpath, self.rinf_LUT_MOY))
    self.rinf_moy = moy[usedData]
    #print "  class KdJamet.read_LUTs_Kd : read ecart file : %s/%s" % (self.lutpath, self.rinf_LUT_ECART)
    ecart = np.loadtxt( "%s/%s" % (self.lutpath, self.rinf_LUT_ECART))
    self.rinf_ecart = ecart[usedData]
    
    # ----------------------------------------------------------------------------------
    #  >= .85 weights settings
    #print "  KdJamet.read_LUTs_Kd : read weights file for Rrs%d/Rrs%d >= .85: %s/%s" % (self.WL[2], self.WL[4], self.lutpath, self.rsup_LUT_POIDS)
    weights = np.loadtxt( "%s/%s" % (self.lutpath, self.rsup_LUT_POIDS), skiprows=1)
    base = 0
    self.rsup_b1 = weights[base:base+self.rsup_NC1,2]
    base = base+self.rsup_NC1
    self.rsup_b2 = weights[base:base+self.rsup_NC2,2]
    base = base+self.rsup_NC2
    self.rsup_b3 = weights[base:base+self.rsup_NS,2]
    base = base+self.rsup_NS
    for i in range(self.rsup_NE):
      self.rsup_w1[i,:] = weights[base:base+self.rsup_NC1,2]
      base = base+self.rsup_NC1
    for i in range(self.rsup_NC1):
      self.rsup_w2[i,:] = weights[base:base+self.rsup_NC2,2]
      base = base+self.rsup_NC2
    self.rsup_w3 = weights[base:base+self.rsup_NC2,2]
    
    usedData = np.array([1,2,3,4,6,8])
    if self.sensor == 'idMODIS':
      usedData = np.array([1,2,3,4,7,9])
    moy = np.loadtxt( "%s/%s" % (self.lutpath, self.rsup_LUT_MOY))
    self.rsup_moy = moy[usedData]
    ecart = np.loadtxt( "%s/%s" % (self.lutpath, self.rsup_LUT_ECART))
    self.rsup_ecart = ecart[usedData]
    
  # -- METHOD compute_allKd ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Kd computing at all Sensor WaveLength (vary for each sensor)
  # Inputs -
  # - indata : array[] [ Rrs443, Rrs490, Rrs510, Rrs555, Rrs665 ]
  # Output -
  # - return Kd[] at Sensor Wave Length
  # --
  def compute_allKd(self, inRrs):
    outKd  = np.zeros(NBWL,dtype=np.float32)
    if inRrs[2] / inRrs[4] < .85 :
      indata = np.zeros(self.rinf_NE,dtype=np.float32)
      indata[0:self.rinf_NE-1] = inRrs[1:self.rinf_NE]
      idLambda = self.rinf_NE-1
      FLAGinf = True
    else:
      indata = np.zeros(self.rsup_NE,dtype=np.float32)
      indata[0:self.rsup_NE-1] = inRrs[1:self.rsup_NE]
      idLambda = self.rsup_NE-1
      FLAGinf = False
    inc = 0
    # compute Kd for each Wave Length of sensor used as iput
    for wl in self.WL:
      indata[idLambda] = wl
      if FLAGinf == True:
        outKd[inc] = self.passe_avant_inf(indata)
      else:
        outKd[inc] = self.passe_avant_sup(indata)
      inc = inc+1
    return outKd
    
  # --  METHOD passe_avant_sup ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Kd computing  for Rrs490/Rrs55 >= .85 at ONE WaveLength
  # Inputs -
  # - indata : array[rinf_NE] [ Rrs443, Rrs490, Rrs510, Rrs55, Lambda ]
  # Output -
  # - return y (Kd) at indata Wave Length {the wave length can be any value between min and max Wave length
  #   used for the sensor  }
  # --
  def passe_avant_sup(self, indata):
    # Normalization
    x = ((2/3.) * (indata[:] - self.rsup_moy[0:self.rsup_NE])) / self.rsup_ecart[0:self.rsup_NE]

    a = np.zeros(self.rsup_NC1,dtype=np.float64)
    for i in range(self.rsup_NC1):
      for j in range(self.rsup_NE):
        a[i] = a[i] + x[j] * self.rsup_w1[j][i]
    a[:] = 1.715905*np.tanh((2./3.)*(a[:] + self.rsup_b1[:]))

    b = np.zeros(self.rsup_NC2,dtype=np.float64)
    for i in range(self.rsup_NC2):
      for j in range(self.rsup_NC1):
        b[i] = b[i] + a[j] * self.rsup_w2[j][i]
    b[:] = 1.715905*np.tanh((2./3.)*(b[:] + self.rsup_b2[:]))

    y = np.sum(b[:] * self.rsup_w3[:])
    # Denormalisation 
    y = 1.5*(y + self.rsup_b3)*self.rsup_ecart[self.rsup_NES-1] + self.rsup_moy[self.rsup_NES-1]
    y = 10.**y
    return(y)
  
  # --   METHOD passe_avant_inf ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Kd computing  for Rrs490/Rrs55 < .85 at ONE WaveLength
  # Inputs 
  # - indata : array[rinf_NE] [ Rrs443, Rrs490, Rrs510, Rrs55, Rrs665, Lambda ]
  # Output -
  # - return y (Kd) at indata Wave Length {the wave length can be any value between min and max Wave length
  #   used for the sensor  }
  # --
  def passe_avant_inf(self, indata):
    # Normalization
    x = ((2/3.) * (indata[:] - self.rinf_moy[0:self.rinf_NE])) / self.rinf_ecart[0:self.rinf_NE]
    a = np.zeros(self.rinf_NC1,dtype=np.float64)
    for i in range(self.rinf_NC1):
      for j in range(self.rinf_NE):
        a[i] = a[i] + x[j] * self.rinf_w1[j][i]
    a[:] = 1.715905*np.tanh((2./3.)*(a[:] + self.rinf_b1[:]))

    b = np.zeros(self.rinf_NC2,dtype=np.float64)
    for i in range(self.rinf_NC2):
      for j in range(self.rinf_NC1):
        b[i] = b[i] + a[j] * self.rinf_w2[j][i]
    b[:] = 1.715905*np.tanh((2./3.)*(b[:] + self.rinf_b2[:]))

    y = np.sum(b[:] * self.rinf_w3[:])
    # Denormalisation 
    y = 1.5*(y + self.rinf_b3)*self.rinf_ecart[self.rinf_NES-1] + self.rinf_moy[self.rinf_NES-1]
    y = 10.**y
    return(y)