diff --git a/src/LIB/Python/Panel_Plot.py b/src/LIB/Python/Panel_Plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..38952c114964e03fc0beb93f0435cc8823947117
--- /dev/null
+++ b/src/LIB/Python/Panel_Plot.py
@@ -0,0 +1,791 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 24 10:49:45 2021
+
+@author: rodierq
+"""
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from matplotlib.colors import ListedColormap
+import numpy as np
+import cartopy.feature as cfeature
+
+class PanelPlot():
+
+ def __init__(self, nb_l, nb_c, Lfigsize, bigtitle, titlepad=40, minmaxpad=1.03, timepad=-0.06, lateralminmaxpad=0.86,
+ labelcolorbarpad=6.0, colorbaraspect=20, colorbarpad=0.04 ):
+ self.bigtitle = bigtitle
+ self.Lfigsize = Lfigsize
+ self.nb_l = nb_l
+ self.nb_c = nb_c
+ self.nb_graph = 0 # New independent graph within the subplot
+ self.titlepad = titlepad # Title pad (vertical shift) from graph
+ self.minmaxpad = minmaxpad # Min/Max print pad (vertical shift)
+ self.timepad = timepad # Time print pad (vertical shift)
+ self.colorbarpad = colorbarpad # Colorbar pad (horizontal shift from graph)
+ self.lateralminmaxpad = lateralminmaxpad
+ self.labelcolorbarpad = labelcolorbarpad # Vertical colorbal label pad
+ self.colorbaraspect = colorbaraspect # Ratio of long to short dimensions of colorbar w.r.t. the figure
+
+ # Initialization of the panel plots
+ self.fig = plt.figure(figsize=(self.Lfigsize[0],self.Lfigsize[1]))
+ self.fig.set_dpi(400)
+ self.fig.suptitle(self.bigtitle,fontsize=16)
+
+ def save_graph(self, iplt, fig):
+ """
+ Create a temporary png file of (sub)-plot(s) which can be converted to PDF
+ """
+ self.iplt = iplt
+ self.fig = fig
+
+ self.iplt+=1
+ self.fig.savefig('tempgraph'+str(self.iplt))
+ return self.iplt
+
+ def draw_Backmap(self,drawCoastLines, ax, projo):
+ """
+ Handle drawing of the background plot (coastlines, departements, grid lines and labels)
+ """
+ self.drawCoastLines = drawCoastLines
+ self.projo = projo
+
+ # Grid lines and labels
+ if 'PlateCarree' in str(projo):
+ gl = ax.gridlines(crs=self.projo, draw_labels=True, linewidth=1, color='gray')
+ gl.xlabels_top = False
+ gl.ylabels_right = False
+
+ # Coastlines
+ if self.drawCoastLines and 'GeoAxes' in str(type(ax)):
+ ax.coastlines(resolution='10m')
+
+ # Countries border
+ ax.add_feature(cfeature.BORDERS)
+ ax.add_feature(cfeature.LAKES, alpha=0.7)
+
+ def addWhitecm(self, Laddwhite, colormap_in, nb_level):
+ if Laddwhite:
+ color_map = cm.get_cmap(colormap_in, 256)
+ newcolor_map = color_map(np.linspace(0, 1, 256))
+ whites = np.array([1, 1, 1, 1]) #RBG code + opacity
+ for i in range(int(256/nb_level)): newcolor_map[:i, :] = whites
+ newcmp = ListedColormap(newcolor_map)
+ return newcmp
+
+
+ def set_Title(self, ax, i, title, Lid_overlap, xlab, ylab):
+ """
+ Handle top title of each graph
+ Parameters :
+ - titlepad : global attribute for vertical shift placement w.r.t the graph
+ """
+ self.ax = ax
+ self.title = title
+ self.Lid_overlap = Lid_overlap
+ self.i = i
+ #self.ax[self.i].set_xlabel("test", fontweight='bold')
+ if not Lid_overlap:
+ self.ax[self.i].set_title(title, pad=self.titlepad)
+ else: # If graph overlap, title is concatenated
+ new_title = self.ax[self.i].get_title() + ' and ' + title
+ self.ax[self.i].set_title(new_title, pad=self.titlepad)
+
+ def set_xlim(self, ax, i, xlim):
+ """
+ Handle x limits plotted if necessary
+ """
+ self.ax = ax
+ self.xlim = xlim
+ self.i = i
+
+ self.ax[self.i].set_xlim(xlim[0],xlim[1])
+
+ def set_ylim(self, ax, i, ylim):
+ """
+ Handle x limits plotted if necessary
+ """
+ self.ax = ax
+ self.ylim = ylim
+ self.i = i
+
+ self.ax[self.i].set_ylim(ylim[0],ylim[1])
+
+ def set_XYaxislab(self, ax, i, xlab, ylab):
+ """
+ Handle x and y axis labels
+ """
+ self.ax = ax
+ self.xlab = xlab
+ self.ylab = ylab
+ self.i = i
+
+ # This handing label is a known issue with GeoAxes of cartopy
+ # https://stackoverflow.com/questions/35479508/cartopy-set-xlabel-set-ylabel-not-ticklabels
+ # https://github.com/SciTools/cartopy/issues/1332
+ if 'GeoAxes' in str(type(self.ax[self.i])):
+ self.ax[self.i].text(-0.11, 0.45, ylab, verticalalignment='top', horizontalalignment='left',
+ rotation='vertical', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=11)
+ self.ax[self.i].text(0.45, -0.06, xlab, verticalalignment='top', horizontalalignment='left',
+ rotation='horizontal', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=11)
+ else:
+ self.ax[self.i].set_xlabel(xlab)
+ self.ax[self.i].set_ylabel(ylab)
+
+ def addLine(self, ax, beg_coord, end_coord, color='black', linewidth=1):
+ self.ax = ax
+ self.beg_coord = beg_coord
+ self.end_coord = end_coord
+ self.color = color
+ self.linewidth = linewidth
+
+ x1, y1 = [self.beg_coord[0],self.end_coord[0]], [self.beg_coord[1], self.end_coord[1]]
+ ax.plot(x1, y1,color=self.color,linewidth=self.linewidth)
+
+ def print_minmax(self, var, title):
+ print(str(title) + " min = " + str(np.nanmin(var)) + " max = " + str(np.nanmax(var)))
+
+ def set_minmaxText(self, ax, i, var, title, Lid_overlap, facconv):
+ """
+ Show min and max value Text in the plot
+ If overlap variable, text is align to the right
+ TODO : handle more than 2 overlap variables
+ """
+ self.ax = ax
+ self.var = var
+ self.i = i
+ self.title = title
+ self.facconv = facconv
+
+ strtext = " min = " + "{:.3e}".format(np.nanmin(var*facconv)) + " max = " + "{:.3e}".format(np.nanmax(var*facconv))
+ if not Lid_overlap:
+ self.ax[self.i].text(0.01, self.minmaxpad, strtext, verticalalignment='top', horizontalalignment='left',
+ transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+ else:
+ self.ax[self.i].text(self.lateralminmaxpad, self.minmaxpad, strtext, verticalalignment='top', horizontalalignment='right',
+ transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+ # Print to help choose min/max value for ticks
+ self.print_minmax(var*facconv, title)
+
+ def showTimeText(self, ax, i, timetxt):
+ """
+ Show time validity
+ """
+ self.ax = ax
+ self.i = i
+ self.timetxt = timetxt
+
+ strtext = "Time = " + timetxt
+ self.ax[self.i].text(0.0, self.timepad, strtext, verticalalignment='top', horizontalalignment='left',
+ transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+
+
+ def psectionV(self, Lxx=[], Lzz=[], Lvar=[], Lxlab=[], Lylab=[], Ltitle=[], Lminval=[], Lmaxval=[],
+ Lstep=[], Lstepticks=[], Lcolormap=[], Lcbarlabel=[], LcolorLine=[],
+ Lfacconv=[], ax=[], Lid_overlap=[], colorbar=True, orog=[], Lxlim=[], Lylim=[], Ltime=[], Lpltype=[], LaddWhite_cm=[]):
+ """
+ Horizontal cross section plot
+ Arguments :
+ - Lxx : List of x or y coordinate variable or time axis
+ - Lzz : List of z coordinates variable
+ - Lvar : List of variables to plot
+ - Lxlab : List of x-axis label
+ - Lylab : List of y-axis label
+ - Lxlim : List of x (min, max) value plotted
+ - Lylim : List of y (min, max) value plotted
+ - Ltime : List of time (validity)
+ - Ltitle : List of sub-title
+ - Lminval: List of minimum value for each colorbar
+ - Lmaxval: List of maximum value for each colorbar
+ - Lstep : List of color-steps for each colorbar
+ - Lstepticks : List of value of labels for each colorbar
+ - Lcolormap : List of colormap
+ - LcolorLine : List of colors for colors arg of contour (color line only)
+ - Lcbarlabel : List of colorbar label legend (units)
+ - Lfacconv : List of factors for unit conversion of each variables
+ - ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
+ - Lid_overlap: List of number index of plot to overlap current variables
+ - Lpltype : List of types of plot 'cf' or 'c'. cf=contourf, c=contour (lines only)
+ - colorbar : show colorbar or not
+ - LaddWhite_cm : List of boolean to add white color to a colormap at the first (low value) tick colorbar
+ """
+ self.ax = ax
+ firstCall = (len(self.ax) == 0)
+
+ # Initialize default value w.r.t to the number of plots
+# D={'Lxlab':Lxlab, 'Lylab':Lylab, 'Ltitle':Ltitle,'Lminval':Lminval, 'Lmaxval':Lmaxval,
+# 'Lstep':Lstep, 'Lstepticks':Lstepticks, 'Lcolormap':Lcolormap, 'Lcbarlabel':Lcbarlabel, 'Lfacconv':Lfacconv, 'Ltime':Ltime,
+# 'LaddWhite_cm':LaddWhite_cm, 'Lpltype':Lpltype}
+# D = initialize_default_val(Lvar, D)
+
+ # Default values
+ if not Lfacconv: Lfacconv = [1.0]*len(Lvar)
+ if not Lcolormap and not LcolorLine: Lcolormap=['gist_rainbow_r']*len(Lvar) # If no color given, a cmap is given
+ if not Lcolormap: LcolorLine=['black']*len(Lvar)
+ if not Lpltype: Lpltype=['cf']*len(Lvar)
+ if not LaddWhite_cm: LaddWhite_cm=[False]*len(Lvar)
+ if not Lylab: Lylab = ['']*len(Lvar)
+ if not Lxlab: Lxlab = ['']*len(Lvar)
+ # Add an extra percentage of the top max value for forcing the colorbar show the true user maximum value (correct a bug)
+ Lmaxval = list(map(lambda x, y: x + 1E-6*y, Lmaxval, Lstep) ) #The extra value is 1E-6 times the step ticks of the colorbar
+
+ # On all variables to plot
+ for i,var in enumerate(Lvar):
+ if firstCall: #1st call
+ iax = i
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1))
+ self.nb_graph+=1
+ elif Lid_overlap != []: #overlapping plot
+ iax = Lid_overlap[i]
+ else: #existing ax with no overlapping (graphd appended to existing panel)
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1))
+ self.nb_graph+=1
+ iax = len(self.ax)-1 # The ax index of the new coming plot is the length of the existant ax -1 for indices matter
+
+ # Colors normalization
+ norm = mpl.colors.Normalize(vmax=Lmaxval[i], vmin=Lminval[i])
+
+ # Print min/max (stout and on plot)
+ self.set_minmaxText(self.ax, iax, var, Ltitle[i], Lid_overlap, Lfacconv[i])
+
+ # Print time validity
+ if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
+
+ # Number of contours level
+ if not Lstep[i]: # Default value of number of steps is 20
+ Lstep[i] = (Lmaxval[i] - Lminval[i])/20
+ Lstepticks[i] = Lstep[i]
+
+ levels_contour = np.arange(Lminval[i],Lmaxval[i],step=Lstep[i])
+
+ # Add White to colormap
+ if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(LaddWhite_cm[i], Lcolormap[i], len(levels_contour))
+
+ # Plot
+ if Lpltype[i]=='c': # Contour
+ if LcolorLine:
+ cf = self.ax[iax].contour(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], colors=LcolorLine[i])
+ else:
+ cf = self.ax[iax].contour(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+ else: # Contourf
+ cf = self.ax[iax].contourf(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+
+ # Title
+ self.set_Title(self.ax, iax, Ltitle[i], Lid_overlap,Lxlab[i], Lylab[i])
+
+ # X/Y Axis label
+ self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ # X/Y Axis limits value
+ if Lxlim:
+ try:
+ self.set_xlim(self.ax, iax, Lxlim[i])
+ except:
+ pass
+ if Lylim:
+ try:
+ self.set_ylim(self.ax, iax, Lylim[i])
+ except:
+ pass
+
+ # Color label on contour-ligne
+ if Lpltype[i]=='c': # Contour
+ self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i]))
+
+ #Filling area under topography
+ if not orog==[]:
+ self.ax[iax].fill_between(Lxx[i][0,:], orog, color='black')
+
+ # Colorbar
+ if colorbar:
+ cb=plt.colorbar(cf, ax=self.ax[iax], fraction=0.031, pad=self.colorbarpad, ticks=np.arange(Lminval[i],Lmaxval[i],Lstepticks[i]), aspect=self.colorbaraspect)
+ cb.ax.set_title(Lcbarlabel[i], pad = self.labelcolorbarpad, loc='left') #This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
+
+ return self.fig
+
+ def pXY_lines(self, Lzz=[], Lvar=[], Lxlab=[], Lylab=[], Ltitle=[], Llinetype=[], Llinewidth=[],
+ Llinecolor=[], Llinelabel=[], Lfacconv=[], ax=[], id_overlap=None, Lxlim=[],
+ Lylim=[], Ltime=[], LaxisColor=[]):
+ """
+ XY (multiple)-lines plot
+ Arguments :
+ - Lzz : List of z coordinates variable [1D]
+ - Lvar : List of variables to plot [1D]
+ - Lxlab : List of x-axis label
+ - Lylab : List of y-axis label
+ - Lxlim : List of x (min, max) value plotted
+ - Lylim : List of y (min, max) value plotted
+ - Ltime : List of time (validity)
+ - Ltitle : List of sub-title
+ - Llinewidth : List of lines thickness
+ - Llinetype : List of line types
+ - Lcolorlines: List of color lines
+ - Llinelabel : List of legend label lines
+ - Lfacconv : List of factors for unit conversion of each variables
+ - ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
+ - Lid_overlap: List of number index of plot to overlap current variables
+ - LaxisColor : List of colors for multiple x-axis overlap
+ """
+ self.ax = ax
+ firstCall = (len(self.ax) == 0)
+ # Defaults value convert to x number of variables list
+ if not Lfacconv: Lfacconv = [1.0]*len(Lvar)
+ if not Llinewidth: Llinewidth = [1.0]*len(Lvar)
+ if not Llinecolor: Llinecolor = ['blue']*len(Lvar)
+ if not Llinetype: Llinetype = ['-']*len(Lvar)
+ if not Llinelabel: Llinelabel = ['']*len(Lvar)
+ if not LaxisColor: LaxisColor = ['black']*len(Lvar)
+ if not Lylab: Lylab = ['']*len(Lvar)
+ if not Lxlab: Lxlab = ['']*len(Lvar)
+
+ if firstCall: # 1st call
+ iax = 0
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,1, label='graph axe x down'))
+ self.nb_graph+=1
+ elif id_overlap: # overlapping plot with a different x-axis
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph, label='graph axe x top',frame_on=False))
+ iax = len(self.ax)-1
+ else: # existing ax with no overlapping (graph appended to existing panel)
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1))
+ self.nb_graph+=1
+ iax = len(self.ax)-1 # The ax index of the new coming plot is the length of the existant ax -1 for indices matter
+
+ # On all variables to plot
+ for i,var in enumerate(Lvar):
+ # Print time validity
+ if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
+
+ # Plot
+ cf = self.ax[iax].plot(var*Lfacconv[i], Lzz[i], color=Llinecolor[i], ls=Llinetype[i],
+ label=Llinelabel[i], linewidth=Llinewidth[i])
+ # Legend
+ #TODO : Handling legend with overlap two axis lines in the same box. For now, placement is by hand
+ if not id_overlap:
+ self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.95))
+ else:
+ self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.90))
+
+ # Title
+ if Ltitle: self.set_Title(self.ax, iax, Ltitle[i], id_overlap,Lxlab[i], Lylab[i])
+
+ # X/Y Axis label
+ if id_overlap:
+ self.ax[iax].xaxis.tick_top()
+ self.ax[iax].xaxis.set_label_position('top')
+ self.ax[iax].set_xlabel(Lxlab[i])
+ self.ax[iax].set_ylabel(Lylab[i])
+ else:
+ self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ self.ax[iax].tick_params(axis='x', colors=LaxisColor[i])
+
+ # X/Y Axis limits value
+ if Lxlim:
+ try:
+ self.set_xlim(self.ax, iax, Lxlim[i])
+ except:
+ pass
+ if Lylim:
+ try:
+ self.set_ylim(self.ax, iax, Lylim[i])
+ except:
+ pass
+ return self.fig
+
+ def psectionH(self, lon=[],lat=[], Lvar=[], Lcarte=[], Llevel=[], Lxlab=[], Lylab=[], Ltitle=[], Lminval=[], Lmaxval=[],
+ Lstep=[], Lstepticks=[], Lcolormap=[], LcolorLine=[], Lcbarlabel=[], Lproj=[], Lfacconv=[], coastLines=True, ax=[],
+ Lid_overlap=[], colorbar=True, Ltime=[], LaddWhite_cm=[], Lpltype=[]):
+ """
+ Horizontal cross section plot
+ Arguments :
+ - lon : longitude 2D array
+ - lat : latitude 2D array
+ - Lvar : List of variables to plot
+ - Lcarte : Zooming [lonmin, lonmax, latmin, latmax]
+ - Llevel : List of k-level value for the section plot (ignored if variable is already 2D)
+ - Lxlab : List of x-axis label
+ - Lylab : List of y-axis label
+ - Ltitle : List of sub-title
+ - Ltime : List of time (validity)
+ - Lminval: List of minimum value for each colorbar
+ - Lmaxval: List of maximum value for each colorbar
+ - Lstep : List of color-steps for each colorbar
+ - Lstepticks : List of value of labels for each colorbar
+ - Lcolormap : List of colormap
+ - LcolorLine : List of colors for colors arg of contour (color line only)
+ - Lcbarlabel : List of colorbar label legend (units)
+ - Lproj : List of ccrs cartopy projection ([] for cartesian coordinates)
+ - Lfacconv : List of factors for unit conversion of each variables
+ - coastLines : Boolean to plot coast lines and grid lines
+ - ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
+ - Lid_overlap : List of number index of plot to overlap current variables
+ - colorbar : show colorbar or not
+ - Lpltype : List of types of plot 'cf' or 'c'. cf=contourf, c=contour (lines only)
+ - LaddWhite_cm : List of boolean to add white color to a colormap at the first (low value) tick colorbar
+ """
+ self.ax = ax
+ firstCall = (len(self.ax) == 0)
+
+ # Initialize default value w.r.t to the number of plots
+# D={'lon':lon, 'lat':lat, 'Lcarte':Lcarte, 'Llevel':Llevel, 'Lxlab':Lxlab, 'Lylab':Lylab, 'Ltitle':Ltitle,'Lminval':Lminval, 'Lmaxval':Lmaxval,
+# 'Lstep':Lstep, 'Lstepticks':Lstepticks, 'Lcolormap':Lcolormap, 'Lcbarlabel':Lcbarlabel, 'Lproj':Lproj, 'Lfacconv':Lfacconv, 'Ltime':Ltime,
+# 'LaddWhite_cm':LaddWhite_cm, 'Lpltype':Lpltype}
+# D = initialize_default_val(Lvar, D)
+
+ # If all plots are not using conversion factor, convert it to List
+ if not Lfacconv: Lfacconv = [1.0]*len(Lvar)
+ if not Lylab: Lylab = ['']*len(Lvar)
+ if not Lxlab: Lxlab = ['']*len(Lvar)
+ if not Lcolormap and not LcolorLine: Lcolormap=['gist_rainbow_r']*len(Lvar) # If no color given, a cmap is given
+ if not Lcolormap: LcolorLine=['black']*len(Lvar)
+ if not Lpltype: Lpltype=['cf']*len(Lvar)
+ if not LaddWhite_cm: LaddWhite_cm=[False]*len(Lvar)
+ # Add an extra percentage of the top max value for forcing the colorbar show the true user maximum value (correct a bug)
+ if Lstep: Lmaxval = list(map(lambda x, y: x + 1E-6*y, Lmaxval, Lstep) ) #The extra value is 1E-6 times the step ticks of the colorbar
+
+ # This following must be after the extra percentage top value addition
+ if not Lstep: Lstep=[None]*len(Lvar)
+ if not Lstepticks: Lstepticks=[None]*len(Lvar)
+
+ # On all variables to plot
+ for i,var in enumerate(Lvar):
+ if firstCall: #1st call
+ iax = i
+ if Lproj:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1, projection = Lproj[i]))
+ else: # Cartesian coordinates
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1))
+ self.nb_graph+=1
+ elif Lid_overlap != []: #overlapping plot
+ iax = Lid_overlap[i]
+ else: #existing ax with no overlapping (graphd appended to existing panel)
+ if Lproj:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1, projection = Lproj[i]))
+ else: # Cartesian coordinates
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1))
+ self.nb_graph+=1
+ iax = len(self.ax)-1 # The ax index of the new coming plot is the length of the existant ax -1 for indices matter
+
+ #Colors normalization
+ norm = mpl.colors.Normalize(vmax=Lmaxval[i], vmin=Lminval[i])
+
+ # Zooming
+ if len(Lcarte) == 4: #zoom
+ self.ax[iax].set_xlim(Lcarte[0], Lcarte[1])
+ self.ax[iax].set_ylim(Lcarte[2], Lcarte[3])
+ # Variable to plot w.r.t dimensions
+ if var.ndim==2:
+ vartoPlot = var[:,:]
+ else:
+ vartoPlot = var[Llevel[i],:,:]
+
+ # Print min/max (stout and on plot)
+ self.set_minmaxText(self.ax, iax, vartoPlot, Ltitle[i], Lid_overlap, Lfacconv[i])
+
+ # Print time validity
+ if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
+
+ # Number of contours level
+ if not Lstep[i]: # Default value of number of steps is 20
+ Lstep[i] = (Lmaxval[i] - Lminval[i])/20
+ Lstepticks[i] = Lstep[i]
+
+ levels_contour = np.arange(Lminval[i],Lmaxval[i],step=Lstep[i])
+
+ # Add White to colormap
+ if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(LaddWhite_cm[i], Lcolormap[i], len(levels_contour))
+
+ # Plot
+ if Lproj:
+ if Lpltype[i]=='c': # Contour
+ if LcolorLine:
+ cf = self.ax[iax].contourf(lon[i], lat[i], vartoPlot*Lfacconv[i], levels=levels_contour,transform=Lproj[i],
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], colors=LcolorLine[i])
+ else:
+ cf = self.ax[iax].contourf(lon[i], lat[i], vartoPlot*Lfacconv[i], levels=levels_contour,transform=Lproj[i],
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+ else:
+ cf = self.ax[iax].contourf(lon[i], lat[i], vartoPlot*Lfacconv[i], levels=levels_contour,transform=Lproj[i],
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+ else: # Cartesian coordinates
+ if Lpltype[i]=='c': # Contour
+ cf = self.ax[iax].contourf(lon[i], lat[i], vartoPlot*Lfacconv[i], levels=levels_contour,
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], colors=LcolorLine[i])
+ else:
+ cf = self.ax[iax].contourf(lon[i], lat[i], vartoPlot*Lfacconv[i], levels=levels_contour,
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+ # Title
+ self.set_Title(self.ax, iax, Ltitle[i], Lid_overlap,Lxlab[i], Lylab[i])
+
+ # Coastlines / Grid lines and labels
+ if Lproj: self.draw_Backmap(coastLines, self.ax[iax], Lproj[i])
+
+ # X/Y Axis
+ self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ # Color label on contour-ligne
+ if Lpltype[i]=='c': # Contour
+ self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i]))
+
+ # Colorbar
+ if colorbar:
+ cb=plt.colorbar(cf, ax=self.ax[iax], fraction=0.031, pad=self.colorbarpad, ticks=np.arange(Lminval[i],Lmaxval[i],Lstepticks[i]), aspect=self.colorbaraspect)
+ cb.ax.set_title(Lcbarlabel[i], pad = self.labelcolorbarpad, loc='left') #This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
+
+ return self.fig
+
+
+ def pvector(self, Lxx=[], Lyy=[], Lvar1=[], Lvar2=[], Lcarte=[], Llevel=[], Lxlab=[], Lylab=[],
+ Ltitle=[], Lwidth=[], Larrowstep=[], Lcolor=[], Llegendval=[], Lcbarlabel=[],
+ Lproj=[], Lfacconv=[], ax=[], coastLines=True, Lid_overlap=[], Ltime=[], Lscale=[],
+ Lylim=[], Lxlim=[]):
+ """
+ Horizontal vectors lines
+ Arguments :
+ - Lxx : List of x or y coordinate variable (lat or ni or nm)
+ - Lyy : List of y coordinates variable (lon or level)
+ - Lvar1 : List of wind-component along x/y or oblic axis (3D for hor. section, 2D for vertical section)
+ - Lvar2 : List of wind-component along y-axis : v-component for horizontal section / w-component for vertical section
+ - Lcarte : Zooming [lonmin, lonmax, latmin, latmax]
+ - Llevel : List of k-level value for the horizontal section plot (ignored if variable is already 2D)
+ - Lxlab : List of x-axis label
+ - Lylab : List of y-axis label
+ - Lxlim : List of x (min, max) value plotted
+ - Lylim : List of y (min, max) value plotted
+ - Ltitle : List of sub-titles
+ - Ltime : List of time (validity)
+ - Lwidth : List of thickness of the arrows
+ - Lscale : List of scale for the length of the arrows (high value <=> small length)
+ - Larrowstep : List of sub-sample (frequency) if too much arrows
+ - Lcolor : List of colors for the arrows (default: black)
+ - Llegendval : List of value for the legend of the default arrow
+ - Lcbarlabel : List of labels for the legend of the default arrow
+ - Lproj : List of ccrs cartopy projection
+ - Lfacconv : List of factors for unit conversion of each variables
+ - coastLines : Boolean to plot coast lines and grid lines
+ - ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
+ - Lid_overlap : List of number index of plot to overlap current variables
+ """
+ self.ax = ax
+ firstCall = (len(self.ax) == 0)
+ # If all plots are not using conversion factor, convert it to List
+ if not Lfacconv: Lfacconv = [1.0]*len(Lvar1)
+ if not Lcolor: Lcolor = ['black']*len(Lvar1)
+ if not Lscale : Lscale = [None]*len(Lvar1)
+ if not Lylab: Lylab = ['']*len(Lvar1)
+ if not Lxlab: Lxlab = ['']*len(Lvar1)
+ # On all variables to plot
+ for i,var1 in enumerate(Lvar1):
+ if firstCall: #1st call
+ iax = i
+ if Lproj:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1, projection = Lproj[i]))
+ else:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1))
+ self.nb_graph+=1
+ elif Lid_overlap != []: #overlapping plot
+ iax = Lid_overlap[i]
+ else: #existing ax with no overlapping (graphd appended to existing panel)
+ if Lproj:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1, projection = Lproj[i]))
+ else:
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1))
+ self.nb_graph+=1
+ iax = len(self.ax)-1 # The ax index of the new coming plot is the length of the existant ax -1 for indices matter
+
+ # Zooming
+ if len(Lcarte) == 4: #zoom
+ self.ax[iax].set_xlim(Lcarte[0], Lcarte[1])
+ self.ax[iax].set_ylim(Lcarte[2], Lcarte[3])
+
+ # Variable to plot w.r.t dimensions
+ if var1.ndim==2:
+ vartoPlot1 = var1[:,:]
+ vartoPlot2 = Lvar2[i][:,:]
+ else: # Variable is 3D : only for horizontal section
+ vartoPlot1 = var1[Llevel[i],:,:]
+ vartoPlot2 = Lvar2[i][Llevel[i],:,:]
+
+ # Print min/max val to help choose colorbar steps
+ self.set_minmaxText(self.ax, iax, np.sqrt(vartoPlot1**2 + vartoPlot2**2), Ltitle[i], Lid_overlap, Lfacconv[i])
+
+ # Print time validity
+ if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
+
+ # Plot
+ axeX = Lxx[i]
+ axeY = Lyy[i]
+ if Lxx[i].ndim == 2:
+ cf = self.ax[iax].quiver(axeX[::Larrowstep[i],::Larrowstep[i]], axeY[::Larrowstep[i],::Larrowstep[i]],
+ vartoPlot1[::Larrowstep[i],::Larrowstep[i]], vartoPlot2[::Larrowstep[i],::Larrowstep[i]],
+ width=Lwidth[i] ,angles='uv', color=Lcolor[i],scale=Lscale[i])
+ else:
+ cf = self.ax[iax].quiver(axeX[::Larrowstep[i]], axeY[::Larrowstep[i]],
+ vartoPlot1[::Larrowstep[i],::Larrowstep[i]], vartoPlot2[::Larrowstep[i],::Larrowstep[i]],
+ width=Lwidth[i] ,angles='uv', color=Lcolor[i], scale=Lscale[i])
+ # Title
+ self.set_Title(self.ax, iax, Ltitle[i], Lid_overlap,Lxlab[i], Lylab[i])
+
+ # X/Y Axis Label
+ self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ # X/Y Axis limits value
+ if Lxlim:
+ try:
+ self.set_xlim(self.ax, iax, Lxlim[i])
+ except:
+ pass
+ if Lylim:
+ try:
+ self.set_ylim(self.ax, iax, Lylim[i])
+ except:
+ pass
+
+ # Coastlines:
+ if Lproj: self.draw_Backmap(coastLines, self.ax[iax], Lproj[i])
+
+ # Arrow legend key
+ qk = self.ax[iax].quiverkey(cf, 1.0, -0.05, Llegendval[i], str(Llegendval[i]) + Lcbarlabel[i], labelpos='E', color='black')
+
+ return self.fig
+
+ def pXY_bar(self, Lbins=[], Lvar=[], Lxlab=[], Lylab=[], Ltitle=[], Lcolor=[], Lwidth=[],
+ Llinecolor=[], Llinewidth=[], Lfacconv=[], ax=[], id_overlap=None, Lxlim=[],
+ Lylim=[], Ltime=[], LaxisColor=[]):
+ """
+ XY Histogram
+ Arguments :
+ - Lbins : List of bins
+ - Lvar : List of the value for each bin
+ - Lxlab : List of x-axis label
+ - Lylab : List of y-axis label
+ - Ltitle : List of sub-title
+ - Lcolor : List of color (or sequence of colors for each value) of the bars
+ - Lwidth : List of width of the bars
+ - Llinecolor : List of line color of the bar edges
+ - Llinewidth : List of lines thickness of the bar edges
+ - Lfacconv : List of factors for unit conversion of each variables
+ - Lxlim : List of x (min, max) value plotted
+ - Lylim : List of y (min, max) value plotted
+ - Ltime : List of time (validity)
+ - ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
+ - Lid_overlap: List of number index of plot to overlap current variables
+ - LaxisColor : List of colors for multiple x-axis overlap
+ """
+ self.ax = ax
+ firstCall = (len(self.ax) == 0)
+ # Defaults value convert to x number of variables list
+ if not Lfacconv: Lfacconv = [1.0]*len(Lvar)
+ if not LaxisColor: LaxisColor = ['black']*len(Lvar)
+ if not Lylab: Lylab = ['']*len(Lvar)
+ if not Lxlab: Lxlab = ['']*len(Lvar)
+ if not Lcolor: Lcolor = ['black']*len(Lvar)
+ if not Lwidth: Lwidth = [1]*len(Lvar)
+ if not Llinecolor: Llinecolor = ['black']*len(Lvar)
+ if not Llinewidth: Llinewidth = [0]*len(Lvar)
+
+ # On all variables to plot
+ for i,var in enumerate(Lvar):
+ if firstCall: # 1st call
+ iax = i
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,i+1, label='graph axe x down'))
+ self.nb_graph+=1
+ elif id_overlap: # overlapping plot with a different x-axis
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph, label='graph axe x top',frame_on=False))
+ iax = len(self.ax)-1
+ else: # existing ax with no overlapping (graph appended to existing panel)
+ self.ax.append(self.fig.add_subplot(self.nb_l,self.nb_c,self.nb_graph+1))
+ self.nb_graph+=1
+ print(len(self.ax))
+ iax = len(self.ax)-1 # The ax index of the new coming plot is the length of the existant ax -1 for indices matter
+
+ # Print time validity
+ if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
+
+ # Bins by labels
+ labels = np.array([str(L) for L in Lbins[i]])
+
+ # Plot
+ cf = self.ax[iax].bar(labels, var*Lfacconv[i], width=Lwidth[i], color=Lcolor[i], linewidth=Llinewidth[i], edgecolor=Llinecolor[i])
+
+ # Legend
+ #TODO : Handling legend with overlap two axis lines in the same box. For now, placement is by hand
+ if not id_overlap:
+ self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.95))
+ else:
+ self.ax[iax].legend(loc='upper right', bbox_to_anchor=(1, 0.90))
+
+ # Title
+ if Ltitle: self.set_Title(self.ax, iax, Ltitle[i], id_overlap,Lxlab[i], Lylab[i])
+
+ # X/Y Axis label
+ if id_overlap:
+ self.ax[iax].xaxis.tick_top()
+ self.ax[iax].xaxis.set_label_position('top')
+ self.ax[iax].set_xlabel(Lxlab[i])
+ self.ax[iax].set_ylabel(Lylab[i])
+ else:
+ self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ self.ax[iax].tick_params(axis='x', colors=LaxisColor[i])
+
+ # X/Y Axis limits value
+ if Lxlim:
+ try:
+ self.set_xlim(self.ax, iax, Lxlim[i])
+ except:
+ pass
+ if Lylim:
+ try:
+ self.set_ylim(self.ax, iax, Lylim[i])
+ except:
+ pass
+
+ return self.fig
+#def initialize_default_val(Lvar, Dparam):
+# #TODO : initialize default value for all parameter of all type of graphs
+# #Returns All the parameters given in Dparam where :
+# "- If no value is found (empty list []) : return the default value (if exist) * number of graph
+# #- If ONE value only is found : return the value copied x times the number of graph
+# #- If the list is complete, nothing is done
+# #CURRENT PROBLEM
+# #The returned value do not change the true referenced variable given as argument
+# # Number of graphs
+# l = len(Lvar)
+#
+# # Add an extra percentage of the top max value for forcing the colorbar show the true user maximum value (correct a bug)
+# #if Dparam['Lstep'] and Dparam['Lmaxval']: Dparam['Lmaxval'] = list(map(lambda x, y: x + 1E-6*y, Dparam['Lmaxval'], Dparam['Lstep']) ) #The extra value is 1E-6 times the step ticks of the colorbar
+# print(Dparam.items())
+#
+# # Read on which parameters initialize the default values
+# for args_t in list(Dparam.items()): # Test on all arguments present, if they are empty list, default values apply for each graph
+# args = list(args_t)
+# print(args)
+# if args[0] == 'Lfacconv' and not args[1]: args[1] = [1.0]*l
+# elif args[0] == 'Lcolormap' and not args[1]: args[1] = ['gist_rainbow_r']*l
+# elif args[0] == 'LcolorLine' and not args[1]: args[1] = ['black']*l
+# elif args[0] == 'Lpltype' and not args[1]: args[1]= ['cf']*l
+# elif args[0] == 'LaddWhite_cm' and not args[1]: args[1] = ['False']*l
+# elif args[0] == 'Lstep' and not args[1]: args[1] = [None]*l # default value filled later
+# elif args[0] == 'Lstepticks' and not args[1]: args[1] = [None]*l
+# Dparam[args[0]] = args[1]
+# print(args)
+#
+# # Check if there is no value for a parameter
+## for args_t in list(Dparam.items()):
+## args = list(args_t)
+## if args[1]
+# # Number of contours level
+## for i in range(l):
+## if 'Lstepticks' in arguments.args and 'Lmaxval' in arguments.args and 'Lminval' in arguments.args:
+## Lstep[i] = (Lmaxval[i] - Lminval[i])/20 # Default value of number of steps is 20
+## Lstepticks[i] = Lstep[i] # Default value is stepticks are the same as steps values
+#
+#
+# return Dparam
\ No newline at end of file
diff --git a/src/LIB/Python/misc_functions.py b/src/LIB/Python/misc_functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..e1e467a767691b327d4030445820299d6f77a3ca
--- /dev/null
+++ b/src/LIB/Python/misc_functions.py
@@ -0,0 +1,135 @@
+#!/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
diff --git a/src/LIB/Python/read_MNHfile.py b/src/LIB/Python/read_MNHfile.py
new file mode 100644
index 0000000000000000000000000000000000000000..4279599fb48fb9c67ffc640c157c3088462d163f
--- /dev/null
+++ b/src/LIB/Python/read_MNHfile.py
@@ -0,0 +1,208 @@
+#!/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