Le programme lidar_ozone_actris_v01.py ====================================== Les grandes étapes du traitement lidar -------------------------------------- Généralités +++++++++++ Tout le programme a été concu autours de grandes étapes de traitement, elles mêmes sous classé en deux grandes parties. En effet, la donnée brute lidar est une grandeur physique compté par de détecteurs. Cette donnée est grandement impacté par l'instrument qui a réalisée la mesure. Le but de la première étape, appelée **prétraitement** consiste à supprimer tous les effets dus à l'instrument. On va donc retrouver la suppression du fond de ciel, la correction de la saturation des détecteurs et la correction du facteur géométrique. Globalement, la grandeur physique en sortie de ce premier traitement est similaire à celle d'entrée. La deuxième étape, appelée courament **traitement** consiste à transformer la grandeur physique mesurée en grandeur géophysique recherché. Par exemple, dans le cas d'un lidar ozone, on va transformer un nombre de photon reçu en concentration de molécule d'ozone. Le **prétraitement** est unique en fonction de l'instrument, en revanche, le **traitement** doit être identique en foncion de la donnée géophysique recherchée. Fond de ciel, bruit du détecteur et signal induit +++++++++++++++++++++++++++++++++++++++++++++++++ On appelle "fond de ciel" 3 phénomènes spécifiques : - **la lumière issue du ciel** : cela peut être la lumière du soleil, de la lune, des étoiles, de la lumière artificielle des villes, et du batiment le plus proche. Globalement, c'est un bruit blanc, dépendant de la longueur d'onde, mais non dépendant de l'altitude. Les variations en fonction du temps sont lentes, mais non nulles. - **le bruit de fond de détection** : de base, le détecteur et la carte d'acquisition ne sont pas parfait, et, même plongés dans le noor total, le détecteur 'compte' quelques photons de temps en temps. C'est en général un bruit blanc, dépendant du détecteur, mais non dépendant de l'altitude. En théorie, ce bruit n'évolue pas dans le temps. Il est cependant intéressant à vérifier de temps en temps, car le vieillissement de celui-ci peut avoir un impact sur cette valeur. - **le signal induit** : c'est un signal qui est dépendant de la lumière reçue. Globalement, c'est photons rentrant dans le détecteur sont captés par celui-ci et relachés plus tard. Ce phénomène peut être important lorsque le signal reçue sur les détecteurs est important. Cette fois ci, le signal n'est pas un bruit blanc car dépendant de l'altitude Dans la partie "fond de ciel", il s'agit de supprimer ces 3 phénomènes en même temps. En théorie, une simple moyenne est suffisant, mais du fait su signal induit, on va souvent utiliser une régression linéaire sur le fond du signal, afin de corriger la valeur de signal induit. Cette correction linéaire doit être faible, une forte pente de la régression indique un problème important de signal induit, soit un autre problème. .. figure:: ./figures/background.png :alt: background example :scale: 50 % :align: center Example of background correction Désaturation ++++++++++++ Le détecteur n'est pas parfait, et il peut être éblouit s'il recoit trop de signal. Cela veut dire qu'il a une valeur seuil de comptage. Cependant, avant d'être totalement éblouie, il y a une zone où il compte "mal" le nombre de photons, et plus particulièrement, où le nombre de photon comptés est infénieur au nombre de photons reçus. Il existe plusieurs équations mathématiques permettant de simuler cette saturation du détecteur, souvent paramétrisables, c'est-à-dire où il faut choisir un paramètre, qui peut évoluer avec le veillissement du détecteur. .. figure:: ./figures/desaturation_etudeMuller.jpg :alt: desaturation example :align: center Etude de la désaturation des signaux par la méthode de Muller Facteur géométrique +++++++++++++++++++ Le facteur géométrique correspond au taux de visualisation du laser par le télescope : lorsqu'il est égal à 0, le laser n'est pas dans le champ de vue du télescope, lorsqu'il est égale à 1, il est totalement, et entre les deux lorsque le laser est partiellement dans le champ de vue du télescope. .. figure:: ./figures/geometrical_factor.png :alt: desaturation example :align: center Exemple de facteur géométrique. Gluing des données ++++++++++++++++++ Cette technique permet de relier des signaux entre eux. En effet, les détecteurs n'était pas parfait, ils ne peuvent travailler que sur une gamme d'altitude. Pour couvrir une large gamme d'altitude, il faut donc utiliser plusieurs détecteurs, et donc relier les signaux par la suite. .. figure:: ./figures/gluing.png :alt: desaturation example :align: center Exemple de gluing des signaux, entre une voie haute (en orange) et une voie basse (en bleu). A gauche, on peut voir les signaux bruts, au milieu les signaux glués, et à droite le type de raccordement effectué (ici, un système en cosinus et sinus carrés. En plus du raccordement, il y a eu un décalage de la voie haute sur la voie basse. Filtrage des signaux ++++++++++++++++++++ Le filtrage des signaux lidars va permettre de réduire le bruit, et de pouvoir aller plus haut en altitude. Comme le bruit du issue du détecteur est dépendant de la puissance du signal, et donc de l'altitude, il est courant d'utiliser un filtrage dont le nombre de point du filtre (et donc sa résolution) dépend de l'altitude. Dans les basses altitudes, lorsque le signal est fort, on filtre peu, mais dans les hautes altitudes, lorsque le signal est faible, on filtre beaucoup plus. Le fait de filtrer va permettre de réduire le bruit au détriment de la résolution. En effet, pour filtrer, il faut utiliser un certain nombre de points, et donc le nombre de points indépendant va diminuer. .. figure:: ./figures/filtering.png :alt: desaturation example :align: center Exemple du filtrage d'un signal. A gauche, le signal filtré (en orange) est moins bruité sur le signal brut (en bleu). Au niveau de l'erreur, on voit que le signal filtrée à une erreur plus faible. Enfin, à droite, on voit la résolution, qui passe de 50 mètres à 4 km d'altitude à presque 400 mètres à 40 km d'altitude. Le traitement aérosols ++++++++++++++++++++++ Il s'agit de récupérer un profil de rétrodiffusion par la méthode Klett. Il s'applique seulement sur les voies Rayleigh Mie, à des longueurs d'onde non perturbées par d'autres molécules (comme la voie ON de l'ozone). On va souvent visualiser le signal en backscatter ratio, c'est-à-dire le rapport entre le signal diffusé par les aérosols et le signal rétrodiffusé total (aérosols et molécules de l'atmosphère). Lorsque ce signal est à 1, il n'y a pas d'aéorosls, et lorsque celui-ci est supérieur à 1, c'est qu'il y a présence d'aérosols. En dessous de 1, il y a un problème dans la mesure. .. figure:: ./figures/baskcatterratio.png :alt: desaturation example :align: center Exemple de backscater ratio Le traitement ozone +++++++++++++++++++ On part de l'équation lidar pour un signal Rayleigh : .. math:: P(\lambda_r, z) = K(\lambda_e, \lambda_r)*\frac{A}{(z-z_0)^2}*\beta(\lambda_e, \lambda_r, z)* exp \left( - \sum_i \left( \tau_{i}(\lambda_e, z) + \tau_{i}(\lambda_r, z)) \right) \right) avec : - :math:`P` : le signal reçu par le détecteur - :math:`\lambda_e` : longueur d'onde d'emission du laser - :math:`\lambda_r` : longueur d'onde de reception du détecteur - :math:`z` : altitude - :math:`K` : constante instrumentale - :math:`A` : angle solide - :math:`beta` : coefficient de diffusion - :math:`\tau_{i}` : transmission de la moléculue i, les molécules principales de l'atmosphère (N2 et O2), les aérosols et les autres constituants minoritaires de l'atmosphère (ozone, SO2, etc) La transmission est l'intégration de l'extinction en fonction de l'altitude, et s'écrit donc ainsi : .. math:: \tau_{i}(\lambda, z) = \int_{z'=z_0}^z \alpha_i(\lambda, z') Soit le signal corrigé de l'altitude :math:`Pz2` : .. math:: Pz2(\lambda_r, z) = \frac{P(\lambda_r, z)}{(z-z_0)^2} On dérive puis on prend le logarithme pour obtenir l'équation suivante : .. math:: \frac{d}{dz}ln \left(Pz2(\lambda_r, z) \right) = \frac{d}{dz}ln \beta(\lambda_e, \lambda_r, z) - \sum_i \left( \alpha_{i}(\lambda_e, z) + \alpha_{i}(\lambda_r, z) \right)) En prenant le signal à deux longueurs d'ondes :math:`\lambda_{on}` et :math:`\lambda_{off}`, on a alors deux équations : .. math:: :label: eq:dial \left\{ \begin{array}{ll} \frac{d}{dz}ln \left(Pz2(\lambda_{on\_r}, z) \right) = \frac{d}{dz}ln \beta(\lambda_{on\_e}, \lambda_{on\_r}, z) - \sum_i \left( \alpha_{i}(\lambda_{on\_e}, z) + \alpha_{i}(\lambda_{on\_r}, z) \right) \\ \frac{d}{dz}ln \left(Pz2(\lambda_{off\_r}, z) \right) = \frac{d}{dz}ln \beta(\lambda_{off\_e}, \lambda_{off\_r}, z) - \sum_i \left( \alpha_{i}(\lambda_{off\_e}, z) + \alpha_{i}(\lambda_{off\_r}, z) \right) \\ \end{array}\right. Afin de simplifier l'équation suivante, nous posons : .. math:: \left\{ \begin{array}{ll} \beta(\lambda_{on\_e}, \lambda_{on\_r}, z) = \beta(\lambda_{on}, z) \\ \beta(\lambda_{off\_e}, \lambda_{off\_r}, z) = \beta(\lambda_{off}, z) \\ \alpha_{i}(\lambda_{on\_e}, z) + \alpha_{i}(\lambda_{on\_r}, z) = \alpha_{i}(\lambda_{on}, z) \\ \alpha_{i}(\lambda_{off\_e}, z) + \alpha_{i}(\lambda_{off\_r}, z) = \alpha_{i}(\lambda_{off}, z) \\ \end{array}\right. On soustrayant les deux équations :eq:`eq:dial`, on obtient alors : .. math:: \frac{d}{dz}ln \left( \frac{Pz2(\lambda_{on}, z)}{Pz2(\lambda_{off}, z)} \right) = \frac{d}{dz}ln \left( \frac{\beta(\lambda_{on}, z)}{\beta(\lambda_{off}, z)}\right) + \sum_i \left( \alpha_i(\lambda_{off}, z) - \alpha_i(\lambda_{on}, z) \right) On peut également écrire :math:`\alpha` de la façon suivante : .. math:: \alpha_i=\sigma_i(\lambda_i, z)*n_i(z) avec : - :math:`\sigma_i` : section efficace du constituant atmosphérique i - :math:`n_i` : concentration du constituant atmosphérique i On peut réécrire une partie de la somme de la façon suivante : .. math:: \alpha_i(\lambda_{off}, z) - \alpha_i(\lambda_{on, z}) = n_i \left( \sigma_{i}(\lambda_{off\_e}, z) + \sigma_{i}(\lambda_{off\_r}, z) - \sigma_{i}(\lambda_{on\_e}, z) - \sigma_{i}(\lambda_{on\_r}, z) \right) Si on pose : .. math:: \Delta \sigma_i (\lambda, z) = \sigma_{i}(\lambda_{off\_e}, z) - \sigma_{i}(\lambda_{off\_r}, z) - \sigma_{i}(\lambda_{on\_e}, z) - \sigma_{i}(\lambda_{on\_r}, z) L'équation devient alors : .. math:: \alpha_i(\lambda_{off}, z) - \alpha_i(\lambda_{on, z}) = n_i(z) \Delta \sigma_i (\lambda, z) L'équation principale devient donc : .. math:: \frac{d}{dz}ln \left( \frac{Pz2(\lambda_{on}, z)}{Pz2(\lambda_{off}, z)} \right) = \frac{d}{dz}ln \left( \frac{\beta(\lambda_{on}, z)}{\beta(\lambda_{off}, z)}\right) + n_{o3}(z) \Delta \sigma_{o3} (\lambda, z) + n_{atm}(z) \Delta \sigma_{atm} (\lambda, z) + \sum_j \left( \alpha_j(\lambda_{off}, z) - \alpha_j(\lambda_{on}, z) \right) D'où la concentration d'ozone peut s'écrire de la façon suivante : .. math:: n_{o3}(z) = \frac{\frac{d}{dz}ln \left( \frac{Pz2(\lambda_{on}, z)}{Pz2(\lambda_{off}, z)} \right) - \frac{d}{dz}ln \left( \frac{\beta(\lambda_{on}, z)}{\beta(\lambda_{off}, z)}\right) - n_{atm}(z) \Delta \sigma_{atm} (\lambda, z) - \sum_j \left( \alpha_j(\lambda_{off}, z) - \alpha_j(\lambda_{on}, z) \right)}{\Delta \sigma_{o3} (\lambda, z) } .. figure:: ./figures/ozone.png :alt: ozone example :align: center Exemple de profil d'ozone Le traitement température +++++++++++++++++++++++++ Fonctionnement de la fonction lidar_actris_v01 ---------------------------------------------- Lancement de la fonction ++++++++++++++++++++++++ La fonction **lidar_ozone_actris_v01** accepte deux arguments : le fichier paramètre et la date à traiter (optionnel). La date est optionnelle, car elle peut se trouver dans le fichier de paramètre également. On le verra plus tard, mais cette fonctionnalité permet de faire tourner le programme avec un fichier paramètre par station. Le fichier paramètre est une chaine de caractère, correspondant au chemin complet du fichier paramètre à utiliser. On peut également mettre en entrée un objet de classe ParametersFromXml (plus d'information concernant cette classe ici : :doc:`xml_parameter_file`). .. code-block:: python def lidar_ozone_actris_v01(path_parameters_file: str | ParametersFromXml | None = None, dt: str | date | None = None) -> Lidar | None: """ Data process for DIAL ozone lidar, ACTRIS. :param path_parameters_file: path of the parameters file (default: None) :type path_parameters_file: str or ParametersFromXml or None :param dt: date to process (default: None) :type dt: str or datetime.date or None :return: the Lidar object :rtype: Lidar or None """ Initialisation ++++++++++++++ Cette partie du code permet de retourner la valeur **None** si la varaible **parameter** en entrée est égale à **None**. Si le programme renvoie un **None** et non une erreure, cela permet de pouvoir relancer le programme de façon automatique, comme on le verra plus tard. .. code-block:: python ######## # Init # ######## if path_parameters_file is None: return None Creation de l'objet Lidar +++++++++++++++++++++++++++++ Tout le programme repose sur la classe **Lidar**. La première étape est donc de l'initialiser. Elle a en entrée les deux paramètres d'entrée de la fonction **lidar_actris_v01**. C'est à partir de cet objet que sera stockées les différentes données (lidars, auxilières, traitements, etc), ainsi que toutes les fonctions de traitement. A noter que les données sont stockées selon le formalisme NetCDF. .. code-block:: python ##################### # 0. Parameter File # ##################### lidar = Lidar(path_parameters_file, dt) Si la date est l'un des paramètres, le programme va enregistrer cette date dans le fichier xml, à l'endroit suivant : .. code-block:: xml Le format iso est le suivant : YYYY-MM-DD. Exemple : 2023-03-28 pour le 28 mars 2023. Importation des données lidars ++++++++++++++++++++++++++++++ Première étape, il faut importer les données lidars dans l'objet **Lidar**. Ce sont ces données qui vont créées les coordonnées du fichier NetCDF : altitude et temps. .. code-block:: python ############################# # 1. Reading Raw lidar Data # ############################# lidar.import_lidar_data() Pour paraméter cette partie, il faut remplir les champs suivants dans le fichier xml : .. code-block:: xml pathname_lidar Teslas 102, 103, 100, 101, 104, 105 Dans ce cas, le programme va chercher à lire le fichier **pathname_lidar**. Il est nécessaire de rajouter le format du fichier (ici, Teslas), ainsi que les signaux types des différentes voies. Les valeurs de ces différents **Signal_Type** sont présenté ici : :doc:`signal_type`. La section **filenames** peut également contenir plusieurs fichiers, s'il y a plusieurs fichiers à lire. Il suffit de séparer les différents noms de fichier par un délimiteur (par defaut : ",", mais peut être remplacé par une option dans le fichier xml, comme ceci : .. code-block:: xml pathname_lidar_1;pathname_lidar_2 Dans le cas d'un fichier dépendent d'une hiérarchie de classement dépendant de la date, le champs **filenames** peut être décomposé en 4 sous-parties : le répertoire de base, un formalisme de répertoire avec la date à l'intérnieur (au format datetime de Python), le prefixe du nom de fichier, le suffixe du nom de fichier, et au milieu, le format de la date dans le nom du fichier. Dans l'exemple suivant, le programme va chercher tous les fichiers dont la structure est la suivante : Configuration_files/Main_files/Lidar/sho240328*.s*. En effet, dans notre cas, la partie date dans le répertoire est vide, donc non utilisé, et le format de la date utilisé dans le fichier est l'année sur 2 digits, le mois sur 2 digits et le jour sur 2 digits. .. code-block:: xml 2024-03-28 Configuration_files/Main_files/Lidar sho *.s* %y%m%d On peut ajouter l'unité de la donnée dans le fichier d'entrée : il y a deux choix possibles, "photocouting" (par défaut) ou "MHz". On le rajoute de la façon suivante : .. code-block:: xml MHz On peut également décider de ne pas sauvegarder toutes les données issues du fichier d'entrée. Pour cela, on va utiliser l'option "Remove_Channel" : il suffit d'ajouter les Signal_Type qu'il ne faudra pas prendre en compte. Exemple : .. code-block:: xml 101 Lorsque le fichier brut contient des données avec des résolutions spatiales différente, il est nécessaire de réaliser une régression afin de mettre toutes les données dans la même résolution. Pour cela, on peut ajouter la ligne suivante. Dans ce cas, on force la résolution à 15 mètres. .. code-block:: xml 15 On peut ajouter un délai au niveau de l'électronique, pour le calcul du vecteur altitude. Pour cela, on utilise le champ **electronic_delay**. Ce délai peut être différent en fonction de la voie. Dans l'exmple suivant, c'est un délai de 100ns pour les 6 premières voies. .. code-block:: xml 100, 100, 100, 100, 100, 100 Importation des données auxilières ++++++++++++++++++++++++++++++++++ Pour fonctionner, le programme a, au minimum, besoin d'un profil de pression et de température, en fonction de l'altitude (composition atmosphérique de l'atmosphère). Ajouter à cela, le programme a besoin pour de l'affichage des composants atmosphériques liées à l'ozone. Cela donne le code suivant : .. code-block:: python ############################ # 2. Import auxiliary data # ############################ lidar.import_atmospheric_component() lidar.import_other_atmospheric_component("ozone") Pour la partie climatologie de l'atmosphère, il faut remplir les champs **Atmosphere_file** du fichier xml. Le programme accepte plusieurs modèles possibles, qui seront joint en fonction de l'altitude. En effet, en régle général, on va utiliser un radiosondage dans les premières couches de l'atmosphère (jusqu'à 30km) qu'on va compléter avec un ou des modèles climatologiques. Il faut donc renseigner les informations de chacun de ces modèles. Pour cela, on va créer des objets **modele_X**, avec X un numéro (ce numéro n'est pas important, mais pour des soucis de lisibilités, il est préréfable de mettre le numéro d'ordre). Pour chaque modèle, il faudra remplir : - le ou les nom des fichiers où se trouvent le ou les fichiers (en utilisant le même modèle que pour le nom du fichier lidar) - le format du ou des fichiers (radiosondage, ncep, arletty) - l'ordre d'importante pour la concaténation - la valeur par défaut - des variables liées au modèle (exemple : le nom de la station pour les fichiers Arletty A noter que, pour le nom des fichiers, une nouvelle variable est apparut : **time_shift**. Il permet de décaler en temps la recherche du fichier. Par exmple, si on souhaite récupérer le radiosondage de la veille, on peut ajouter -1 comme valeur pour ce champ. Voici un exemple de code xml : .. code-block:: xml Configuration_files/Main_files/Atmosphere ht, te .nmc %y%m%d 0 ncep 1 nan OHP Configuration_files/Main_files/Atmosphere/MAP85_PRE.DAT, Configuration_files/Main_files/Atmosphere/MAP85_TEMn.DAT MAP85 1 nan Le programme a besoin des sections efficaces de l'ozone pour fonctionner. Pour les ajouter, il faut remplir le champ **Cross_section_file** du fichier xml : .. code-block:: xml Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_218K.dat, Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_228K.dat, Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_243K.dat, Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_273K.dat, Configuration_files/Main_files/Ozone_cross_section/O3_CRS_BDM_295K.dat Bass And Paur files Le programme a besoin d'un fichier climatologie de l'ozone (pour afficher des courbes). Il faut donc remplir le champ **Climatology_ozone_file** du fichier xml : .. code-block:: xml Configuration_files/Climatology_files/climatology_ozone_ohp.nc netcdf Création d'un signal lidar simulé +++++++++++++++++++++++++++++++++ Pour les figures, et valider les données, on va comparer les données en fonction d'un signal lidar théorique, à chaque étape de traitement. Pour cela, il est nécessaire de le calculer. On utilise alors la commande suivante : .. code-block:: python ########################### # Lidar Signal Simulation # ########################### lidar.create_simulated_lidar() Pré-traitement des données ++++++++++++++++++++++++++ Le programme est fait pour réaliser 3 prétraitements : calcul du fond de ciel, désaturation et correction du facteur géométrique. Pour appliquer ces différents paramètres, il faut utiliser la fonction **apply_process** de la classe Lidar. Cela donne le code suivant : .. code-block:: python ############### # Pre-Process # ############### lidar.background() lidar.saturation() lidar.geometrical_factor() lidar.apply_preprocess() Exemple de programme xml pour le paramètrage de ces fonctions : .. code-block:: xml pelon, pelon, pelon, pelon, pelon, pelon 3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9 77.895,79.84,80.7487,70.1884,60,60 Weighted linear regression, Weighted linear regression, Average, Average, Average, Average 100220,92270,87620,92270,65870,60020 149870,138000,149870,149870,149870,149870 1,1,0,0,0,0 False,False,False,False,False,False mrd Traitement aérosols +++++++++++++++++++ Afin de pouvoir visualiser la présence ou non d'aérosols, et de pouvoir corriger des aérosols le profil ozone, le programme calcule la rétrodiffusion due aux aérosols. .. code-block:: python #################### # Aerosols Process # #################### lidar.aerosols_process() Exemple de programme xml pour le paramètrage de ces fonctions : .. code-block:: xml 40 25 0 Klett Calcul du profil ozone par la méthode DIAL ++++++++++++++++++++++++++++++++++++++++++ Afin de calculer l'ozone, il est nécessaire de filtrer le signal (dans le cas de la méthode DIAL, en le dérivant), de liées les mesures hautes et basses, de calculer l'ozone à partir des signaux dérivés, puis de lier les voies Raymeigh et les voies Raman. .. code-block:: python ################# # Ozone Process # ################# lidar.filtering() lidar.glue(name_param_tag="MergingSlope") lidar.ozone_process() lidar.glue() Exemple de programme xml pour le paramètrage de ces fonctions : .. code-block:: xml Savitzky-Golay dial poly 3,3,4,4,5,5 120,120,120,120,60,60 35,35,30,30,40,40 300,300,300,300,200,200 km df, ir, origin 100, 101 102, 103 linear, linear 0, 0 17.12,17.12 18.02,18.02 132 136 linear 0 9.02 10.07 2 Domaine de validité +++++++++++++++++++ Calcul des colonnes partielles d'ozone ++++++++++++++++++++++++++++++++++++++ Ajout d'un facteur géométrique ++++++++++++++++++++++++++++++ Export données ++++++++++++++ Lancement du programme en bash ------------------------------ Afin de pouvoir lancer le programme avec un script bash, on programme a été crée : run_dial.py (:doc:`run_dial`)