From e260407967ae349231d90452a37fd6a1c0cb03c5 Mon Sep 17 00:00:00 2001 From: aurelien stcherbinine Date: Wed, 13 May 2020 18:19:18 +0200 Subject: [PATCH] update to version 1.4.2 --- README.md | 2 +- docs/doc.md | 22 +- docs/doc_omega_data.md | 125 ++++++++++- docs/doc_omega_plots.md | 134 ++++++++++- docs/doc_useful_functions.md | 19 ++ omegapy/omega_data.py | 420 ++++++++++++++++++++++++++--------- omegapy/omega_plots.py | 291 +++++++++++++++++++++++- omegapy/useful_functions.py | 34 ++- setup.py | 2 +- 9 files changed, 917 insertions(+), 132 deletions(-) diff --git a/README.md b/README.md index 5a6e0f5..717a495 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -![version](https://img.shields.io/badge/version-1.4-blue) +![version](https://img.shields.io/badge/version-1.4.2-blue) ![pythonversion](https://img.shields.io/badge/Python-3.7-blue) ![idlversion](https://img.shields.io/badge/IDL-8.3-blue) diff --git a/docs/doc.md b/docs/doc.md index f00d899..10df02d 100644 --- a/docs/doc.md +++ b/docs/doc.md @@ -21,7 +21,7 @@ Using IDL routines containing in omegapy/omega_routines/*. `import_list_obs_csv(filename)` -`corr_therm(omega)` +`corr_therm(omega, npool=1)` `corr_therm2(omega)` @@ -29,9 +29,9 @@ Using IDL routines containing in omegapy/omega_routines/*. `corr_atm2(omega)` -`corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True)` +`corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True, npool=1)` -`corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True)` +`corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True, npool=1)` `set_omega_bin_path(new_path)` @@ -47,6 +47,12 @@ Using IDL routines containing in omegapy/omega_routines/*. `update_cube_quality(obs_name='ORB*.pkl', folder='auto', version=Version, base_folder=_omega_py_path)` +`test_cube(obs)` + +`compute_list_good_observations(savfilename='liste_good_obs.csv', folder='../data/OMEGA/liste_obs', security=True)` + +`utc_to_my(dt)` + ## [`omegapy.omega_plots`](doc_omega_plots.md) @@ -62,7 +68,13 @@ Display of OMEGAdata cubes. `show_ibd_v2(omega, ibd, cmap='viridis', vmin=None, vmax=None, alpha=None, title='auto', cb_title = 'IBD', lonlim=(None, None), latlim=(None, None), Nfig=None, polar=False, cbar=True)` -`show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', data_list=None, plot=True, grid=True, out=False, **kwargs)` +`show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', data_list=None, mask_list=None, plot=True, grid=True, out=False, negatives_longitudes=False, **kwargs)` + +`save_map_omega_list(omega_list, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, lam=1.085, data_list=None, data_desc='', mask_list=None, sav_filename='auto', base_folder='../data/OMEGA/sav_map_list_v2/', ext='')` + +`load_map_omega_list(filename)` + +`show_omega_list_v2_man(data, grid_lat, grid_lon, infos, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', grid=True, negatives_longitudes=False, **kwargs)` ## [`omegapy.useful_functions`](doc_useful_functions.md) @@ -81,6 +93,8 @@ Some useful generic functions. `load_pickle(filename, disp=True)` +`test_security_overwrite(path)` + `reg_lin(X, Y)` `planck(lam, T)` diff --git a/docs/doc_omega_data.md b/docs/doc_omega_data.md index 12969d1..2fd861c 100644 --- a/docs/doc_omega_data.md +++ b/docs/doc_omega_data.md @@ -19,9 +19,11 @@ Using IDL routines containing in omegapy/omega_routines/*. `load_omega_list(basename, disp=True)` +`load_omega_list2(liste_obs, therm_corr=True, atm_corr=True, **kwargs)` + `import_list_obs_csv(filename)` -`corr_therm(omega)` +`corr_therm(omega, npool=1)` `corr_therm2(omega)` @@ -29,9 +31,9 @@ Using IDL routines containing in omegapy/omega_routines/*. `corr_atm2(omega)` -`corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True)` +`corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True, npool=1)` -`corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True)` +`corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, security=True, overwrite=True, compress=True, npool=1)` `set_omega_bin_path(new_path)` @@ -47,6 +49,12 @@ Using IDL routines containing in omegapy/omega_routines/*. `update_cube_quality(obs_name='ORB*.pkl', folder='auto', version=_Version, base_folder=_omega_py_path)` +`test_cube(obs)` + +`compute_list_good_observations(savfilename='liste_good_obs.csv', folder='../data/OMEGA/liste_obs', security=True)` + +`utc_to_my(dt)` + ### OMEGAdata class ~~~python @@ -335,6 +343,31 @@ omegapy.omega_data.load_omega_list(basename, disp=True): The array of loaded objects of OMEGA/MEx observation. ~~~ +~~~python +omegapy.omega_data.load_omega_list2(liste_obs, therm_corr=True, atm_corr=True, **kwargs) + Load a list of saved OMEGAdata objects, using load_omega(). + + Parameters + ========== + liste_obs : array of str + List of OMEGA/MEx observations ID. + therm_corr : bool or None, optional (default None) + | True -> Only results with thermal correction. + | False -> Only results without thermal correction. + | None -> Both with and without thermal correction. + atm_corr : bool or None, optional (default None) + | True -> Only results with atmospheric correction. + | False -> Only results without atmospheric correction. + | None -> Both with and without atmospheric correction. + **kwargs: + Optional arguments for autoload_omega(). + + Returns + ======= + omega_list : list of OMEGAdata objects + The list of loaded objects of OMEGA/MEx observation. +~~~ + ~~~python omegapy.omega_data.import_list_obs_csv(filename): Import a list of observations ID from a csv file generated by JMars. @@ -354,13 +387,18 @@ omegapy.omega_data.import_list_obs_csv(filename): Correction thermique 1 - Détermitation réflectance à 5µm à partir de celle à 2.4µm puis CN -- Méthode historique (Calvin & Erard 1997) ~~~python -omegapy.omega_data.corr_therm(omega): +omegapy.omega_data.corr_therm(omega, npool=1): Remove the thermal component in the OMEGA hyperspectral cube. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== omega : OMEGAdata The OMEGA observation data. + npool : int, optional (default 1) + Number of parallelized worker process to use. Returns ======= @@ -428,10 +466,13 @@ omegapy.omega_data.corr_atm2(omega): ### Correction & Saving Import an OMEGA/MEx observation, apply a thermal and atmospheric correction (M1) and save the OMEGAdata object at each step. -~~~ +~~~python omegapy.omega_data.corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, - overwrite=True, compress=True): + overwrite=True, compress=True, npool=1): Correction and saving of OMEGA/MEx observations. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== @@ -452,12 +493,17 @@ omegapy.omega_data.corr_save_omega(obsname, folder='auto', base_folder=_omega_py compress : bool, optional (default True) If True, the radiance cube after correction is removed (i.e. set to None) in order to reduce the size of the saved file. + npool : int, optional (default 1) + Number of parallelized worker process to use. ~~~ ~~~ omegapy.omega_data.corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, - security=True, overwrite=True, compress=True): + security=True, overwrite=True, compress=True, npool=1): Correction and saving of a list of OMEGA/MEx observations. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== @@ -478,6 +524,8 @@ omegapy.omega_data.corr_save_omega_list(liste_obs, folder='auto', base_folder=_o compress : bool, optional (default True) If True, the radiance cube after correction is removed (i.e. set to None) in order to reduce the size of the saved file. + npool : int, optional (default 1) + Number of parallelized worker process to use. ~~~ ~~~python @@ -537,7 +585,7 @@ omegapy.omega_data.get_names(omega_list): ~~~python omegapy.omega_data.get_ls(omega_list): - """Return the array of the Solar longitude of each OMEGA/MEx observation in omega_list. + Return the array of the Solar longitude of each OMEGA/MEx observation in omega_list. Parameters ========== @@ -552,7 +600,7 @@ omegapy.omega_data.get_ls(omega_list): ### Update quality -~~python +~~~Python omegapy.omega_data.update_cube_quality(obs_name='ORB*.pkl', folder='auto', version=_Version, base_folder=_omega_py_path): Update the quality attribute of previously saved OMEGAdata objects. @@ -571,3 +619,62 @@ omegapy.omega_data.update_cube_quality(obs_name='ORB*.pkl', folder='auto', versi The base folder path. ~~~ + +### Testing observations quality +~~~python +omegapy.omega_data.test_cube(obs): + Test the quality of an OMEGA/MEx observation from the header informations + witout open it. + + Parameters + ========== + obs : str + The name of the OMEGA observation. + + Returns + ======= + test_quality : bool + | True -> Accepted observation. + | False -> Rejected observation. +~~~ + +~~~python +omegapy.omega_data.compute_list_good_observations(savfilename='liste_good_obs.csv', + folder='../data/OMEGA/liste_obs', security=True): + Scan the available OMEGA/MEx data cubes and list the observations considered as + good quality. + The results are saved in the specified csv file. + + Parameters + ========== + savfilename : str + The name of the csv file to save the data. + folder : str + The name of the folder where the saved file will be located. + Final saved file path = folder + savfilename + security : bool, optional (default True) + Enable / disable checking before overwriting a file. + | True -> Check if the target file already exists before overwriting on it. + And if is the case, you will be asked for a confirmation. + | False -> Didn't care about the already existing files. +~~~ + +### Conversion UTC to MY +~~~python +omegapy.omega_data.utc_to_my(dt): + Convert a UTC datetime to the corresponding Martian Year (MY). + + Martian Years are numbered according to the calendar proposed by R. Todd Clancy + (Clancy et al., Journal of Geophys. Res 105, p 9553, 2000): Martian Year 1 begins + (at a time such that Ls=0) on April 11th, 1955. + + Parameters + ========== + dt : datetime.datetime + The UTC datetime object. + + Returns + ======= + my : int + The corresponding Martian Year. +~~~ diff --git a/docs/doc_omega_plots.md b/docs/doc_omega_plots.md index 783f17e..a64a116 100644 --- a/docs/doc_omega_plots.md +++ b/docs/doc_omega_plots.md @@ -14,7 +14,13 @@ Display of OMEGAdata cubes. `show_ibd_v2(omega, ibd, cmap='viridis', vmin=None, vmax=None, alpha=None, title='auto', cb_title = 'IBD', lonlim=(None, None), latlim=(None, None), Nfig=None, polar=False, cbar=True, grid=True)` -`show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', data_list=None, plot=True, grid=True, out=False, **kwargs)` +`show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', data_list=None, mask_list=None, plot=True, grid=True, out=False, negatives_longitudes=False, **kwargs)` + +`save_map_omega_list(omega_list, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, lam=1.085, data_list=None, data_desc='', mask_list=None, sav_filename='auto', base_folder='../data/OMEGA/sav_map_list_v2/', ext='')` + +`load_map_omega_list(filename)` + +`show_omega_list_v2_man(data, grid_lat, grid_lon, infos, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', grid=True, negatives_longitudes=False, **kwargs)` ### Display cube @@ -241,7 +247,8 @@ omegapy.omega_plots.show_ibd_v2(omega, ibd, cmap='viridis', vmin=None, vmax=None omegapy.omega_plots.show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', - data_list=None, plot=True, grid=True, out=False, **kwargs): + data_list=None, mask_list=None, plot=True, grid=True, out=False, + negatives_longitudes=False, **kwargs): Display an composite map from a list OMEGA/MEx observations, sampled on a new lat/lon grid. Parameters @@ -281,12 +288,21 @@ omegapy.omega_plots.show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_m data_list : 3D array or None, optional (default None) 1D array of the same dimension of `omega_list` containing alternative maps (2D arrays), in the **same order** than the observations of `omega_list`. + mask_list : 3D array + 1D array of the same dimension of `omega_list` containing the masks to remove the + corrupted pixels of each observaiton, in the **same order** than the observations of + `omega_list`. + Each mask is a 2D array, filled with 1 for good pixels and NaN for bad ones. plot : bool, optional (default True) If True -> Diplay the final figure. grid : bool, optional (default True) Enable the display of the lat/lon grid. out : bool, optional (default False) If True -> Return output. + negatives_longitudes : bool, optional (default False) + Argument for non-polar plots. + | True -> longitudes between 0° and 360°. + | False -> longitudus between -180° and 180°. **kwargs: Optional arguments for the plt.pcolormesh() function. @@ -304,3 +320,117 @@ omegapy.omega_plots.show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_m The array indicating which observations have been used to fill each grid position. ~~~ +### Save & restore composite map of several OMEGA observations, sample on a lat/lon grid +~~~python +omegapy.omega_plots.save_map_omega_list(omega_list, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, + pas_lat=0.1, pas_lon=0.1, lam=1.085, data_list=None, data_desc='', + mask_list=None, sav_filename='auto', + base_folder='../data/OMEGA/sav_map_list_v2/', ext=''): + Save the output of the omega_plots.show_omega_list_v2() function with the requested + parameters as a dictionary. + + Parameters + ========== + omega_list : array of OMEGAdata + The list of OMEGA/MEx observations. + lat_min : float, optional (default -90) + The minimal latitude of the grid. + lat_max : float, optional (default 90) + The maximum latitude of the grid. + lon_min : float, optional (default 0) + The minimal longitude of the grid. + lon_max : float, optional (default 360) + The maximal longitude of the grid. + pas_lat : float, optional (default 0.1) + The latitude intervals of the grid. + pas_lon : float, optional (default 0.1) + The longitude intervals of the grid. + lam : float, optional (default 1.085) + The selected wavelength (in µm). + data_list : 3D array or None, optional (default None) + 1D array of the same dimension of `omega_list` containing alternative maps (2D arrays), + in the **same order** than the observations of `omega_list`. + data_desc : str, optional (default '') + Description of the data contained in data_list (if used). + mask_list : 3D array + 1D array of the same dimension of `omega_list` containing the masks to remove the + corrupted pixels of each observaiton, in the **same order** than the observations of + `omega_list`. + Each mask is a 2D array, filled with 1 for good pixels and NaN for bad ones. + sav_filename : str, optional (default 'auto') + The saving file name. + | If 'auto' -> Automatically generated. + base_folder : str, optional (default '../data/OMEGA/sav_map_list_v2/') + The base folder to save the data. + ext : str, optional (default '') + Extension to add at the end of the filename (useful in case of automatic generation). +~~~ + +~~~python +omegapy.omega_plots.load_map_omega_list(filename): + Load and return the result of omega_plots.show_omega_list_v2() previously saved + with save_map_omega_list(). + + Parameters + ========== + filename : str + The file path. + + Returns + ======= + data : 2D array + The omega reflectance at lam, sampled on the new lat/lon grid. + mask : 2D array + The array indicating where the new grid has been filled by the OMEGA data. + grid lat : 2D array + The new latitude grid. + grid lon : 2D array + The new longitude grid. + mask_obs : 2D array of str + The array indicating which observations have been used to fill each grid position. + infos : dict + The informations about the computation of the data. +~~~ + +~~~python +omegapy.omega_plots.show_omega_list_v2_man(data, grid_lat, grid_lon, infos, cmap='Greys_r', vmin=None, vmax=None, + title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', + grid=True, negatives_longitudes=False, **kwargs): + Display an composite map from a list OMEGA/MEx observations, previously sampled on + a new lat/lon grid with show_omega_list_v2() and saved with save_map_omega_list(). + + Parameters + ========== + data : 2D array + The omega reflectance at lam, sampled on the new lat/lon grid. + grid lat : 2D array + The new latitude grid. + grid lon : 2D array + The new longitude grid. + infos : dict + The informations about the computation of the data. + cmap : str, optional (default 'Greys_r') + The matplotlib colormap. + vmin : float or None, optional (default None) + The lower bound of the coloscale. + vmax : float or None, optional (default None) + The upper bound of the colorscale. + title : str, optional (default 'auto') + The title of the figure. + Nfig : int or str or None, optional (default None) + The target figure ID. + polar : bool, optional (default False) + If True -> Use a polar projection for the plot. + cbar : bool, optional (default True) + If True -> Diplay the colorbar. + cb_title : str, optional (default 'auto') + The title of the colorbar. + grid : bool, optional (default True) + Enable the display of the lat/lon grid. + negatives_longitudes : bool, optional (default False) + Argument for non-polar plots. + | True -> longitudes between 0° and 360°. + | False -> longitudus between -180° and 180°. + **kwargs: + Optional arguments for the plt.pcolormesh() function. +~~~ diff --git a/docs/doc_useful_functions.md b/docs/doc_useful_functions.md index 5c2e864..23ee607 100644 --- a/docs/doc_useful_functions.md +++ b/docs/doc_useful_functions.md @@ -16,6 +16,8 @@ Some useful generic functions. `load_pickle(filename, disp=True)` +`test_security_overwrite(path)` + `reg_lin(X, Y)` `planck(lam, T)` @@ -142,6 +144,23 @@ omegapy.useful_functions.load_pickle(filename, disp=True): The loaded object. ~~~ +~~~python +omegapy.useful_functions.test_security_overwrite(path): + Test if a file already exists, and if yes ask the user if he wants to + ovewrite it or not. + + Parameters + ========== + path : str + The target file path. + + Returns + ======= + overwrite : bool + | True -> No existent file, or overwriting allowed. + | False -> Existent file, no overwriting. +~~~ + ### Curve fitting ~~~python omegapy.useful_functions.reg_lin(X, Y): diff --git a/omegapy/omega_data.py b/omegapy/omega_data.py index 64b4aa5..9c5b66a 100644 --- a/omegapy/omega_data.py +++ b/omegapy/omega_data.py @@ -3,7 +3,7 @@ ## omega_data.py ## Created by Aurélien STCHERBININE -## Last modified by Aurélien STCHERBININE : 21/04/2020 +## Last modified by Aurélien STCHERBININE : 13/05/2020 ##---------------------------------------------------------------------------------------- """Importation of OMEGA observations in the OMEGAdata class. @@ -21,11 +21,13 @@ from scipy.optimize import curve_fit from scipy.optimize import minimize from scipy.io import readsav -from datetime import datetime +import datetime import pickle import os import glob import pandas as pd +import multiprocessing as mp +import itertools # Local from . import useful_functions as uf @@ -191,7 +193,7 @@ def __init__(self, obs='', empty=False, data_path=_omega_bin_path): self.emer = np.array([[]]) self.inci = np.array([[]]) self.specmars = np.array([]) - self.utc = datetime.now() + self.utc = datetime.datetime.now() self.orbit = None self.surf_temp = np.array([[]]) self.ic = {'V' : np.arange(265, 333), @@ -281,7 +283,7 @@ def __init__(self, obs='', empty=False, data_path=_omega_bin_path): cube_rf2 = np.moveaxis(cube_rf, [0,1,2], [0,2,1]) # Observation UTC date & time Y, M, D, h, m, s = np.average(utc[:,:6], axis=0).astype(np.int64) - utc_dt = datetime(Y, M, D, h, m, s) + utc_dt = datetime.datetime(Y, M, D, h, m, s) # Longitude pixels grid ny, nx = lon.shape lon_px = np.moveaxis(geocube[:,13:17,:], [0,1,2], [0,2,1]) * 1e-4 @@ -775,7 +777,7 @@ def autosave_omega(omega, folder='auto', base_folder=_omega_py_path, security=Tr target_path = os.path.join(base_folder, folder, savname) # Testing existent file if security: - write = test_security_overwrite(target_path) + write = uf.test_security_overwrite(target_path) else: write = True # Sauvegarde pickle @@ -852,39 +854,44 @@ def autoload_omega(obs_name, folder='auto', version=_Version, base_folder=_omega ##---------------------------------------------------------------------------------------- ## Correction thermique 1 - Détermitation réflectance à 5µm à partir de celle à 2.4µm ## puis CN -- Méthode historique (Calvin & Erard 1997) -def corr_therm_sp(omega, x, y, disp=True): +_omega_tmp = None +def _corr_therm_sp(args): """Remove the thermal component in an OMEGA spectrum, using the historical method based on the reflectance determination using reference spectra from Calvin & Erard (1997). + + Important note : this function is dedicated to internal use in the corr_therm() function, + as it use non-local objects related to the multiprocessing. Parameters ========== - omega : OMEGAdata - The OMEGA observation data. - x : int - The x-coordinate of the pixel. - y : int - The y-coordinate of the pixel. - disp : bool, optional (default True) - If True display the fitted temperature/reflectance in the console. + args : tuple of 3 elements + | x : int + | The x-coordinate of the pixel. + | y : int + | The y-coordinate of the pixel. + | disp : bool, optional (default True) + | If True display the fitted temperature/reflectance in the console. Returns ======= - lam : 1D array - The wavelength array (in µm). sp_rf_corr : 1D array The reflectance spectrum, corrected from the thermal component. T_fit : float The retrieved surface temperature (in K). + x : int + The x-coordinate of the pixel. + y : int + The y-coordinate of the pixel. """ - # Test correction - if omega.therm_corr: - print("\033[1;33mThermal correction already applied\033[0m") - return omega.lam, omega.sp_rf[y, x] + # Extraction arguments + x, y, disp = args # Extraction données + global _omega_tmp + omega = _omega_tmp lam = omega.lam + sp_sol = omega.specmars sp_rf = omega.cube_rf[y, x] sp_i = omega.cube_i[y, x] - sp_sol = omega.specmars ecl = np.cos(omega.inci[y, x] * np.pi/180) # spectels #97-#112 des spectres de ref <-> 2.3-2.5µm fref = os.path.join(package_path, 'OMEGA_dataref/refclair_sombr_omega_CL.dat') # from Erard and Calvin (1997) @@ -924,15 +931,21 @@ def simu_sp_5microns(i_lams, T): sp_rf_corr = deepcopy(sp_rf) sp_rf_corr.fill(np.nan) T_fit = np.nan - return lam, sp_rf_corr, T_fit + # Output + return sp_rf_corr, T_fit, x, y -def corr_therm(omega): +def corr_therm(omega, npool=1): """Remove the thermal component in the OMEGA hyperspectral cube. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== omega : OMEGAdata The OMEGA observation data. + npool : int, optional (default 1) + Number of parallelized worker process to use. Returns ======= @@ -945,19 +958,34 @@ def corr_therm(omega): print("\033[1;33mThermal correction already applied\033[0m") return deepcopy(omega) # Initialisation - omega2 = deepcopy(omega) + global _omega_tmp + _omega_tmp = deepcopy(omega) omega_corr = deepcopy(omega) - ny, nx, nlam = omega2.cube_i.shape - # Correction spectres - for x in tqdm(range(nx)): - for y in tqdm(range(ny)): - sp_rf_corr, surf_temp = corr_therm_sp(omega2, x, y, disp=False)[1:] - omega_corr.cube_rf[y,x] = sp_rf_corr - omega_corr.surf_temp[y,x] = surf_temp - # Sortie + ny, nx, nlam = omega.cube_i.shape + rf_corr = np.zeros((ny,nx,nlam), dtype=np.float64) + surf_temp = np.zeros((ny,nx), dtype=np.float64) + # Itérateur + it = [(x, y, False) for x, y in itertools.product(range(nx), range(ny))] + # Correction thermique + # chunksize = len(it) // npool # Approx size of each process + chunksize = 1 + # pool = mp.Pool(npool) + with mp.Pool(npool) as pool: + for res in tqdm(pool.imap_unordered(_corr_therm_sp, it, chunksize), total=len(it)): + sp_rf_corr, T_fit, x, y = res + rf_corr[y,x] = sp_rf_corr + surf_temp[y,x] = T_fit + pool.close() + _omega_tmp = None + omega_corr.cube_rf = rf_corr + omega_corr.surf_temp = surf_temp + # Update infos omega_corr.therm_corr = True - omega_corr.therm_corr_infos['datetime'] = datetime.now() + omega_corr.therm_corr_infos['datetime'] = datetime.datetime.now() omega_corr.therm_corr_infos['method'] = '(M1) Calvin & Erard' + # Sortie + # tfin = time.time() + # print('Duration : {0:.0f} min {1:.2f} sec'.format((tfin-tini)//60, (tfin-tini)%60)) return omega_corr ##---------------------------------------------------------------------------------------- @@ -1053,7 +1081,7 @@ def corr_therm2(omega): omega_corr.surf_temp[y,x] = surf_temp # Sortie omega_corr.therm_corr = True - omega_corr.therm_corr_infos['datetime'] = datetime.now() + omega_corr.therm_corr_infos['datetime'] = datetime.datetime.now() omega_corr.therm_corr_infos['method'] = '(M2) Simultaneous refl & temp' return omega_corr @@ -1127,7 +1155,7 @@ def corr_atm(omega): omega_corr.cube_rf[y,x] = sp_rf_corr # Sortie omega_corr.atm_corr = True - omega_corr.atm_corr_infos['datetime'] = datetime.now() + omega_corr.atm_corr_infos['datetime'] = datetime.datetime.now() omega_corr.atm_corr_infos['method'] = 'M1 : same reflectance level at 1.93μm and 2.01μm' return omega_corr @@ -1213,7 +1241,7 @@ def corr_atm2(omega): omega_corr.cube_rf[y,x] = cube_rf[y,x] * tr_atm**(-expo) # Sortie omega_corr.atm_corr = True - omega_corr.atm_corr_infos['datetime'] = datetime.now() + omega_corr.atm_corr_infos['datetime'] = datetime.datetime.now() omega_corr.atm_corr_infos['method'] = 'M2 : flattest specra between 1.97µm and 2.00µm' return omega_corr @@ -1273,73 +1301,13 @@ def corr_mode_128(omega): return omega_corr ##---------------------------------------------------------------------------------------- -## Correction cube OMEGA -def corr_omega(omega): - """Remove the thermal and atmospheric component in the OMEGA hyperspectral cube. - - Parameters - ========== - omega : OMEGAdata - The OMEGA observation data. - - Returns - ======= - omega_corr : OMEGAdata - The input OMEGA observation, where the reflectance is corrected from - the thermal and atmospheric component. - """ - # Initialisation - omega2 = deepcopy(omega) - omega_corr = deepcopy(omega) - ny, nx, nlam = omega2.cube_i.shape - # Correction spectres - for x in tqdm(range(nx)): - for y in tqdm(range(ny)): - # Correction thermique - lam, sp_rf_corr = corr_therm_sp(omega2, x, y, disp=False) - # Correction atmosphérique - # sp_rf_corr2 = corr_atm_sp() - omega_corr.cube_rf[y,x] = sp_rf_corr - # Sortie - omega_corr.therm_corr = True - omega_corr.atm_corr = True - return omega_corr - -##---------------------------------------------------------------------------------------- -## Sauvegarde -def test_security_overwrite(path): - """Test if a file already exists, and if yes ask the user if he wants to - ovewrite it or not. - - Parameters - ========== - path : str - The target file path. - - Returns - ======= - overwrite : bool - | True -> No existent file, or overwriting allowed. - | False -> Existent file, no overwriting. - """ - erase = 'n' - if glob.glob(path) != []: - try: - erase = input('Do you really want to erase and replace `' + path + - '` ? (y/N) ') - except KeyboardInterrupt: - erase = 'n' - if erase != 'y' : - print("\033[1mFile preserved\033[0m") - return False - else: - return True - else: - return True - +## Correction & sauvegarde def corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security=True, - overwrite=True, compress=True): + overwrite=True, compress=True, npool=1): """Correction and saving of OMEGA/MEx observations. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== @@ -1360,6 +1328,8 @@ def corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security compress : bool, optional (default True) If True, the radiance cube after correction is removed (i.e. set to None) in order to reduce the size of the saved file. + npool : int, optional (default 1) + Number of parallelized worker process to use. """ if folder == 'auto': folder = 'v' + str(_Version) @@ -1373,11 +1343,11 @@ def corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security else: exists = False if security: - overwrite = test_security_overwrite(basename.format('*')) + overwrite = uf.test_security_overwrite(basename.format('*')) if (not exists) or (exists and overwrite): save_omega(omega, folder=folder, base_folder=base_folder) print('\n\033[01mThermal correction\033[0m') - omega_corr = corr_therm(omega) + omega_corr = corr_therm(omega, npool) if compress: omega_corr.cube_i = None save_omega(omega_corr, folder=folder, base_folder=base_folder, suff='corr_therm') @@ -1388,8 +1358,11 @@ def corr_save_omega(obsname, folder='auto', base_folder=_omega_py_path, security print('\n\033[01;34mExistent files preserved for {0} - v{1}\033[0m\n'.format(name, _Version)) def corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, - security=True, overwrite=True, compress=True): + security=True, overwrite=True, compress=True, npool=1): """Correction and saving of a list of OMEGA/MEx observations. + + Parallelization is implemented using the multiprocessing module. The number of + process to run is controlled by the npool argument. Parameters ========== @@ -1410,13 +1383,15 @@ def corr_save_omega_list(liste_obs, folder='auto', base_folder=_omega_py_path, compress : bool, optional (default True) If True, the radiance cube after correction is removed (i.e. set to None) in order to reduce the size of the saved file. + npool : int, optional (default 1) + Number of parallelized worker process to use. """ N = len(liste_obs) if folder == 'auto': folder = 'v' + str(_Version) for i, obsname in enumerate(liste_obs): print('\n\033[01mComputing observation {0} / {1} : {2}\033[0m\n'.format(i+1, N, obsname)) - corr_save_omega(obsname, folder, base_folder, security, overwrite, compress) + corr_save_omega(obsname, folder, base_folder, security, overwrite, compress, npool) print("\n\033[01;32m Done\033[0m\n") ##---------------------------------------------------------------------------------------- @@ -1578,6 +1553,241 @@ def update_cube_quality(obs_name='ORB*.pkl', folder='auto', version=_Version, save_omega(omega, fname, '', '', '', '', False) print('\033[1m{0} files updated\033[0m'.format(len(fnames))) +##---------------------------------------------------------------------------------------- +## Importation liste OMEGAdata avec filtrage automatisé +def load_omega_list2(liste_obs, therm_corr=True, atm_corr=True, **kwargs): + """Load a list of saved OMEGAdata objects, using load_omega(). + + Parameters + ========== + liste_obs : array of str + List of OMEGA/MEx observations ID. + therm_corr : bool or None, optional (default None) + | True -> Only results with thermal correction. + | False -> Only results without thermal correction. + | None -> Both with and without thermal correction. + atm_corr : bool or None, optional (default None) + | True -> Only results with atmospheric correction. + | False -> Only results without atmospheric correction. + | None -> Both with and without atmospheric correction. + **kwargs: + Optional arguments for autoload_omega(). + + Returns + ======= + omega_list : list of OMEGAdata objects + The list of loaded objects of OMEGA/MEx observation. + """ + omega_list = [] + Nabs = 0 + OBC = readsav('../data/OMEGA_dataref/OBC_OMEGA_OCT2017.sav') + good_orbits_OBC = np.array(OBC['good_orbits'][0], dtype=int) + for i, obsname in enumerate(tqdm(liste_obs)): + omega = autoload_omega(obsname, therm_corr=therm_corr, atm_corr=atm_corr, disp=False, + **kwargs) + if omega is None: + Nabs += 1 + continue + if not omega.orbit in good_orbits_OBC: + continue + if omega.quality == 0: + continue + if omega.target != 'MARS': + continue + if omega.mode_channel != 1: + continue + if omega.data_quality == 0: + continue + if omega.point_mode == 'N/A': + continue + if (int(omega.name[-1]) == 0) and (omega.npixel == 64) and (omega.bits_per_data == 1): + continue + # if omega.npixel == 16: + # continue + omega_list.append(omega) + Ntot = len(liste_obs) + Nacc = len(omega_list) + Nrej = Ntot - Nacc - Nabs + print('\n\033[1m{0} observations in list_obs\n'.format(Ntot) + + '{0} loaded, {1} rejected, {2} not found\033[0m\n'.format(Nacc, Nrej, Nabs)) + return omega_list + +def test_cube(obs): + """Test the quality of an OMEGA/MEx observation from the header informations + witout open it. + + Parameters + ========== + obs : str + The name of the OMEGA observation. + + Returns + ======= + test_quality : bool + | True -> Accepted observation. + | False -> Rejected observation. + """ + # Recherhe nom de fichier + data_path = _omega_bin_path + obs_name = uf.myglob(os.path.join(data_path, '*' + obs + '*.QUB')) + if obs_name is None: + print("\033[1;33mAborted\033[0m") + return False + nomfic0 = obs_name[obs_name.rfind('/')+1:-4] # Récupération nom + numCube = int(nomfic0[-1]) + # Lecture header fichier .QUB + hd_qub = _read_header(obs_name[:-4] + '.QUB') + summation = np.int64(hd_qub['DOWNTRACK_SUMMING']) + bits_per_data = np.float64(hd_qub['INST_CMPRS_RATE']) + data_quality = np.int64(hd_qub['DATA_QUALITY_ID']) + mode_channel_tmp = hd_qub['COMMAND_DESC'][34:36] + if mode_channel_tmp == 'EF': + mode_channel = 1 + elif mode_channel_tmp == '80': + mode_channel = 2 + elif mode_channel_tmp == 'C7': + mode_channel = 3 + else: + mode_channel = mode_channel_tmp + # Lecture header fichier .NAV + if glob.glob(obs_name[:-4] + '.NAV') == []: + return False # Pas de fichier .NAV + hd_nav = _read_header(obs_name[:-4] + '.NAV') + npixel, npara, nscan = np.array(hd_nav['CORE_ITEMS'][1:-1].split(','), dtype=np.int64) + point_mode = hd_nav['SPACECRAFT_POINTING_MODE'][1:-1] + target = hd_nav['TARGET_NAME'] + # Test si cube OK + if target != 'MARS': + return False + elif mode_channel != 1: + return False + elif data_quality == 0: + return False + elif point_mode == 'N/A': + return False + elif (numCube == 0) and (npixel == 64) and (bits_per_data == 1): + return False + else: + return True + +##---------------------------------------------------------------------------------------- +## List available good observations & generate csv file +def compute_list_good_observations(savfilename='liste_good_obs.csv', + folder='../data/OMEGA/liste_obs', security=True): + """Scan the available OMEGA/MEx data cubes and list the observations considered as + good quality. + The results are saved in the specified csv file. + + Parameters + ========== + savfilename : str + The name of the csv file to save the data. + folder : str + The name of the folder where the saved file will be located. + Final saved file path = folder + savfilename + security : bool, optional (default True) + Enable / disable checking before overwriting a file. + | True -> Check if the target file already exists before overwriting on it. + And if is the case, you will be asked for a confirmation. + | False -> Didn't care about the already existing files. + """ + # Test existence fichier de sauvegarde + sav_file_path = os.path.join(folder, savfilename) + if security: + test_overwrite = uf.test_security_overwrite(sav_file_path) + if not test_overwrite: + return None + # Liste observations disponibles + bin_obs_list = glob.glob(os.path.join(_omega_bin_path, 'ORB*.QUB')) + bin_obs_list.sort() + # Initialisation + gobs = open(sav_file_path, 'w', encoding='utf-8') + gobs.write('obsname, Ls [°], lat_min [°], lat_max [°], lon_min [°], lon_max [°], ' + + 'UTC date/time, Npixel, Nscan\n') + Nacc = 0 + # Test qualité de chaque observation + for obs_name in tqdm(bin_obs_list): + nomfic0 = obs_name[obs_name.rfind('/')+1:-4] # Récupération nom + numCube = nomfic0[-1] + # Lecture header fichier .QUB + hd_qub = _read_header(obs_name[:-4] + '.QUB') + summation = np.int64(hd_qub['DOWNTRACK_SUMMING']) + bits_per_data = np.float64(hd_qub['INST_CMPRS_RATE']) + data_quality = np.int64(hd_qub['DATA_QUALITY_ID']) + mode_channel_tmp = hd_qub['COMMAND_DESC'][34:36] + if mode_channel_tmp == 'EF': + mode_channel = 1 + elif mode_channel_tmp == '80': + mode_channel = 2 + elif mode_channel_tmp == 'C7': + mode_channel = 3 + else: + mode_channel = mode_channel_tmp + # Lecture header fichier .NAV + if glob.glob(obs_name[:-4] + '.NAV') == []: + continue + hd_nav = _read_header(obs_name[:-4] + '.NAV') + npixel, npara, nscan = np.array(hd_nav['CORE_ITEMS'][1:-1].split(','), dtype=np.int64) + point_mode = hd_nav['SPACECRAFT_POINTING_MODE'][1:-1] + target = hd_nav['TARGET_NAME'] + test = True + # Test si cube OK + if target != 'MARS': + continue + elif mode_channel != 1: + continue + elif data_quality == 0: + continue + elif point_mode == 'N/A': + continue + elif (numCube == '0') and (npixel == 64) and (bits_per_data == 1): + continue + else: + # Si OK -> sauvegarde des infos dans le fichier + gobs.write(('{obsname:s}, {ls:s}, {lat_min:s}, {lat_max:s}, {lon_min:s},' + +'{lon_max:s}, {utc:s}, {npixel:d}, {nscan:d}\n').format(obsname = nomfic0, + ls = hd_nav['SOLAR_LONGITUDE'], + lat_min = hd_nav['MINIMUM_LATITUDE'], + lat_max = hd_nav['MAXIMUM_LATITUDE'], + lon_min = hd_nav['WESTERNMOST_LONGITUDE'], + lon_max = hd_nav['EASTERNMOST_LONGITUDE'], + utc = hd_nav['START_TIME'][:16], + npixel = npixel, + nscan = nscan)) + Nacc += 1 + gobs.close() + # Résultats + Ntot = len(bin_obs_list) + Nrej = Ntot - Nacc + print('\n\033[1m{0} observations found\n'.format(Ntot) + + '{0} accepted, {1} rejected\033[0m'.format(Nacc, Nrej)) + print('\n\033[01;34mResults saved in \033[0;03m' + sav_file_path + '\033[0m') + +##---------------------------------------------------------------------------------------- +## UTC to MY +def utc_to_my(dt): + """Convert a UTC datetime to the corresponding Martian Year (MY). + + Martian Years are numbered according to the calendar proposed by R. Todd Clancy + (Clancy et al., Journal of Geophys. Res 105, p 9553, 2000): Martian Year 1 begins + (at a time such that Ls=0) on April 11th, 1955. + + Parameters + ========== + dt : datetime.datetime + The UTC datetime object. + + Returns + ======= + my : int + The corresponding Martian Year. + """ + datetime_my1 = datetime.datetime(1955, 4, 11) # Start MY 1 + my_sol_duration = 668.6 # Nb of Martian sols during a MY + sol_sec_duration = 88775.245 # Duration of a sol in seconds + my = int( (dt - datetime_my1).total_seconds() // (my_sol_duration * sol_sec_duration)) + 1 + return my + ##---------------------------------------------------------------------------------------- ## End of code ##---------------------------------------------------------------------------------------- diff --git a/omegapy/omega_plots.py b/omegapy/omega_plots.py index 7fffe6f..3a96740 100644 --- a/omegapy/omega_plots.py +++ b/omegapy/omega_plots.py @@ -3,7 +3,7 @@ ## omega_plots.py ## Created by Aurélien STCHERBININE -## Last modified by Aurélien STCHERBININE : 03/04/2020 +## Last modified by Aurélien STCHERBININE : 12/05/2020 ##---------------------------------------------------------------------------------------- """Display of OMEGAdata cubes. @@ -18,6 +18,8 @@ from matplotlib.widgets import Slider, Button, RadioButtons from copy import deepcopy from tqdm import tqdm +import datetime +import os # Local from . import useful_functions as uf from . import omega_data as od @@ -33,7 +35,7 @@ mpl.rcParams['keymap.forward'].remove('right') # Name of the current file -py_file = 'omega_plots.py' +_py_file = 'omega_plots.py' ##---------------------------------------------------------------------------------------- ## Affichage cube @@ -826,10 +828,33 @@ def check_list_data_omega(omega_list, data_list, disp=True): if disp: print("\033[01;32mCompatibility between omega_list and data_list OK\033[0m") +def check_list_mask_omega(omega_list, mask_list, disp=True): + """Check the compatibility between mask_list and the list of OMEGA/MEx observations. + Raise ValueError if uncompatibility. + + Parameters + ========== + omega_list : array of OMEGAdata + List of OMEGA/MEx observations. + mask_list : 3D array + List of masks to remove the corrupted pixels of each OMEGA/MEx observation. + disp : bool, optional (default True) + Enable the display of the result of the test. + """ + if len(omega_list) != len(mask_list): + raise ValueError("omega_list and mask_list must have the same size") + else: + for i in range(len(omega_list)): + if omega_list[i].lat.shape != mask_list[i].shape: + raise ValueError("The shapes of items {0} of omega_list and mask_list does not match.".format(i)) + if disp: + print("\033[01;32mCompatibility between omega_list and mask_list OK\033[0m") + def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, pas_lat=0.1, pas_lon=0.1, cmap='Greys_r', vmin=None, vmax=None, title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', - data_list=None, plot=True, grid=True, out=False, **kwargs): + data_list=None, mask_list=None, plot=True, grid=True, out=False, + negatives_longitudes=False, **kwargs): """Display an composite map from a list OMEGA/MEx observations, sampled on a new lat/lon grid. Parameters @@ -869,12 +894,21 @@ def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0 data_list : 3D array or None, optional (default None) 1D array of the same dimension of `omega_list` containing alternative maps (2D arrays), in the **same order** than the observations of `omega_list`. + mask_list : 3D array + 1D array of the same dimension of `omega_list` containing the masks to remove the + corrupted pixels of each observaiton, in the **same order** than the observations of + `omega_list`. + Each mask is a 2D array, filled with 1 for good pixels and NaN for bad ones. plot : bool, optional (default True) If True -> Diplay the final figure. grid : bool, optional (default True) Enable the display of the lat/lon grid. out : bool, optional (default False) If True -> Return output. + negatives_longitudes : bool, optional (default False) + Argument for non-polar plots. + | True -> longitudes between 0° and 360°. + | False -> longitudus between -180° and 180°. **kwargs: Optional arguments for the plt.pcolormesh() function. @@ -899,10 +933,16 @@ def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0 data, mask = np.zeros((Nlon, Nlat)), np.zeros((Nlon, Nlat)) mask_obs = np.ndarray((Nlon, Nlat), dtype=object) mask_obs.fill('') + if not (mask_list is None): + check_list_mask_omega(omega_list, mask_list, disp=True) if data_list is None: - for omega in tqdm(omega_list): + for i, omega in enumerate(tqdm(omega_list)): i_lam = uf.where_closer(lam, omega.lam) - data0, mask0 = proj_grid(omega, omega.cube_rf[:,:,i_lam], lat_min, lat_max, + if mask_list is None: + data_tmp = omega.cube_rf[:,:,i_lam] # Reflectance without mask + else: + data_tmp = omega.cube_rf[:,:,i_lam] * mask_list[i] # Reflectance with mask + data0, mask0 = proj_grid(omega, data_tmp, lat_min, lat_max, lon_min, lon_max, pas_lat, pas_lon)[:2] data += np.nan_to_num(data0) # Conversion NaN -> 0 pour somme des images mask += mask0 @@ -910,7 +950,11 @@ def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0 else: check_list_data_omega(omega_list, data_list, disp=True) for i, omega in enumerate(tqdm(omega_list)): - data0, mask0 = proj_grid(omega, data_list[i], lat_min, lat_max, + if mask_list is None: + data_tmp = data_list[i] # Data without mask + else: + data_tmp = data_list[i] * mask_list[i] # Data with mask + data0, mask0 = proj_grid(omega, data_tmp, lat_min, lat_max, lon_min, lon_max, pas_lat, pas_lon)[:2] data += np.nan_to_num(data0) # Conversion NaN -> 0 pour somme des images mask += mask0 @@ -930,15 +974,29 @@ def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0 ax.set_yticklabels([]) # remove the latitude values in the plot ax.set_theta_offset(-np.pi/2) # longitude origin at the bottom plt.xlim(0, 2*np.pi) - if (lat_max == 90) & (lat_min >= 0): + if np.abs(lat_max) >= np.abs(lat_min): latlim = (lat_max, lat_min) else: latlim = (lat_min, lat_max) plt.ylim(latlim) else: - plt.pcolormesh(grid_lon, grid_lat, data2, cmap=cmap, vmin=vmin, - vmax=vmax, **kwargs) + if negatives_longitudes and (lon_max > 180): + n_neg_lon = np.sum(grid_lon[:,0] > 180) # nb of negative longitudes (>180°) + i_lon180 = np.where(grid_lon[:,0] > 180)[0][0] # 1st index of lon > 180° + grid_lon_nl = deepcopy(grid_lon) # new longitude grid [-180°, 180°] + grid_lon_nl[:n_neg_lon] = grid_lon[i_lon180-1:-1] - 360 + grid_lon_nl[n_neg_lon:] = grid_lon[:i_lon180] + data2_nl = deepcopy(data2) # new data array + data2_nl[:n_neg_lon] = data2[i_lon180-1:] + data2_nl[n_neg_lon:] = data2[:i_lon180-1] + plt.pcolormesh(grid_lon_nl, grid_lat, data2_nl, cmap=cmap, vmin=vmin, + vmax=vmax, **kwargs) + lon_min, lon_max = grid_lon_nl[[0,-1], 0] # new longitude bounds + else: + plt.pcolormesh(grid_lon, grid_lat, data2, cmap=cmap, vmin=vmin, + vmax=vmax, **kwargs) plt.gca().axis('equal') + plt.gca().set_adjustable('box') plt.xlabel('Longitude [°]') plt.ylabel('Latitude [°]') plt.xlim(lon_min, lon_max) @@ -971,6 +1029,221 @@ def show_omega_list_v2(omega_list, lam=1.085, lat_min=-90, lat_max=90, lon_min=0 mask2 = np.clip(mask, 0, 1) return data2, mask2, grid_lat, grid_lon, mask_obs +##---------------------------------------------------------------------------------------- +## Sauvegarde résultats +def save_map_omega_list(omega_list, lat_min=-90, lat_max=90, lon_min=0, lon_max=360, + pas_lat=0.1, pas_lon=0.1, lam=1.085, data_list=None, data_desc='', + mask_list=None, sav_filename='auto', + base_folder='../data/OMEGA/sav_map_list_v2/', ext=''): + """Save the output of the omega_plots.show_omega_list_v2() function with the requested + parameters as a dictionary. + + Parameters + ========== + omega_list : array of OMEGAdata + The list of OMEGA/MEx observations. + lat_min : float, optional (default -90) + The minimal latitude of the grid. + lat_max : float, optional (default 90) + The maximum latitude of the grid. + lon_min : float, optional (default 0) + The minimal longitude of the grid. + lon_max : float, optional (default 360) + The maximal longitude of the grid. + pas_lat : float, optional (default 0.1) + The latitude intervals of the grid. + pas_lon : float, optional (default 0.1) + The longitude intervals of the grid. + lam : float, optional (default 1.085) + The selected wavelength (in µm). + data_list : 3D array or None, optional (default None) + 1D array of the same dimension of `omega_list` containing alternative maps (2D arrays), + in the **same order** than the observations of `omega_list`. + data_desc : str, optional (default '') + Description of the data contained in data_list (if used). + mask_list : 3D array + 1D array of the same dimension of `omega_list` containing the masks to remove the + corrupted pixels of each observaiton, in the **same order** than the observations of + `omega_list`. + Each mask is a 2D array, filled with 1 for good pixels and NaN for bad ones. + sav_filename : str, optional (default 'auto') + The saving file name. + | If 'auto' -> Automatically generated. + base_folder : str, optional (default '../data/OMEGA/sav_map_list_v2/') + The base folder to save the data. + ext : str, optional (default '') + Extension to add at the end of the filename (useful in case of automatic generation). + """ + # Initialization filename + if sav_filename == 'auto': + sav_filename = ('res_show_omega_list_v2__lat{0:0>2d}-{1:0>2d}_pas{2:0}_' + + 'lon{3:0>3d}-{4:0>3d}_pas{5:0}__{6:s}.pkl').format( + lat_min, lat_max, pas_lat, lon_min, lon_max, pas_lon, ext) + sav_filename2 = os.path.join(base_folder, sav_filename) + if data_list is None: + data_desc = 'Reflectance at λ = {0:0} µm'.format(lam) + elif data_desc == '': + data_desc = 'unknown input data' + # Compute the data sampling + data, mask, grid_lat, grid_lon, mask_obs = show_omega_list_v2(omega_list, + lam, lat_min, lat_max, lon_min, lon_max, pas_lat, pas_lon, + data_list=data_list, mask_list=mask_list, plot=False, out=True) + # Sav file + input_params = { + 'omega_list' : od.get_names(omega_list), + 'lat_min' : lat_min, + 'lat_max' : lat_max, + 'lon_min' : lon_min, + 'lon_max' : lon_max, + 'pas_lat' : pas_lat, + 'pas_lon' : pas_lon, + 'data' : data_desc, + 'filename': sav_filename, + 'datetime': datetime.datetime.now().strftime('%d/%m/%Y %H:%M')} + save_file = { + 'data' : data, + 'mask' : mask, + 'grid_lat' : grid_lat, + 'grid_lon' : grid_lon, + 'mask_obs' : mask_obs, + 'infos' : input_params + } + uf.save_pickle(save_file, sav_filename2, True) + +def load_map_omega_list(filename): + """Load and return the result of omega_plots.show_omega_list_v2() previously saved + with save_map_omega_list(). + + Parameters + ========== + filename : str + The file path. + + Returns + ======= + data : 2D array + The omega reflectance at lam, sampled on the new lat/lon grid. + mask : 2D array + The array indicating where the new grid has been filled by the OMEGA data. + grid lat : 2D array + The new latitude grid. + grid lon : 2D array + The new longitude grid. + mask_obs : 2D array of str + The array indicating which observations have been used to fill each grid position. + infos : dict + The informations about the computation of the data. + """ + loaded_dict = uf.load_pickle(filename, True) + data, mask, grid_lat, grid_lon, mask_obs, infos = loaded_dict.values() + return data, mask, grid_lat, grid_lon, mask_obs, infos + +def show_omega_list_v2_man(data, grid_lat, grid_lon, infos, cmap='Greys_r', vmin=None, vmax=None, + title='auto', Nfig=None, polar=False, cbar=True, cb_title='auto', + grid=True, negatives_longitudes=False, **kwargs): + """Display an composite map from a list OMEGA/MEx observations, previously sampled on + a new lat/lon grid with show_omega_list_v2() and saved with save_map_omega_list(). + + Parameters + ========== + data : 2D array + The omega reflectance at lam, sampled on the new lat/lon grid. + grid lat : 2D array + The new latitude grid. + grid lon : 2D array + The new longitude grid. + infos : dict + The informations about the computation of the data. + cmap : str, optional (default 'Greys_r') + The matplotlib colormap. + vmin : float or None, optional (default None) + The lower bound of the coloscale. + vmax : float or None, optional (default None) + The upper bound of the colorscale. + title : str, optional (default 'auto') + The title of the figure. + Nfig : int or str or None, optional (default None) + The target figure ID. + polar : bool, optional (default False) + If True -> Use a polar projection for the plot. + cbar : bool, optional (default True) + If True -> Diplay the colorbar. + cb_title : str, optional (default 'auto') + The title of the colorbar. + grid : bool, optional (default True) + Enable the display of the lat/lon grid. + negatives_longitudes : bool, optional (default False) + Argument for non-polar plots. + | True -> longitudes between 0° and 360°. + | False -> longitudus between -180° and 180°. + **kwargs: + Optional arguments for the plt.pcolormesh() function. + """ + lat_min, lat_max = infos['lat_min'], infos['lat_max'] + lon_min, lon_max = infos['lon_min'], infos['lon_max'] + if title == 'auto': + title = 'Composite map from OMEGA/MEx observations' + fig = plt.figure(Nfig) + Nfig = fig.number # get the actual figure number if Nfig=None + if polar: + ax = plt.axes(polar=True) + plt.pcolormesh(grid_lon*np.pi/180, grid_lat, data, cmap=cmap, + vmin=vmin, vmax=vmax, **kwargs) + ax.set_yticklabels([]) # remove the latitude values in the plot + ax.set_theta_offset(-np.pi/2) # longitude origin at the bottom + plt.xlim(0, 2*np.pi) + # if (lat_max == 90) & (lat_min >= 0): + if np.abs(lat_max) >= np.abs(lat_min): + latlim = (lat_max, lat_min) + else: + latlim = (lat_min, lat_max) + plt.ylim(latlim) + else: + if negatives_longitudes and (lon_max > 180): + n_neg_lon = np.sum(grid_lon[:,0] > 180) # nb of negative longitudes (>180°) + i_lon180 = np.where(grid_lon[:,0] > 180)[0][0] # 1st index of lon > 180° + grid_lon_nl = deepcopy(grid_lon) # new longitude grid [-180°, 180°] + grid_lon_nl[:n_neg_lon] = grid_lon[i_lon180-1:-1] - 360 + grid_lon_nl[n_neg_lon:] = grid_lon[:i_lon180] + data_nl = deepcopy(data) # new data array + data_nl[:n_neg_lon] = data[i_lon180-1:] + data_nl[n_neg_lon:] = data[:i_lon180-1] + plt.pcolormesh(grid_lon_nl, grid_lat, data_nl, cmap=cmap, vmin=vmin, + vmax=vmax, **kwargs) + lon_min, lon_max = grid_lon_nl[[0,-1], 0] # new longitude bounds + else: + plt.pcolormesh(grid_lon, grid_lat, data, cmap=cmap, vmin=vmin, + vmax=vmax, **kwargs) + plt.gca().axis('equal') + plt.gca().set_adjustable('box') + plt.xlabel('Longitude [°]') + plt.ylabel('Latitude [°]') + plt.xlim(lon_min, lon_max) + plt.ylim(lat_min, lat_max) + if cbar: + if cb_title == 'auto': + cb_title = infos['data'] + cb = plt.colorbar() + cb.set_label(cb_title) + if grid: + ax = plt.figure(Nfig).get_axes()[0] + lonlim = ax.get_xlim() + latlim = ax.get_ylim() + lon_sgn = np.sign(lonlim[1] - lonlim[0]) + lat_sgn = np.sign(latlim[1] - latlim[0]) + lon_grid = np.arange(np.round(lonlim[0]/10)*10, np.round(lonlim[1]/10)*10+lon_sgn, + 10 * lon_sgn) # 10° grid in longitude + lat_grid = np.arange(np.round(latlim[0]/10)*10, np.round(latlim[1]/10)*10+lat_sgn, + 10 * lat_sgn) # 10° grid in latitude + plt.grid() + if polar: + ax.set_rticks(lat_grid) + else: + ax.set_xticks(lon_grid) + ax.set_yticks(lat_grid) + plt.title(title) + plt.tight_layout() + ##---------------------------------------------------------------------------------------- ## End of code ##---------------------------------------------------------------------------------------- diff --git a/omegapy/useful_functions.py b/omegapy/useful_functions.py index f27b320..877b049 100644 --- a/omegapy/useful_functions.py +++ b/omegapy/useful_functions.py @@ -3,7 +3,7 @@ ## useful_functions.py ## Created by Aurélien STCHERBININE -## Last modified by Aurélien STCHERBININE : 21/04/2020 +## Last modified by Aurélien STCHERBININE : 24/04/2020 ##----------------------------------------------------------------------------------- """Useful generics functions. @@ -388,6 +388,38 @@ def load_pickle(filename, disp=True): print('\033[03m' + filename2 + '\033[0;01;34m loaded\033[0m') return obj +##----------------------------------------------------------------------------------- +## Test existence avant sauvegarde +def test_security_overwrite(path): + """Test if a file already exists, and if yes ask the user if he wants to + ovewrite it or not. + + Parameters + ========== + path : str + The target file path. + + Returns + ======= + overwrite : bool + | True -> No existent file, or overwriting allowed. + | False -> Existent file, no overwriting. + """ + erase = 'n' + if glob.glob(path) != []: + try: + erase = input('Do you really want to erase and replace \033[3m' + path + + '\033[0m ? (y/N) ') + except KeyboardInterrupt: + erase = 'n' + if erase != 'y' : + print("\033[1mFile preserved\033[0m") + return False + else: + return True + else: + return True + ##----------------------------------------------------------------------------------- ## End of code ##----------------------------------------------------------------------------------- diff --git a/setup.py b/setup.py index c759552..a704ad8 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ setuptools.setup( name='omegapy', - version='1.4.1', + version='1.4.2', author='Aurélien Stcherbinine', author_email='aurelien.stcherbinine@ias.u-psud.fr', description='Python tools for OMEGA/MEx observations analysis',