KdJametPython/class_KdJamet.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy as np
import sys
import os

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

class KdJamet:
  #
  # 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_LUT_PATH not set, add it to your environment variables"
      sys.exit()
    
    if Sensor.has_key(idSensor):
      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" % (idSensor)
  
  # - METHODE init_variables ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Perceptron parameters setting and array initialization for each Sensor
  # --
  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 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"

    #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)
    
    # =======================================================================================================
    # MODIS
    # =======================================================================================================
    if self.sensor == 'idMODIS':
      self.WL = np.array([412, 443, 488, 531, 547, 667])
    # =======================================================================================================
    # SEAWIFS
    # =======================================================================================================
    if self.sensor == 'idSEAWIFS':
      self.WL = np.array([412, 443, 490, 510, 555, 670])
  
  # -- METHODE read_LUTs_Kd ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # perceptron weights LUTs reading
  # LUTS are expected to be located in environment variable : $KD_LUTS_PATH
  # --
  def read_LUTs_Kd(self):
    print "class KdJamet.read_LUTs_Kd Method"
    # ----------------------------------------------------------------------------------
    #  < .85 weights settings
    print "  class 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 "  class 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 "  class 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]
    
  # -- METHODE compute_allKd ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # Kd computing at all Sensor used WaveLength
  # Inputs -
  # - indata : array[] [ Rrs443, Rrs490, Rrs510, Rrs55, 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
    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
    
  # --  METHODE 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
  # --
  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)
  
  # --   METHODE 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
  # --
  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)