- Timestamp:
- May 21, 2025, 4:01:25 PM (4 weeks ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py
r3673 r3771 4 4 5 5 import os 6 import sys 6 7 import numpy as np 7 8 from netCDF4 import Dataset … … 9 10 from scipy import interpolate 10 11 11 ######################### 12 ### Parameters to fill in 13 folder_path = "starts" 14 base_name = "restartpem" 15 dz = 0.1 # Discrization step of the reference grid for the elevation [m] 16 infofile = 'info_PEM.txt' 17 ######################### 18 19 20 ############################################################################################## 12 ### Function to get inputs 13 def get_user_inputs(): 14 folder_path = input("Enter the folder path containing the NetCDF files (press the Enter key for default [starts]): ").strip() 15 if not folder_path: 16 folder_path = "starts" 17 while not os.path.isdir(folder_path): 18 print("Invalid folder path. Please try again.") 19 folder_path = input("Enter the folder path containing the NetCDF files (press the Enter key for default [starts]): ").strip() 20 if not folder_path: 21 folder_path = "starts" 22 23 base_name = input("Enter the base name of the NetCDF files (press the Enter key for default [restartpem]): ").strip() 24 if not base_name: 25 base_name = "restartpem" 26 27 infofile = input("Enter the name of the PEM info file (press the Enter key for default [info_PEM.txt]): ").strip() 28 if not infofile: 29 infofile = "info_PEM.txt" 30 while not os.path.isfile(infofile): 31 print("Invalid file path. Please try again.") 32 infofile = input("Enter the name of the PEM info file (press the Enter key for default [info_PEM.txt]): ").strip() 33 if not infofile: 34 infofile = "info_PEM.txt" 35 36 return folder_path, base_name, infofile 37 21 38 ### Function to read the "startpem.nc" files and process their stratification data 22 39 def process_files(folder_path,base_name): … … 35 52 # Process each file and collect data 36 53 datasets = [] 37 max_top_elevation = 0 54 min_base_elevation = -56943.759374550937 # Base elevation of the deepest subsurface layer 55 max_top_elevation = 0. 38 56 max_nb_str = 0 39 ngrid = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['physical_points'].size # ngrid is the same for all files 40 nslope = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['nslope'].size # nslope is the same for all files 41 longitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['longitude'][:] 42 latitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['latitude'][:] 57 with Dataset(os.path.join(folder_path, base_name + "1.nc"), 'r') as ds: 58 ngrid = ds.dimensions['physical_points'].size # ngrid is the same for all files 59 nslope = ds.dimensions['nslope'].size # nslope is the same for all files 60 longitude = ds.variables['longitude'][:] 61 latitude = ds.variables['latitude'][:] 62 43 63 for i in range(1,nfile + 1): 44 64 file_path = os.path.join(folder_path,base_name + str(i) + ".nc") … … 53 73 for k in range(1,nslope + 1): 54 74 slope_var_name = f"stratif_slope{k:02d}_top_elevation" 75 min_base_elevation = min(min_base_elevation,np.min(ds.variables[slope_var_name][:])) 55 76 max_top_elevation = max(max_top_elevation,np.max(ds.variables[slope_var_name][:])) 56 77 57 print(f"> number of files = {nfile}") 58 print(f"> ngrid = {ngrid}") 59 print(f"> nslope = {nslope}") 60 print(f"> max(nb_str_max) = {max_nb_str}") 61 print(f"> max(top_elevation) = {max_top_elevation}") 78 print(f"> number of files = {nfile}") 79 print(f"> ngrid = {ngrid}") 80 print(f"> nslope = {nslope}") 81 print(f"> max(nb_str_max) = {max_nb_str}") 82 print(f"> min(base_elevation) = {min_base_elevation}") 83 print(f"> max(top_elevation) = {max_top_elevation}") 62 84 63 85 # Concatenate stratif variables with dimension 'nb_str_max' along the "Time" dimension … … 68 90 stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str)) 69 91 stratif_pore = np.zeros((ngrid,nfile,nslope,max_nb_str)) 92 stratif_poreice = np.zeros((ngrid,nfile,nslope,max_nb_str)) 70 93 for var_name in datasets[0].variables: 71 94 if 'top_elevation' in var_name: … … 76 99 stratif_heights[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 77 100 print(f"Processed variable: {var_name}") 78 elif ' co2ice_volfrac' in var_name:101 elif 'h_co2ice' in var_name: 79 102 for i in range(0,ngrid): 80 103 for j in range(0,nfile): … … 83 106 stratif_co2ice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 84 107 print(f"Processed variable: {var_name}") 85 elif 'h 2oice_volfrac' in var_name:108 elif 'h_h2oice' in var_name: 86 109 for i in range(0,ngrid): 87 110 for j in range(0,nfile): … … 90 113 stratif_h2oice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 91 114 print(f"Processed variable: {var_name}") 92 elif ' dust_volfrac' in var_name:115 elif 'h_dust' in var_name: 93 116 for i in range(0,ngrid): 94 117 for j in range(0,nfile): … … 97 120 stratif_dust[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 98 121 print(f"Processed variable: {var_name}") 99 elif ' pore_volfrac' in var_name:122 elif 'h_pore' in var_name: 100 123 for i in range(0,ngrid): 101 124 for j in range(0,nfile): … … 103 126 if f'slope{k + 1:02d}' in var_name: 104 127 stratif_pore[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 128 print(f"Processed variable: {var_name}") 129 elif 'icepore_volfrac' in var_name: 130 for i in range(0,ngrid): 131 for j in range(0,nfile): 132 for k in range(0,nslope): 133 if f'slope{k + 1:02d}' in var_name: 134 stratif_poreice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i] 105 135 print(f"Processed variable: {var_name}") 106 136 … … 110 140 111 141 stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_pore] 112 return stratif_data, max_top_elevation, longitude, latitude 142 143 while True: 144 try: 145 dz = float(input("Enter the discretization step of the reference grid for the elevation [m]: ").strip()) 146 if dz <= 0: 147 print("Discretization step must be strictly positive!") 148 continue 149 if dz > max_top_elevation: 150 print("Discretization step is higher than the maximum top elevation: please provide a correct value!") 151 continue 152 break 153 except ValueError: 154 print("Invalid value.") 155 return stratif_data, min_base_elevation, max_top_elevation, longitude, latitude, dz 113 156 114 157 ### Function to interpolate the stratification data on a reference grid 115 def interpolate_data(stratif_data,m ax_top_elevation,dz):158 def interpolate_data(stratif_data,min_base_elevation,max_top_elevation,dz): 116 159 # Define the reference ref_grid 117 ref_grid = np.arange( 0,max_top_elevation,dz)160 ref_grid = np.arange(min_base_elevation,max_top_elevation,dz) 118 161 print(f"> Number of ref_grid points = {len(ref_grid)}") 119 162 … … 143 186 date_time = [] 144 187 for line in file: 145 data = list(map( int,line.split()))188 data = list(map(float,line.split())) 146 189 data_lines.append(data) 147 190 date_time.append(data[0]) … … 150 193 151 194 ### Processing 152 stratif_data, max_top_elevation, longitude, latitude = process_files(folder_path,base_name) 153 if dz > max_top_elevation: 154 print('The discretization step is higher than the maximum top elevation: please provide a correct value!') 155 exit() 156 ref_grid, gridded_stratif_data = interpolate_data(stratif_data,max_top_elevation,dz) 195 folder_path, base_name, infofile = get_user_inputs() 196 stratif_data, min_base_elevation, max_top_elevation, longitude, latitude, dz = process_files(folder_path,base_name) 197 ref_grid, gridded_stratif_data = interpolate_data(stratif_data,min_base_elevation,max_top_elevation,dz) 157 198 date_time = read_infofile(infofile) 158 199 … … 168 209 for ax in axes.flat: 169 210 time_mesh, elevation_mesh = np.meshgrid(date_time,ref_grid) 170 #im = ax.imshow(np.transpose(gridded_stratif_data[iprop][ 0,:,0,:]),aspect = 'auto',cmap = 'viridis',origin = 'lower',extent = [date_time[0],date_time[-1],ref_grid[0],ref_grid[-1]],vmin = 0,vmax = 1)211 #im = ax.imshow(np.transpose(gridded_stratif_data[iprop][igr,:,isl,:]),aspect = 'auto',cmap = cmap,origin = 'lower',extent = [date_time[0],date_time[-1],ref_grid[0],ref_grid[-1]],vmin = 0,vmax = 1) 171 212 im = ax.pcolormesh(time_mesh,elevation_mesh,np.transpose(gridded_stratif_data[iprop][igr,:,isl,:]),cmap = cmap,shading = 'auto',vmin = 0,vmax = 1) 172 213 ax.set_title(subtitle[iprop])
Note: See TracChangeset
for help on using the changeset viewer.
