123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335 |
- # -*- 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)
|