Quentin 26/03/2021: add a new python3.7 library to plot KTEST

src/LIB/Python/misc_functions.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 25 11:25:31 2021
+
+@author: rodierq
+"""
+import copy
+from scipy.interpolate import RectBivariateSpline
+import numpy as np
+import math
+
+def convert_date(datesince, time_in_sec):
+  print(type(datesince))
+  print(type(time_in_sec))
+  return str(time_in_sec) + datesince[:33]
+
+class mean_operator():
+    def MYM(self,var):
+        ny = var.shape[1]
+        out = copy.deepcopy(var)
+        for j in range(ny-1):
+            out[:,j,:] = (var[:,j,:] + var[:,j+1,:])*0.5
+        return out
+    
+    def MXM(self,var):
+        nx = var.shape[2]
+        out = copy.deepcopy(var)
+        for i in range(nx-1):
+            out[:,:,i] = (var[:,:,i] + var[:,:,i+1])*0.5
+        return out
+    
+    def MZM(self,var):
+        nz = var.shape[0]
+        out = copy.deepcopy(var)
+        for k in range(nz-1):
+            out[k,:,:] = (var[k,:,:] + var[k+1,:,:])*0.5
+        return out
+
+def windvec_verti_proj(u, v, lvl, angle):
+  """
+     Compute the projected horizontal wind vector on an axis with a given angle w.r.t. the x/ni axes (West-East)
+     Arguments :
+            - u : U-wind component
+            - v : V-wind component
+            - angle (radian) of the new axe w.r.t the x/ni axes (West-East). angle = 0 for (z,x) sections, angle=pi/2 for (z,y) sections
+     Returns :
+            - a 3D wind component projected on the axe to be used with Panel_Plot.vector as Lvar1
+  """ 
+  out = copy.deepcopy(u)
+  for k in range(len(lvl)):
+      out[k,:,:] = u[k,:,:]*math.cos(angle) + v[k,:,:]*math.sin(angle)
+  return out
+
+def oblique_proj(var, ni, nj, lvl, i_beg, j_beg, i_end, j_end):
+  """
+     Compute an oblique projection of a 3D variable w.r.t. its axes
+     Arguments :
+            - var          : 3D var to project (example THT)
+            - ni           : 1D x-axe of the 3D dimension
+            - nj           : 1D y-axe of the 3D dimension
+            - level        : 1D level variable (level or level_w)
+            - i_beg, j_beg : coordinate of the begin point of the new axe
+            - i_end, j_end : coordinate of the end point of the new axe
+     Returns :
+            - a 2D (z,m) variable projected on the oblique axe
+            - a 1D m new axe (distance from the beggining point)
+            - the angle (radian) of the new axe w.r.t the x/ni axes (West-East)
+  """
+  dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) #  Distance de la section oblique  m
+  out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
+  axe_m = np.zeros(int(dist_seg)+1) #Axe des abscisses qui sera tracé selon la coupe
+  axe_m_coord = [] #Coordonnées x,y des points qui composent l'axe
+  axe_m_coord.append( (ni[i_beg],nj[j_beg]) ) #Le premier point est celui donné par l'utilisateur
+  for m in range(int(dist_seg)): #Discrétisation selon distance de la coupe / int(distance_de_la_coupe)
+    axe_m_coord.append( (axe_m_coord[0][0] + (ni[i_end]-ni[i_beg])/(int(dist_seg))*(m+1), 
+                         axe_m_coord[0][1] + (nj[j_end]-nj[j_beg])/(int(dist_seg))*(m+1) ))
+    axe_m[m+1] = np.sqrt((ni[i_beg]-axe_m_coord[m+1][0])**2 + (nj[j_beg]-axe_m_coord[m+1][1])**2)
+  
+  for k in range(len(lvl)):
+    a=RectBivariateSpline(ni, nj,var[k,:,:],kx=1,ky=1) #Interpolation par niveau à l'ordre 1 pour éviter des valeurs négatives de champs strictement > 0
+    for m in range(int(dist_seg)+1):
+      out_var[k,m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1]) # La fonction ev de RectBivariate retourne la valeur la plus proche du point considéré
+      
+  angle_proj = math.acos((ni[i_end]-ni[i_beg])/axe_m[-1])
+  return angle_proj, out_var, axe_m
+
+def comp_altitude1DVar(oneVar2D, orography, ztop, level, n_xory):
+  """
+     Compute and returns an altitude and x or y grid mesh variable in 2D following the topography in 1D
+     To be used with 2D simulations
+     Arguments :
+            - oneVar2D  : a random netCDF 2D var (example UT, THT)
+            - orography : 1D orography (ZS)
+            - ztop      : scalar of the top height of the model (ZTOP)
+            - level     : 1D level variable (level or level_w)
+            - n_xory:   : 1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
+     Returns :
+            - a 2D altitude variable with topography taken into account
+            - a 2D directionnal variable
+  """
+  n_xory_2D = copy.deepcopy(oneVar2D)
+  altitude  =  copy.deepcopy(oneVar2D)
+  for i in range(len(level)):
+    n_xory_2D[i,:] = n_xory
+  for j in range(len(n_xory)):
+    for k in range(len(level)):
+        altitude[k,j] = orography[j] + level[k]*((ztop-orography[j])/ztop)
+  return altitude, n_xory_2D
+
+def comp_altitude2DVar(oneVar3D, orography, ztop, level, n_y, n_x):
+  """
+     Compute and returns an altitude and x or y grid mesh variable in 3D following the topography in 2D
+     To be used with 3D simulations in cartesian coordinates
+     Arguments :
+            - oneVar3D  : a random netCDF 3D var (example UT, THT)
+            - orography : 2D orography (ZS)
+            - ztop      : scalar of the top height of the model (ZTOP)
+            - level     : 1D level variable (level or level_w)
+            - n_xory:   : 1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
+     Returns :
+            - a 3D altitude variable with topography taken into account
+            - a 3D directionnal variable
+  """
+  n_x3D = copy.deepcopy(oneVar3D)
+  n_y3D = copy.deepcopy(oneVar3D)
+  altitude  =  copy.deepcopy(oneVar3D)
+  for i in range(len(level)):
+    n_y3D[i,:] = n_y
+    n_x3D[i,:] = n_x
+  for i in range(oneVar3D.shape[2]):
+    for j in range(oneVar3D.shape[1]):
+      for k in range(len(level)):
+          altitude[k,j,i] = orography[j,i] + level[k]*((ztop-orography[j,i])/ztop)
+  return altitude, n_x3D, n_y3D

src/LIB/Python/read_MNHfile.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Feb 22 10:29:13 2021
+
+@author: rodierq
+"""
+
+import netCDF4 as nc
+import numpy as np
+
+def read_netcdf(LnameFiles, Dvar_input, path='.', removeHALO=True):
+  Dvar = {}
+  for i,nameFiles in enumerate(LnameFiles):
+    f_nb = 'f' + str(i+1)
+    print('Reading file ' + f_nb)
+    print(path + nameFiles)
+    if '000' in nameFiles[-6:-3]: #time series file (diachronic)
+        theFile = nc.Dataset(path + nameFiles,'r')
+        if theFile['MASDEV'][0] <= 54:
+            read_TIMESfiles_54(theFile, f_nb, Dvar_input, Dvar)
+        else : # 55 groups variables
+            read_TIMESfiles_55(theFile, f_nb, Dvar_input, Dvar, removeHALO)
+        theFile.close()
+    else:
+        read_BACKUPfile(nameFiles, f_nb, Dvar_input, Dvar, path=path, removeHALO=removeHALO)
+  return Dvar #Return the dic of [files][variables]
+
+def read_BACKUPfile(nameF, ifile, Dvar_input, Dvar_output, path='.', removeHALO=True):
+  theFile = nc.Dataset(path + nameF,'r')
+  Dvar_output[ifile] = {} #initialize dic for each files 
+       
+  #  Reading date since beginning of the model run
+  Dvar_output[ifile]['date'] =  theFile.variables['time'].units
+  Dvar_output[ifile]['time'] =  theFile.variables['time'][0]
+
+  for var in Dvar_input[ifile]: #For each files
+    #  Read variables
+    n_dim = theFile.variables[var].ndim
+    name_dim = theFile.variables[var].dimensions
+ #   print(var + "  n_dim  = " + str(n_dim))
+    
+    if (n_dim ==0) or (n_dim == 1 and 'time' in name_dim):  #Scalaires or Variable time
+      Dvar_output[ifile][var] = theFile.variables[var][0].data
+    else:      
+        if removeHALO:
+          if n_dim == 1:
+            Dvar_output[ifile][var] = theFile.variables[var][1:-1] #Variables 1D
+          elif n_dim == 2:
+            if theFile.variables[var].shape[0] == 1 and 'size' in name_dim[1]: #Variables 2D with the second dimension not a coordinate (--> list of strings : chemicals)
+                Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
+            elif theFile.variables[var].shape[0] == 1: #Variables 2D with first dim = 0
+                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1] #Variables 2D
+            else:
+              Dvar_output[ifile][var] = theFile.variables[var][1:-1,1:-1] #Variables 2D
+          elif n_dim == 5: #Variables time, sizeXX, level, nj, ni (ex: chemical budgets in 3D)
+              Dvar_output[ifile][var] = theFile.variables[var][0, :, 1:-1,:1:-1,1:-1]
+          elif n_dim == 4 and 'time' in name_dim and ('level' in name_dim or 'level_w' in name_dim): # time,z,y,x
+              if theFile.variables[var].shape[1] == 1: #Variables 4D time,z,y,x with time=0 z=0
+                Dvar_output[ifile][var] = theFile.variables[var][0,0,1:-1,1:-1] #Variables 2D y/x
+              elif theFile.variables[var].shape[2] == 1: #Variable 2D  (0,zz,0,xx)
+                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,0,1:-1] #Variables 2D z/y
+              elif theFile.variables[var].shape[3] == 1: #Variable 2D  (0,zz,yy,0)
+                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1,0] #Variables 2D z/x
+              ## ATTENTION VARIABLE 1D codé en 4D non faite
+              else: #Variable 3D simple
+                Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1,1:-1] #Variables time + 3D     
+          elif n_dim == 4 and 'time' in name_dim and 'level' not in name_dim and 'level_w' not in name_dim: # time,nb_something,y,x
+               Dvar_output[ifile][var] = theFile.variables[var][0,:,1:-1,1:-1] #Variables 2D y/x
+          elif n_dim == 3 and 'time' in name_dim: # time, y, x 
+            Dvar_output[ifile][var] = theFile.variables[var][0,1:-1,1:-1]
+          else :
+            Dvar_output[ifile][var] = theFile.variables[var][1:-1,1:-1,1:-1]  #Variables 3D
+        else:
+          if n_dim == 1:
+            Dvar_output[ifile][var] = theFile.variables[var][:] #Variables 1D  
+          elif n_dim == 2:
+            if theFile.variables[var].shape[0] == 1 and 'size' in name_dim[1]: #Variables 2D with the second dimension not a coordinate (--> list of strings : chemicals)
+              Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
+            elif theFile.variables[var].shape[0] == 1: #Variables 2D with first dim = 0
+              Dvar_output[ifile][var] = theFile.variables[var][0,:] #Variables 2D
+            else:
+              Dvar_output[ifile][var] = theFile.variables[var][:,:] #Variables 2D
+          elif n_dim == 5: #Variables time, sizeXX, level, nj, ni (ex: chemical budgets in 3D)
+              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,:,:]
+          elif n_dim == 4: # time,z,y,x
+            if theFile.variables[var].shape[1] == 1: #Variables 4D time,z,y,x with time=0 z=0
+              Dvar_output[ifile][var] = theFile.variables[var][0,0,:,:] #Variables 2D y/x
+            elif theFile.variables[var].shape[2] == 1: #Variable 2D  (0,zz,0,xx)
+              Dvar_output[ifile][var] = theFile.variables[var][0,:,0,:] #Variables 2D z/y
+            elif theFile.variables[var].shape[3] == 1: #Variable 2D  (0,zz,yy,0)
+              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,0] #Variables 2D z/x
+            ## ATTENTION VARIABLE 1D codé en 4D non faite
+            else: #Variable 3D simple
+              Dvar_output[ifile][var] = theFile.variables[var][0,:,:,:] #Variables time + 3D
+          elif n_dim ==3 and name_dim in var.dimensions: # time, y, x
+            Dvar_output[ifile][var] = theFile.variables[var][0,:,:]
+          else:
+            Dvar_output[ifile][var] = theFile.variables[var][:,:,:]  #Variables 3D
+        #  For all variables except scalars, change Fill_Value to NaN
+        Dvar_output[ifile][var]= np.where(Dvar_output[ifile][var] != -99999.0, Dvar_output[ifile][var], np.nan)
+        Dvar_output[ifile][var]= np.where(Dvar_output[ifile][var] != 999.0, Dvar_output[ifile][var], np.nan)
+
+  theFile.close()
+  return Dvar_output #Return the dic of [files][variables]
+
+def read_TIMESfiles_54(theFile, ifile, Dvar_input, Dvar_output):
+    Dvar_output[ifile] = {} #initialize dic for each files 
+
+    #  Level variable is automatically read without the Halo
+    Dvar_output[ifile]['level'] = theFile.variables['level'][1:-1]
+    
+    #  Time variable is automatically read (time since begging of the run) from the 1st variable of the asked variable to read
+    suffix, name_first_var = remove_PROC(Dvar_input[ifile][0])
+    try: #  It is possible that there is only one value (one time) in the .000 file, as such time series are not suitable and the following line can't be executed. The time variable is then not written
+        increment = theFile.variables[name_first_var+'___DATIM'][1,-1] - theFile.variables[name_first_var+'___DATIM'][0,-1] #-1 is the last entry of the date (current UTC time in seconds)
+        length_time = theFile.variables[name_first_var+'___DATIM'].shape[0]
+        Dvar_output[ifile]['time'] = np.arange(increment,increment*(length_time+1),increment)
+    except:
+        pass
+    
+    for var in Dvar_input[ifile]: #For each files
+        suffix, var_name = remove_PROC(var)
+        n_dim = theFile.variables[var].ndim
+        name_dim = theFile.variables[var].dimensions
+
+        #  First, test if the variable is a dimension/coordinate variable
+        if (n_dim ==0):  #  Scalaires variable
+             Dvar_output[ifile][var] = theFile.variables[var][0].data
+             pass
+        elif n_dim == 1:
+            Dvar_output[ifile][var_name] = theFile.variables[var][1:-1]  #  By default, the Halo is always removed because is not in the other variables in any .000 variable
+            pass
+        elif n_dim == 2:
+            Lsize1 = list_size1(n_dim, name_dim)
+            if Lsize1 == [True, False]: Dvar_output[ifile][var_name] = theFile.variables[var][0,1:-1] 
+            pass
+        
+        Lsize1 = list_size1(n_dim, name_dim)
+        if Lsize1 == [True, False, False, True, True]: Dvar_output[ifile][var_name] = theFile.variables[var][0,:,:,0,0].T # Need to be transposed here
+        if Lsize1 == [True, True, False, True, False]: Dvar_output[ifile][var_name] = theFile.variables[var][0,0,:,0,:]
+
+    return Dvar_output #Return the dic of [files][variables]
+
+def read_TIMESfiles_55(theFile, ifile, Dvar_input, Dvar_output, removeHALO=False):
+    Dvar_output[ifile] = {} #initialize dic for each files 
+      
+    for var in Dvar_input[ifile]: #For each var
+        suffix, var_name = remove_PROC(var)
+
+        try: #  NetCDF4 Variables
+            n_dim = theFile.variables[var_name].ndim
+            #  First, test if the variable is a dimension/coordinate variable
+            if (n_dim ==0):  #  Scalaires variable
+                Dvar_output[ifile][var_name] = theFile.variables[var_name][0].data
+            else:
+                if(removeHALO):
+                    if n_dim == 1:
+                        Dvar_output[ifile][var_name] = theFile.variables[var_name][1:-1]
+                    elif n_dim == 2:
+                        Dvar_output[ifile][var_name] = theFile.variables[var_name][1:-1,1:-1] 
+                    else: 
+                        raise NameError("Lecture des variables de dimension sup a 2 pas encore implementees pour fichiers .000")
+                else:
+                    if n_dim == 1:
+                        Dvar_output[ifile][var_name] = theFile.variables[var_name][:]
+                    elif n_dim == 2:
+                        Dvar_output[ifile][var_name] = theFile.variables[var_name][:,:] 
+                    else: 
+                        raise NameError("Lecture des variables de dimension sup a 2 pas encore implementees pour fichiers .000")
+        except KeyError: # NetCDF4 Group
+            if theFile.groups[var_name].type == 'TLES' : #LES type
+                #  Build the true name of the groups.variables :
+                whites = ' '*(17 - len('(cart)') - len(var_name)) #  TODO : a adapter selon le type de la variable cart ou autres
+                Dvar_output[ifile][var] = theFile.groups[var].variables[var + whites + '(cart)'][:,:].T
+            elif theFile.groups[var_name].type == 'CART':  #  Budget CART type
+                shapeVar = theFile.groups[var_name].variables[suffix].shape
+                Ltosqueeze=[] #  Build a tuple with the number of the axis which are 0 dimensions to be removed by np.squeeze
+                if shapeVar[0]==1: Ltosqueeze.append(0)
+                if shapeVar[1]==1: Ltosqueeze.append(1)
+                if shapeVar[2]==1: Ltosqueeze.append(2)
+                if shapeVar[3]==1: Ltosqueeze.append(3)
+                Ltosqueeze=tuple(Ltosqueeze)
+                Dvar_output[ifile][var_name] = np.squeeze(theFile.groups[var_name].variables[suffix][:,:,:,:], axis=Ltosqueeze) 
+            else:
+                raise NameError("Type de groups variables not implemented in read_MNHfile.py")
+    return Dvar_output #Return the dic of [files][variables]
+
+def list_size1(n_dim, named_dim):
+    Lsize1 = []
+    for i in range(n_dim):
+        if 'size1' == named_dim[i]:
+            Lsize1.append(True)
+        else:
+            Lsize1.append(False)
+    return Lsize1
+    
+def remove_PROC(var):
+    if '___PROC' in var:
+        var_name = var[:-8]
+        suffix = "" # No need of suffix for MNHVERSION < 550 (suffix is for NetCDF4 group
+    elif  '___ENDF' in var or '___INIF' in var or '___AVEF' in var:
+        var_name = var[:-7]
+        suffix = var[-4:]
+    else:
+        var_name = var
+        suffix = ''
+    return suffix, var_name