Index: trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3792) +++ trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3809) @@ -1,7 +1,7 @@ #!/usr/bin/env python3 ####################################################################################################### -### Python script to output the stratification data over time from the "restartpem#.nc" files files ### +### Python script to output stratification data over time from "restartpem#.nc" files ### +### and to plot orbital parameters from "obl_ecc_lsp.asc" ### ####################################################################################################### - import os @@ -47,5 +47,15 @@ ).strip() or "info_PEM.txt" - return folder_path, base_name, infofile + orbfile = input( + "Enter the name of the orbital parameters ASCII file " + "(press Enter for default [obl_ecc_lsp.asc]): " + ).strip() or "obl_ecc_lsp.asc" + while not os.path.isfile(orbfile): + print(f" » \"{orbfile}\" does not exist or is not a file.") + orbfile = input( + "Enter a valid orbital parameters ASCII filename (press Enter for default [obl_ecc_lsp.asc]): " + ).strip() or "info_PEM.txt" + + return folder_path, base_name, infofile, orbfile @@ -195,5 +205,5 @@ """ Given raw_prop_arrays for 'co2_ice', 'h2o_ice', 'dust', 'pore' (in meters), - normalize each set of strata so that the sum of those four = 1 per strata. + normalize each set of strata so that the sum of those four = 1 per cell. Returns: - frac_arrays: dict mapping same keys -> 4D arrays of fractions (0..1). @@ -228,11 +238,14 @@ """ Reads "info_PEM.txt". Expects: - - First line: parameters (ignored). - - Each subsequent line: floats where first value is timestamp. - Returns: 1D numpy array of timestamps. + - First line: parameters where the 3rd value is martian_to_earth conversion factor. + - Each subsequent line: floats where first value is simulation timestamp (in Mars years). + Returns: + - date_time: 1D numpy array of timestamps (Mars years) + - martian_to_earth: float conversion factor """ date_time = [] with open(file_name, 'r') as fp: - fp.readline() + first = fp.readline().split() + martian_to_earth = float(first[2]) for line in fp: parts = line.strip().split() @@ -243,5 +256,5 @@ except ValueError: continue - return np.array(date_time, dtype=np.float64) + return np.array(date_time, dtype=np.float64), martian_to_earth @@ -293,40 +306,23 @@ ): """ - Build a reference grid and interpolate strata fractions (0..1) onto it. - - Also returns a 'top_index' array of shape (ngrid, ntime, nslope) that - indicates, for each (ig, t_idx, isl), the number of ref_grid levels - covered by the topmost valid stratum. - - Args: - - heights_data: list of lists where heights_data[t][isl] is a 2D array - (ngrid, n_strata_current) of top_elevation values. - - prop_arrays: dict mapping each property_name to a 4D array of shape - (ngrid, ntime, nslope, max_nb_str) holding fractions [0..1]. - - min_base_for_interp: float; if exclude_sub=True, this is 0.0. - - max_top_elev: float - - dz: float - - exclude_sub: bool. If True, ignore strata with elevation < 0. + Build a reference elevation grid and interpolate strata fractions onto it. Returns: - ref_grid: 1D array of elevations (nz,) - - gridded_data: dict mapping each property_name to a 4D array of shape - (ngrid, ntime, nslope, nz) with interpolated fractions. - - top_index: 3D array (ngrid, ntime, nslope) of ints: number of levels - of ref_grid covered by the topmost stratum. - """ - # Build ref_grid, ensuring at least two points if surface-only and dz > max_top_elev + - gridded_data: dict mapping each property_name to 4D array + (ngrid, ntime, nslope, nz) with interpolated fractions. + - top_index: 3D array (ngrid, ntime, nslope) of ints: + number of levels covered by the topmost stratum. + """ if exclude_sub and (dz > max_top_elev): ref_grid = np.array([0.0, max_top_elev], dtype=np.float32) else: - ref_grid = np.arange(min_base_for_interp, max_top_elev + dz / 2, dz) + ref_grid = np.arange(min_base_for_interp, max_top_elev + dz/2, dz) nz = len(ref_grid) print(f"> Number of reference grid points = {nz}") - # Dimensions sample_prop = next(iter(prop_arrays.values())) - ngrid, ntime, nslope, max_nb_str = sample_prop.shape[0:4] - - # Prepare outputs + ngrid, ntime, nslope, max_nb_str = sample_prop.shape + gridded_data = { prop: np.full((ngrid, ntime, nslope, nz), -1.0, dtype=np.float32) @@ -342,6 +338,5 @@ continue - raw_h = h_mat[ig, :] # (n_strata_current,) - # Create h_all of length max_nb_str, fill with NaN + raw_h = h_mat[ig, :] h_all = np.full((max_nb_str,), np.nan, dtype=np.float32) n_strata_current = raw_h.shape[0] @@ -359,20 +354,15 @@ h_valid = h_all[valid_mask] top_h = np.max(h_valid) - - # Find i_zmax = number of ref_grid levels z <= top_h i_zmax = np.searchsorted(ref_grid, top_h, side='right') top_index[ig, t_idx, isl] = i_zmax - if i_zmax == 0: - # top_h < ref_grid[0], skip interpolation continue for prop, arr in prop_arrays.items(): - prop_profile_all = arr[ig, t_idx, isl, :] # (max_nb_str,) - prop_profile = prop_profile_all[valid_mask] # (n_valid_strata,) + prop_profile_all = arr[ig, t_idx, isl, :] + prop_profile = prop_profile_all[valid_mask] if prop_profile.size == 0: continue - # Step‐wise interpolation (kind='next') f_interp = interp1d( h_valid, @@ -382,6 +372,4 @@ fill_value=-1.0 ) - - # Evaluate for ref_grid[0:i_zmax] gridded_data[prop][ig, t_idx, isl, :i_zmax] = f_interp(ref_grid[:i_zmax]) @@ -399,20 +387,10 @@ ): """ - For each grid point (ig) and slope (isl), generate a 2×2 figure: + For each grid point and slope, generate a 2×2 figure of: - CO2 ice fraction - H2O ice fraction - Dust fraction - Pore fraction - - Fractions are in [0..1]. Values < 0 (fill) are masked. - Using top_index, any elevation above the last stratum is forced to NaN (white). - - Additionally, draw horizontal violet line segments at each stratum top elevation - only over the interval [date_time[t_idx], date_time[t_idx+1]] where that stratum - exists at time t_idx. This way, boundaries appear only where the strata exist. - """ - import numpy as np - import matplotlib.pyplot as plt - + """ prop_names = ['co2_ice', 'h2o_ice', 'dust', 'pore'] titles = ["CO2 ice", "H2O ice", "Dust", "Pore"] @@ -426,7 +404,4 @@ if exclude_sub: positive_indices = np.where(ref_grid >= 0.0)[0] - if positive_indices.size == 0: - print("Warning: no positive elevations in ref_grid → nothing to display.") - return sub_ref_grid = ref_grid[positive_indices] else: @@ -438,16 +413,12 @@ fig, axes = plt.subplots(2, 2, figsize=(10, 8)) fig.suptitle( - f"Content variation over time for (Grid Point {ig+1}, Slope {isl+1})", + f"Content variation over time for (Grid point {ig+1}, Slope {isl+1})", fontsize=14 ) - # For each time step t_idx, gather this stratum's valid tops - # and draw a line segment from date_time[t_idx] to date_time[t_idx+1]. - # We'll skip t_idx = ntime - 1 since no next point. - # Precompute, for each t_idx, the array of valid top elevations: + # Precompute valid stratum tops per time valid_tops_per_time = [] for t_idx in range(ntime): - raw_h = heights_data[t_idx][isl][ig, :] # (n_strata_current,) - # Exclude NaNs or zeros + raw_h = heights_data[t_idx][isl][ig, :] h_all = raw_h[~np.isnan(raw_h)] if exclude_sub: @@ -457,12 +428,10 @@ for idx, prop in enumerate(prop_names): ax = axes.flat[idx] - data_3d = gridded_data[prop][ig, :, isl, :] # shape (ntime, nz) - mat_full = data_3d.T # shape (nz, ntime) - mat = mat_full[positive_indices, :].copy() # (nz_pos, ntime) - - # Mask fill values (< 0) as NaN + data_3d = gridded_data[prop][ig, :, isl, :] + mat_full = data_3d.T + mat = mat_full[positive_indices, :].copy() mat[mat < 0.0] = np.nan - # Mask everything above the top stratum using top_index + # Mask above top stratum for t_idx in range(ntime): i_zmax = top_index[ig, t_idx, isl] @@ -473,5 +442,4 @@ mat[count_z:, t_idx] = np.nan - # Draw pcolormesh im = ax.pcolormesh( date_time, @@ -484,40 +452,10 @@ ) ax.set_title(titles[idx], fontsize=12) - ax.set_xlabel("Time (y)") + ax.set_xlabel("Time (Mars years)") ax.set_ylabel("Elevation (m)") - # Draw horizontal violet segments only where strata exist - for t_idx in range(ntime - 1): - h_vals = valid_tops_per_time[t_idx] - if h_vals.size == 0: - continue - t_left = date_time[t_idx] - t_right = date_time[t_idx + 1] - for h in h_vals: - # Only draw if h falls within sub_ref_grid - if h < sub_ref_grid[0] or h > sub_ref_grid[-1]: - continue - ax.hlines( - y=h, - xmin=t_left, - xmax=t_right, - color='violet', - linewidth=0.7, - linestyle='-' - ) - - # Reserve extra space on the right for the colorbar fig.subplots_adjust(right=0.88) - - # Place a single shared colorbar in its own axes cbar_ax = fig.add_axes([0.90, 0.15, 0.02, 0.7]) - fig.colorbar( - im, - cax=cbar_ax, - orientation='vertical', - label="Content" - ) - - # Tight layout excluding the region we reserved (0.88) + fig.colorbar(im, cax=cbar_ax, orientation='vertical', label="Content") fig.tight_layout(rect=[0, 0, 0.88, 1.0]) @@ -526,91 +464,180 @@ ) fig.savefig(fname, dpi=150) - plt.show() - plt.close(fig) + + +def plot_strata_count_and_total_height(heights_data, date_time, output_folder="."): + """ + For each grid point and slope, plot: + - Number of strata vs time + - Total deposit height vs time + """ + ntime = len(heights_data) + nslope = len(heights_data[0]) + ngrid = heights_data[0][0].shape[0] + + for ig in range(ngrid): + for isl in range(nslope): + n_strata_t = np.zeros(ntime, dtype=int) + total_height_t = np.zeros(ntime, dtype=float) + + for t_idx in range(ntime): + h_mat = heights_data[t_idx][isl] + raw_h = h_mat[ig, :] + valid_mask = (~np.isnan(raw_h)) & (raw_h != 0.0) + if np.any(valid_mask): + h_valid = raw_h[valid_mask] + n_strata_t[t_idx] = h_valid.size + total_height_t[t_idx] = np.max(h_valid) + + fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True) + fig.suptitle( + f"Strata count & total height over time for (Grid point {ig+1}, Slope {isl+1})", + fontsize=14 + ) + + axes[0].plot(date_time, n_strata_t, marker='+', linestyle='-') + axes[0].set_ylabel("Number of strata") + axes[0].grid(True) + + axes[1].plot(date_time, total_height_t, marker='+', linestyle='-') + axes[1].set_xlabel("Time (Mars years)") + axes[1].set_ylabel("Total height (m)") + axes[1].grid(True) + + fig.tight_layout(rect=[0, 0, 1, 0.95]) + fname = os.path.join( + output_folder, f"strata_count_height_ig{ig+1}_is{isl+1}.png" + ) + fig.savefig(fname, dpi=150) + + +def read_orbital_data(orb_file, martian_to_earth): + """ + Read the .asc file containing obliquity, eccentricity and Ls p. + Columns: + 0 = time in thousand Martian years + 1 = obliquity (deg) + 2 = eccentricity + 3 = Ls p (deg) + Converts times to Earth years. + """ + data = np.loadtxt(orb_file) + dates_mka = data[:, 0] + dates_yr = dates_mka * 1e3 / martian_to_earth + obliquity = data[:, 1] + eccentricity = data[:, 2] + lsp = data[:, 3] + return dates_yr, obliquity, eccentricity, lsp + + +def plot_orbital_parameters(infofile, orb_file, date_time, output_folder="."): + """ + Plot the evolution of obliquity, eccentricity and Ls p + versus simulated time. + """ + # Read conversion factor from infofile + _, martian_to_earth = read_infofile(infofile) + + # Read orbital data + dates_yr, obl, ecc, lsp = read_orbital_data(orb_file, martian_to_earth) + + # Interpolate orbital parameters at simulation dates (date_time) + obl_interp = interp1d(dates_yr, obl, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time) + ecc_interp = interp1d(dates_yr, ecc, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time) + lsp_interp = interp1d(dates_yr, lsp, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time) + + # Plot + fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True) + fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14) + + axes[0].plot(date_time, obl_interp, 'r+', linestyle='-') + axes[0].set_ylabel("Obliquity (°)") + axes[0].grid(True) + + axes[1].plot(date_time, ecc_interp, 'b+', linestyle='-') + axes[1].set_ylabel("Eccentricity") + axes[1].grid(True) + + axes[2].plot(date_time, lsp_interp, 'g+', linestyle='-') + axes[2].set_ylabel("Ls p (°)") + axes[2].set_xlabel("Time (Mars years)") + axes[2].grid(True) + + plt.tight_layout(rect=[0, 0, 1, 0.96]) + fname = os.path.join(output_folder, "orbital_parameters.png") + fig.savefig(fname, dpi=150) def main(): # 1) Get user inputs - folder_path, base_name, infofile = get_user_inputs() + folder_path, base_name, infofile, orbfile = get_user_inputs() # 2) List and verify NetCDF files files = list_netcdf_files(folder_path, base_name) if not files: - print(f"No NetCDF files named \"{base_name}#.nc\" found in \"{folder_path}\". Exiting.") + print(f"No NetCDF files named \"{base_name}#.nc\" found in \"{folder_path}\".") sys.exit(1) - nfile = len(files) - print(f"> Found {nfile} NetCDF file(s) matching \"{base_name}#.nc\".") - - # 3) Open one sample to get ngrid, nslope, lon/lat + print(f"> Found {len(files)} NetCDF file(s).") + + # 3) Open one sample to get grid dimensions & coordinates sample_file = files[0] ngrid, nslope, longitude, latitude = open_sample_dataset(sample_file) - print(f"> ngrid = {ngrid}") - print(f"> nslope = {nslope}") - - # 4) Scan all files to collect variable info + global min/max elevations - var_info, max_nb_str, min_base_elev, max_top_elev = collect_stratification_variables( - files, base_name - ) - print(f"> max(nb_str_max) = {max_nb_str}") - print(f"> min(base_elevation) = {min_base_elev:.3f}") - print(f"> max(top_elevation) = {max_top_elev:.3f}") - - # 5) Open all datasets for extraction + print(f"> ngrid = {ngrid}, nslope = {nslope}") + + # 4) Collect variable info + global min/max elevations + var_info, max_nb_str, min_base_elev, max_top_elev = collect_stratification_variables(files, base_name) + print(f"> max strata per slope = {max_nb_str}") + print(f"> min base elev = {min_base_elev:.3f} m, max top elev = {max_top_elev:.3f} m") + + # 5) Load full datasets datasets = load_full_datasets(files) - # 6) Extract raw stratification data - heights_data, raw_prop_arrays, ntime = extract_stratification_data( - datasets, var_info, ngrid, nslope, max_nb_str - ) - - # 7) Close all datasets + # 6) Extract stratification data + heights_data, raw_prop_arrays, ntime = extract_stratification_data(datasets, var_info, ngrid, nslope, max_nb_str) + + # 7) Close datasets for ds in datasets: ds.close() - # 8) Normalize raw prop arrays to volume fractions + # 8) Normalize to fractions frac_arrays = normalize_to_fractions(raw_prop_arrays) - # 9) Ask whether to show subsurface + # 9) Ask whether to include subsurface show_subsurface = get_yes_no_input("Show subsurface layers?") exclude_sub = not show_subsurface - if exclude_sub: min_base_for_interp = 0.0 - print("> Will interpolate only elevations >= 0 m (surface strata).") + print("> Interpolating only elevations >= 0 m (surface strata).") else: min_base_for_interp = min_base_elev - print(f"> Will interpolate full depth (min base = {min_base_elev:.3f} m).") - - # 10) Prompt for discretization step + print(f"> Interpolating full depth down to {min_base_elev:.3f} m.") + + # 10) Prompt discretization step dz = prompt_discretization_step(max_top_elev) - # 11) Build reference grid and interpolate (returns top_index as well) + # 11) Build reference grid and interpolate ref_grid, gridded_data, top_index = interpolate_data_on_refgrid( - heights_data, - frac_arrays, - min_base_for_interp, - max_top_elev, - dz, - exclude_sub=exclude_sub + heights_data, frac_arrays, min_base_for_interp, max_top_elev, dz, exclude_sub=exclude_sub ) - # 12) Read time stamps from "info_PEM.txt" - date_time = read_infofile(infofile) + # 12) Read timestamps and conversion factor from infofile + date_time, martian_to_earth = read_infofile(infofile) if date_time.size != ntime: - print( - "Warning: number of timestamps does not match number of NetCDF files " - f"({date_time.size} vs {ntime})." - ) - - # 13) Plot and save figures (passing top_index and heights_data) + print(f"Warning: {date_time.size} timestamps vs {ntime} NetCDF files.") + + # 13) Plot stratification over time plot_stratification_over_time( - gridded_data, - ref_grid, - top_index, - heights_data, - date_time, - exclude_sub=exclude_sub, - output_folder="." + gridded_data, ref_grid, top_index, heights_data, date_time, + exclude_sub=exclude_sub, output_folder="." ) + + # 14) Plot strata count & total height vs time + plot_strata_count_and_total_height(heights_data, date_time, output_folder=".") + + # 15) Plot orbital parameters + plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".") + + # 16) Show all figures + plt.show()