misc_functions.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
MNH_LIC for details. version 1.

@author: 07/2021 Quentin Rodier
"""
import copy
from scipy.interpolate import RectBivariateSpline
import numpy as np
import math

def convert_date(datesince, 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, level, angle):
    """Compute the projected horizontal wind vector on an axis with a given angle w.r.t. the x/ni axes (West-East)
    
    Parameters
    ----------
    u : array 3D
        U-wind component
    
    v : array 3D
        V-wind component

    level : array 1D
        level dimension array
        
    angle : float
        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
    -------
    
    projected_wind : array 3D
        a 3D wind component projected on the axe to be used with Panel_Plot.pvector as Lvar1
    """ 
    projected_wind = copy.deepcopy(u)
    for k in range(len(level)):
        projected_wind[k,:,:] = u[k,:,:]*math.cos(angle) + v[k,:,:]*math.sin(angle)
    return projected_wind

def oblique_proj(var, ni, nj, lvl, i_beg, j_beg, i_end, j_end):
    """Compute an oblique projection of a variable w.r.t. its axes
    
    Parameters
    ----------
    var : array 3D or 2D
        the variable to project (e.g. THT, ZS)
    
    ni : array 1D
        1D x-axis of the 3D dimension
    
    nj : array 1D
        1D y-axis of the 3D dimension
    
    level : array 1D
        1D z-axe of the 3D dimension
    
    i_beg, j_beg : int
        coordinate of the begin point of the new axe
    
    i_end, j_end : int
        coordinate of the end point of the new axe       
    
     Returns
     -------
    angle_proj : float
        the angle (radian) of the new axe w.r.t the x/ni axes (West-East)
         
    out_var : array 2D or 1D
        a 2D (z,m) or 1D (m) variable projected on the oblique axe
    
    axe_m : array 1D
        a 1D m new axe (distance from the beggining point)
    
    """
    dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) # Distance de la section oblique  m
    if var.ndim ==3:
        out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
    else:  # 2D
        out_var = np.zeros(int(dist_seg)+1) # Initialisation du nouveau champs projeté dans la coupe (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)

    if var.ndim ==3: # 3D variables to project
        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é
    else:  # 2D variables to project
        a=RectBivariateSpline(ni, nj,var[:,:],kx=1,ky=1)
        for m in range(int(dist_seg)+1):
            out_var[m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1])
    
    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
    
    Parameters
    ----------
    oneVar2D : array 2D
        a 2D array (e.g. UT, THT)
        
    orography : array 1D
        1D orography (ZS)
        
    ztop : real
        scalar of the top height of the model (ZTOP)
        
    level : array 1D
        1D level variable (level or level_w)
        
    n_xory : array 1D  
        1D directionnal grid variable (ni_u, nj_u, ni_v or nj_v)
        
    Returns
    -------
    altitude
        a 2D altitude variable with topography taken into account
        
    n_xory_2D
        a 2D directionnal variable duplicated from n_xory
    """
    n_xory_2D = copy.deepcopy(oneVar2D)
    altitude  =  copy.deepcopy(oneVar2D)
    
    for k in range(len(level)):
        n_xory_2D[k,:] = 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
    
    Parameters
    ----------
    oneVar3D : array 3D
        a 3D array (e.g. UT, THT)
        
    orography : array 2D
        2D orography (ZS)
        
    ztop : real
        scalar of the top height of the model (ZTOP)
        
    level : array 1D
        1D level variable (level or level_w)
        
    n_x : array 1D  
        1D directionnal grid variable along i (ni_u, or ni_v)
        
    n_y : array 1D  
        1D directionnal grid variable along j (nj_u, or nj_v)
        
    Returns
    -------
    altitude
        a 3D altitude variable with topography taken into account
        
    n_x3D
        a 3D directionnal variable duplicated from n_x
        
    n_y3D
        a 3D directionnal variable duplicated from n_y
    """
    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]):
            if ztop==0:
                altitude[:,i,j] = level[:]
            else:
                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