TD lecture de fichier NKE

# Historique
* [X] Mise en place sous gitlab.osupytheas.fr
* [X] Ajout d'un debut de script TD7 (a remplir) + ajout des fichiers 2 à 10.txt (différentes series temporelles)
* [X] Ajout des notebooks TD2 à TD8 + ajout des fichiers
* [X] Ajout de requirements.txt pour gestion des dépendances
* [X] Ajout du fichier .py et du notebook jupyter .ipynb pour TD1
* [X] Création d'un git pour Stat

README.md

-Ajout d'un debut de script TD7 (a remplir) + ajout des fichiers 2 à 10.txt (différentes series temporelles)
+TD lecture de fichier NKE
 # Historique
+* [X] Mise en place sous gitlab.osupytheas.fr
+* [X] Ajout d'un debut de script TD7 (a remplir) + ajout des fichiers 2 à 10.txt (différentes series temporelles)
 * [X] Ajout des notebooks TD2 à TD8 + ajout des fichiers
 * [X] Ajout de requirements.txt pour gestion des dépendances
 * [X] Ajout du fichier .py et du notebook jupyter .ipynb pour TD1

requirements.txt

 cycler==0.10.0
 kiwisolver==1.1.0
-matplotlib==3.1.3
-numpy==1.18.1
-pandas==1.0.1
+matplotlib==3.3.4
+numpy==1.19.5
+pandas==1.1.5
+Pillow==8.1.0
 pyparsing==2.4.6
 python-dateutil==2.8.1
 pytz==2019.3

sources_python/script_TD_NKE_pression.py

+#!/usr/bin/python
+# coding: utf-8
+
+#################################################"    
+##                 MODULES                  #####"
+#################################################"
+
+import matplotlib.pyplot as plt
+from datetime import timedelta
+from datetime import datetime
+import matplotlib.dates as dt
+import numpy.matlib as matlib
+import scipy.io as scio
+from scipy import stats
+import pandas as pd
+import numpy as np
+import pylab as py
+import csv as csv
+import os as os
+import warnings
+import pickle
+import sys
+###############################################
+warnings.filterwarnings("ignore")
+
+
+
+                                                                      #      
+                                                                     ###                  
+                                                                    # # #
+                                                                      #
+                                                                      #
+##########################################################################################################################################################  
+                                                                      #
+                                                                      #
+                                                                    # # #
+                                                                     ###         
+                                                                      #
+###############################
+## CLASS TDNKE:
+##          
+###################################################################################   
+class TDNKE():
+###################################################################################   
+###############################    
+    def _load(self):
+        """ Load a ascii file from SP2T NKE"""
+        # Change comma by dot
+        [adr, ext]=os.path.splitext(os.path.basename(self.filename))
+        dateparse = lambda X: datetime.strptime(X,'%d/%m/%Y\t%H:%M:%S:%f')
+        self.mydata = pd.read_table(self.filename,header=4,
+                            decimal='.',
+                            usecols=[0,1,2, 3],
+                            names=["Temperature", "level","thedate","thetime"],
+                            parse_dates={"Date": ["thedate","thetime"]},
+                            date_parser=dateparse,
+                            delim_whitespace = True,
+                            #delimiter='\s+',
+                            #skiprows=100202,
+                            #nrows=2537035,
+                            skip_blank_lines=True)
+        self.mydata.set_index('Date', inplace=True)
+        self.Level_unit='m'
+        print ( '-----------------------------------------------')
+        print ('Lecture ', adr, '............OK')
+        print ('-----------------------------------------------')
+
+###################################################################################   
+    def _savepickle(self):
+        """ Save everything in a pickle"""
+        filenamepickle=self.filename[:-3]+'pkl'
+        DataPickle = open(filenamepickle, "wb")
+        pickle.dump(self, DataPickle)
+        DataPickle.close()
+        message='{:<10}{:.<30}{:.>40}'.format('Save:', filenamepickle, 'Ok')
+        print ('---------------------------------------------------------------')
+        print (message)
+        print ('---------------------------------------------------------------')
+
+################################################
+    def _spectralanalysis(self, Tf='No'):
+            # Pour Pandas
+            ### Resize the data into matrix
+            meandate,X=self._setdata(self.mydata.index,self.mydata.level.values,self.burst)            
+            a,N=np.shape(X)
+            # Mean water Depth
+            prof=np.nanmean(X,axis=1)
+            prof2=np.rot90(np.tile(prof, (N,1)))    
+            X2=X-prof2
+            #remove the linear trend
+            X2=self._detrend(X2)  
+
+            #  Spectral analysis
+            Pxx=np.abs(np.fft.rfft(X2,N)/(N/2))
+            a,N2=np.shape(Pxx)
+            freqs=np.linspace(0,self.Sample_Frequency/2,N2)
+            df=freqs[1]-freqs[0]
+
+            if Tf=='yes':
+                # Correction with transfer function (linear theory)
+                Tff=func_lineartheory(Pxx, freqs, prof, altmes)
+                Pxx =Pxx/Tff
+            
+            Pxx2=0.5*Pxx**2/df        
+            M=self._waveheight(Pxx2, freqs) #Hs
+            FFT=pd.DataFrame(M,index=['Hswind', 'Hs2','Hsinfra','HsVLF','Tp','Tm_wind', 'Tm_swell', 'Tm_infra','Tm_VLF'])
+            FFT=FFT.T
+            FFT['Date']=dt.num2date(meandate)
+            self.FFT=FFT.set_index('Date')
+            self.spectre=Pxx2
+            message='{:.<40}{:.>20}'.format('Process: FFT spectral analysis', 'Ok')   
+            print (message)
+
+#################################################################################################################
+    def _setdata(self,t, x, burst, nbr=None):
+        #Pour Pandas
+        ''' Set the x data vector into a X data matrix of burst x nbrburst size '''
+        
+        measurements=len(x)# Nbr of data
+        if nbr==None:
+            nbr=int(np.ceil(float(measurements)/burst)*burst)   # Nbr of data for a vector of burst x nbrburst    
+        x_reshape=np.empty([nbr])
+        x_reshape.fill(np.nan)
+        
+        x_reshape[0:measurements, ]=x # Fill the vector of x data
+        X=x_reshape.reshape([int(nbr/burst),burst]) # Reshape into matrix
+        # Calculate first time of each burst
+        t_reshape=np.empty([nbr])
+        t_reshape.fill(np.nan)
+        t=dt.date2num(t.values)
+        t_reshape[0:measurements]=t
+        t_burst=t_reshape.reshape([int(nbr/burst),burst])
+        return t_burst[:,0],X
+###################################################################################   
+    def _detrend(self,Y, ax=1):
+        ''' from a matrix Y, detrend the data along the axis ax'''
+        # find linear regression line, subtract off data to detrend
+        
+        [a,b]=Y.shape
+
+        
+        if ax==0:
+            x=np.array(range(0,a))
+            DY=[]
+            for i in range(b):
+                Yl=Y[:,i]
+                not_nan_ind = ~np.isnan(Yl)
+                m, b, r_val, p_val, std_err = stats.linregress(x[not_nan_ind],Yl[not_nan_ind])
+                detrend_y = Yl - (m*x + b)
+                DY.append(detrend_y)
+                
+        if ax==1:
+            x=np.array(range(0,b))
+            DY=[]
+            for i in range(a):
+                Yl=Y[i,:]
+                not_nan_ind = ~np.isnan(Yl)
+                m, b, r_val, p_val, std_err = stats.linregress(x[not_nan_ind],Yl[not_nan_ind])
+                detrend_y = Yl - (m*x + b)
+                DY.append(detrend_y)
+
+           
+        DY=np.array(DY)
+        message='{:.<40}{:.>20}'.format('Process: detrend', 'Ok')
+        print(message)
+        return DY
+
+################################################
+    def _waveheight(self,S, f, cof=0.5, wind= 0.3, infra=0.05, lowinfra=0.004, VLF=0.001):
+
+        """" From a spectrum, Calculate differents wave height"""
+        
+        # parameters
+        df=f[1]-f[0]
+        p=[0,1, 2, -1, -2]
+
+        #integration
+        S_swell=np.copy(S)
+        S_infra=np.copy(S)
+        S_wind=np.copy(S)
+        S_VLF=np.copy(S)    
+        S_wind[:,np.where(np.logical_or(f>cof,f<wind))]=0
+        S_swell[:,np.where(np.logical_or(f>wind,f<infra))]=0
+        S_infra[:,np.where(np.logical_or(f>infra,f<lowinfra))]=0
+        S_VLF[:,np.where(np.logical_or(f>lowinfra,f<VLF))]=0
+        # moment calculation
+        m0_wind=np.nansum(S_wind, axis=1)*df
+        m0_swell=np.nansum(S_swell, axis=1)*df
+        m0_infra=np.nansum(S_infra, axis=1)*df
+        m0_VLF=np.nansum(S_VLF, axis=1)*df
+
+        m1_wind=np.nansum(f**p[1]*S_wind*df, axis=1)
+        m2_wind=np.nansum(f**p[2]*S_wind*df, axis=1)
+        m_1_wind=np.nansum(f**p[3]*S_wind*df, axis=1)
+        m_2_wind=np.nansum(f**p[4]*S_wind*df, axis=1)
+
+        
+        m1_swell=np.nansum(f**p[1]*S_swell*df, axis=1)
+        m2_swell=np.nansum(f**p[2]*S_swell*df, axis=1)
+        m_1_swell=np.nansum(f**p[3]*S_swell*df, axis=1)
+        m_2_swell=np.nansum(f**p[4]*S_swell*df, axis=1)
+        
+        m1_infra=np.nansum(f**p[1]*S_infra*df, axis=1)
+        m2_infra=np.nansum(f**p[2]*S_infra*df, axis=1)
+        m_1_infra=np.nansum(f**p[3]*S_infra*df, axis=1)
+        m_2_infra=np.nansum(f**p[4]*S_infra*df, axis=1)
+
+        m1_VLF=np.nansum(f**p[1]*S_VLF*df, axis=1)
+        m2_VLF=np.nansum(f**p[2]*S_VLF*df, axis=1)
+        m_1_VLF=np.nansum(f**p[3]*S_VLF*df, axis=1)
+        m_2_VLF=np.nansum(f**p[4]*S_VLF*df, axis=1)
+
+        #Height
+        Hs_wind=4*np.sqrt(m0_wind)       
+        Hs_swell=4*np.sqrt(m0_swell)   
+        Hs_infra=4*np.sqrt(m0_infra)
+        Hs_VLF=4*np.sqrt(m0_VLF)
+        # Period
+        position=np.nanargmax(S_swell, axis=1)
+        fmax=f[position]
+        Tp=1/fmax
+
+        Tm1_wind=(m1_wind/m0_wind)**-(1/p[1])
+        Tm2_wind=(m2_wind/m0_wind)**-(1/p[2])
+        Tm_1_wind=(m_1_wind/m0_wind)**-(1/p[3])
+        Tm_2_wind=(m_2_wind/m0_wind)**-(1/p[4])
+        Tm_wind=[Tm1_wind, Tm2_wind, Tm_1_wind, Tm_2_wind]
+        
+        Tm1_swell=(m1_swell/m0_swell)**-(1/p[1])
+        Tm2_swell=(m2_swell/m0_swell)**-(1/p[2])
+        Tm_1_swell=(m_1_swell/m0_swell)**-(1/p[3])
+        Tm_2_swell=(m_2_swell/m0_swell)**-(1/p[4])
+        Tm_swell=[Tm1_swell, Tm2_swell, Tm_1_swell, Tm_2_swell]
+
+        Tm1_infra=(m1_infra/m0_infra)**-(1/p[1])
+        Tm2_infra=(m2_infra/m0_infra)**-(1/p[2])
+        Tm_1_infra=(m_1_infra/m0_infra)**-(1/p[3])
+        Tm_2_infra=(m_2_infra/m0_infra)**-(1/p[4])
+        Tm_infra=[Tm1_infra, Tm2_infra, Tm_1_infra, Tm_2_infra]
+
+        
+        Tm1_VLF=(m1_VLF/m0_VLF)**-(1/p[1])
+        Tm2_VLF=(m2_VLF/m0_VLF)**-(1/p[2])
+        Tm_1_VLF=(m_1_VLF/m0_VLF)**-(1/p[3])
+        Tm_2_VLF=(m_2_VLF/m0_VLF)**-(1/p[4])
+        Tm_VLF=[Tm1_VLF, Tm2_VLF, Tm_1_VLF, Tm_2_VLF]
+            
+        message='{:.<40}{:.>20}'.format('Process: Wave height calculation', 'Ok')
+        print (message)
+        return Hs_wind, Hs_swell, Hs_infra, Hs_VLF, 1/fmax, Tm_wind, Tm_swell, Tm_infra, Tm_VLF
+
+        
+        
+#################################################################################################################
+    def fig2D(self,t,y,ly='labely',figname='figure.png', datef='%d %Hh'):
+        """ Create a 2D plot = time x variable (Hs, Pression, depth...)"""
+       
+        # Figure creation        
+        taille = (10,5) # 10*tdpi x 4*tdpi -> 1200x800 si tdpi =200
+        tdpi=200
+        fig1=plt.figure(figsize=taille)
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')    
+        lx=r''+'Time'
+        ly=r''+ly
+        plt.plot(t,y, 'r.')
+        plt.xlabel(lx, fontsize=16)
+        plt.ylabel(ly, fontsize=16)
+        plt.grid()
+        date_format = dt.DateFormatter(datef)
+        plt.autoscale(enable=True, axis='both', tight=True)
+        plt.gca().xaxis.set_major_formatter(date_format)
+      
+        #--------------------------------------------------
+        #   Save the figure
+        #--------------------------------------------------
+        fig1.savefig(figname,dpi=tdpi)
+        message='Figure 2D:{:<55},{:.>} dpi{:.>10}'.format(figname , tdpi,'Ok')
+        print (message)   
+
+
+                                                                      #      
+                                                                     ###                  
+                                                                    # # #
+                                                                      #
+                                                                      #
+##########################################################################################################################################################  
+                                                                      #
+                                                                      #
+                                                                    # # #
+                                                                     ###         
+                                                                      #
+##########################################################"    
+##                 PYTHON SCRIPT                     #####"
+##########################################################"
+
+
+#############################
+### DEF VARIABLES
+############################# 
+deg = u'\xb0'
+data_NKE=TDNKE()
+# Choose the working directory
+os.chdir('../../../workspace/marseille')
+data_NKE.filename='Sonde standard_30098_20190309_102905.txt'
+data_NKE.Sample_Frequency=4
+data_NKE.burst=2400 # Nombre de point pour l'analyse spectrale
+##########################################################  
+###              OPEN FILE                             ###
+##########################################################  
+# READ .TXT FILE    
+data_NKE._load() # Pressure in m
+########################################################## 
+###          SPECTRAL ANALYSIS                       #####
+##########################################################
+data_NKE._spectralanalysis()
+########################################################## 
+###                  OUTPUT                          #####
+##########################################################
+# FIGURES
+
+figname=data_NKE.filename[:-4] + '_Pressure.png'
+data_NKE.fig2D(data_NKE.mydata.index, data_NKE.mydata.level,'Pressure (bar)', figname, '%y/%m/%d')
+
+figname=data_NKE.filename[:-4] + '_Temperature.png'
+data_NKE.fig2D(data_NKE.mydata.index,data_NKE.mydata.Temperature,'Temperature (degree C)', figname, '%y/%m/%d ')
+
+figname=data_NKE.filename[:-4] + '_Hs2_fft.png'
+data_NKE.fig2D(data_NKE.FFT.index,data_NKE.FFT.Hs2,"Significant wave height (m)",  figname, '%y/%m/%d')
+
+figname=data_NKE.filename[:-4]+ '_Tp.png'
+data_NKE.fig2D(data_NKE.FFT.index,data_NKE.FFT.Tp,"Peak period(s) fft method",  figname, '%y/%m/%d')
+
+data_NKE._savepickle()
+
+                                                                      #      
+                                                                     ###                  
+                                                                    # # #
+                                                                      #
+                                                                      #
+##########################################################################################################################################################
+#import pdb; pdb.set_trace()