Mercurial > public > eldamwl_plotting / file comparison
comparison: eldamwl_plot.py
eldamwl_plot.py
- changeset 0
- fce4cae19357
- child 2
- 763a7ac22031
equal
deleted
inserted
replaced
| -1:000000000000 | 0:fce4cae19357 |
|---|---|
| 1 from typing import List, Any | |
| 2 | |
| 3 import xarray as xr | |
| 4 from config import * | |
| 5 import plotly.graph_objects as go | |
| 6 from plotly.subplots import make_subplots | |
| 7 import netCDF4 as nc | |
| 8 from datetime import datetime | |
| 9 import numpy as np | |
| 10 import sys | |
| 11 from os.path import exists | |
| 12 | |
| 13 | |
| 14 def read_variable_in_group(nc_file_path, group_name, variable_name): | |
| 15 try: | |
| 16 with nc.Dataset(nc_file_path, 'r') as dataset: | |
| 17 # Find the group | |
| 18 group = dataset.groups[group_name] | |
| 19 | |
| 20 # Access the variable | |
| 21 try: | |
| 22 variable = group.variables[variable_name] | |
| 23 dimensions = variable.dimensions | |
| 24 # print("Dimensions:", dimensions) | |
| 25 # variable_metadata = { variable.name, variable.long_name, variable.units, variable.ndim} | |
| 26 | |
| 27 # variable_metadata.add = variable.long_name | |
| 28 variable_metadata = \ | |
| 29 {'name': variable.name.capitalize(), | |
| 30 'long_name': variable.long_name.capitalize(), | |
| 31 'units': variable.units, | |
| 32 'ndim': variable.ndim} | |
| 33 | |
| 34 # Create an index array for slicing the variable | |
| 35 index_array = [slice(None)] * variable.ndim | |
| 36 # Access the variable data using the index array | |
| 37 variable_data = variable[index_array] | |
| 38 | |
| 39 # Fetch the dimension values | |
| 40 dimension_values = {} | |
| 41 for dim_name in dimensions: | |
| 42 dim_values = group.variables[dim_name][:] | |
| 43 dimension_values[dim_name] = dim_values | |
| 44 # Adjust the index array for slicing | |
| 45 index_array[variable.dimensions.index(dim_name)] = 0 | |
| 46 | |
| 47 # Now you have the variable data and its associated dimensions in a dictionary | |
| 48 result = {'data': variable_data, 'dimensions': dimension_values, 'metadata': variable_metadata} | |
| 49 | |
| 50 # Print or use the result as needed | |
| 51 # print(result) | |
| 52 return result | |
| 53 | |
| 54 except KeyError: | |
| 55 print("The specified variable does not exist in the group.") | |
| 56 | |
| 57 except FileNotFoundError: | |
| 58 print(f"The file '{nc_file_path}' does not exist.") | |
| 59 except Exception as e: | |
| 60 print(f"The group {e} does not exist.") | |
| 61 | |
| 62 | |
| 63 def read_variable_in_group2(nc_file_path, group_name, variable_name): | |
| 64 try: | |
| 65 with nc.Dataset(nc_file_path, 'r') as dataset: | |
| 66 # Find the group | |
| 67 group = dataset.groups[group_name] | |
| 68 | |
| 69 # Access the variable | |
| 70 try: | |
| 71 variable = group.variables[variable_name] | |
| 72 dimensions = variable.dimensions | |
| 73 # print("Dimensions:", dimensions) | |
| 74 # variable_metadata = { variable.name, variable.long_name, variable.units, variable.ndim} | |
| 75 | |
| 76 # variable_metadata.add = variable.long_name | |
| 77 variable_metadata = \ | |
| 78 {'name': variable.name.capitalize(), | |
| 79 'long_name': variable.long_name.capitalize(), | |
| 80 'units': variable.units, | |
| 81 'ndim': variable.ndim} | |
| 82 | |
| 83 # Create an index array for slicing the variable | |
| 84 index_array = [slice(None)] * variable.ndim | |
| 85 # Access the variable data using the index array | |
| 86 variable_data = variable[index_array] | |
| 87 | |
| 88 # Fetch the dimension values | |
| 89 dimension_values = {} | |
| 90 for dim_name in dimensions: | |
| 91 dim_values = group.variables[dim_name][:] | |
| 92 dimension_values[dim_name] = dim_values | |
| 93 # Adjust the index array for slicing | |
| 94 index_array[variable.dimensions.index(dim_name)] = 0 | |
| 95 | |
| 96 # Now you have the variable data and its associated dimensions in a dictionary | |
| 97 # result = {'data': variable_data, 'dimensions': dimension_values, 'metadata': variable_metadata} | |
| 98 | |
| 99 # Print or use the result as needed | |
| 100 # print(result) | |
| 101 return variable_data, dimension_values, variable_metadata | |
| 102 | |
| 103 except KeyError: | |
| 104 print(f"The variable {variable_name} does not exist in the {group_name} group.") | |
| 105 except FileNotFoundError: | |
| 106 print(f"The file {nc_file_path} does not exist.") | |
| 107 except Exception as e: | |
| 108 print(f"The group {e} does not exist.") | |
| 109 | |
| 110 | |
| 111 def plot_vertical_profiles_plotly(variables, data, dimensions, y_array, positive_error_array, negative_error_array, | |
| 112 xaxis_label, yaxis_label, title, timeinrows): | |
| 113 """ | |
| 114 ToDo Complete | |
| 115 Plot vertical profiles using plotly | |
| 116 | |
| 117 Parameters: | |
| 118 - | |
| 119 | |
| 120 Returns: | |
| 121 - fig: html | |
| 122 """ | |
| 123 | |
| 124 # Set number of columns, number of rows, and subtitles if needed | |
| 125 ncols = len(variables) | |
| 126 subtitles = None | |
| 127 if timeinrows is True: | |
| 128 nrows = len(dimensions['time']) | |
| 129 if nrows > 1: | |
| 130 for t in range(nrows): | |
| 131 for v in range(ncols): | |
| 132 if subtitles is None: | |
| 133 subtitles = [str(datetime.fromtimestamp(dimensions['time'][t]))] | |
| 134 else: | |
| 135 subtitles += [str(datetime.fromtimestamp(dimensions['time'][t]))] | |
| 136 else: | |
| 137 nrows = 1 | |
| 138 subtitles = '' | |
| 139 | |
| 140 # Create figure | |
| 141 fig = make_subplots(rows=nrows, cols=ncols, | |
| 142 subplot_titles=subtitles, | |
| 143 shared_yaxes=True) | |
| 144 | |
| 145 # Define colours | |
| 146 clrs = [['rgb(0,35,255)', 'rgb(0,20,150)', 'rgb(0,170,250)'], # UV | |
| 147 ['rgb(90,210,50)', 'rgb(40,90,25)', 'rgb(120,250,80)'], # VIS | |
| 148 ['rgb(250,30,30)', 'rgb(200,25,30)', 'rgb(250,125,0)']] # IR | |
| 149 clrs = np.array(clrs) | |
| 150 | |
| 151 # Plot | |
| 152 for v in range(len(variables)): | |
| 153 for t in range(len(dimensions['time'])): | |
| 154 for w in range(len(dimensions['wavelength'])): | |
| 155 tn = str(datetime.fromtimestamp(dimensions['time'][t])) | |
| 156 wn = str(dimensions['wavelength'][w]) | |
| 157 if timeinrows == True: | |
| 158 lgnd = wn + ' nm' | |
| 159 rown = t + 1 | |
| 160 colorint = 1 | |
| 161 colorind = 0 | |
| 162 subtitle = str(datetime.fromtimestamp(dimensions['time'][t])) | |
| 163 else: | |
| 164 lgnd = wn + ' nm - ' + tn | |
| 165 rown = 1 | |
| 166 colorind = t | |
| 167 colorint = 1 - (0.4 * t) | |
| 168 | |
| 169 # colour definition | |
| 170 if abs(dimensions['wavelength'][w] - 1064) < 1: | |
| 171 clr = clrs[2, colorind] | |
| 172 # clr = 'rgb(235,70,20)' # 'rgba(215,27,96,0.4)' | |
| 173 elif abs(dimensions['wavelength'][w] - 532) < 1: | |
| 174 clr = clrs[1, colorind] | |
| 175 elif abs(dimensions['wavelength'][w] - 355) < 1: | |
| 176 clr = clrs[0, colorind] | |
| 177 else: | |
| 178 clr = 'rgb(0,0,0)' | |
| 179 | |
| 180 data_values = data[variables[v]][w, t, :] | |
| 181 error_values = positive_error_array[variables[v]][w, t, :] | |
| 182 vertical_levels = y_array["data"][t, :] | |
| 183 | |
| 184 # Add trace | |
| 185 fig.add_trace(go.Scatter( | |
| 186 x=data_values, | |
| 187 y=vertical_levels, | |
| 188 error_x=dict(array=error_values, width=0, thickness=0.1), | |
| 189 # ToDo Include distinction between positive and negative errors | |
| 190 mode='lines', | |
| 191 line=dict(width=2, color=clr), | |
| 192 name=lgnd), | |
| 193 row=rown, col=v + 1, | |
| 194 ) | |
| 195 | |
| 196 # ToDo Check formatting options (showlegend, mapbox) | |
| 197 fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1) | |
| 198 fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown, | |
| 199 col=v + 1) | |
| 200 | |
| 201 fig.update_xaxes(ticks="outside") | |
| 202 fig.update_yaxes(ticks="outside", row=rown, col=v + 1) # ToDo check tickformat (i.e. tickformat=":.2e") | |
| 203 if v == 0: | |
| 204 fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1) | |
| 205 | |
| 206 # Set axis labels and title | |
| 207 fig.update_layout( | |
| 208 yaxis_title=yaxis_label, | |
| 209 title=title, | |
| 210 template="seaborn", # You can change the template as per your preference | |
| 211 # ['ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation', | |
| 212 # 'xgridoff', 'ygridoff', 'gridon', 'none'] | |
| 213 height=max(600, 450 * rown), width=min(max(400 * len(variables), 750), 1400), | |
| 214 ) | |
| 215 | |
| 216 return fig | |
| 217 | |
| 218 | |
| 219 # Files for testing | |
| 220 # netcdf_file = 'hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc' # Backscatter and extinction, no error | |
| 221 # netcdf_file = 'hpb_012_0000627_202308240200_202308240400_20230824hpb0200_ELDAmwl_v0.0.1.nc' # Extinction with error | |
| 222 | |
| 223 | |
| 224 try: | |
| 225 netcdf_file = sys.argv[1] | |
| 226 except: | |
| 227 print('Syntax is incorrect, you need to pass the filename of the file to be plotted.\n' | |
| 228 'Example: ./eldamwl_plot.py hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc') | |
| 229 exit() | |
| 230 | |
| 231 netcdf_file_path = netcdf_path + netcdf_file | |
| 232 | |
| 233 file_exists = exists(netcdf_file_path) | |
| 234 if not file_exists: | |
| 235 print('The file that you provided does not exist.') | |
| 236 exit() | |
| 237 | |
| 238 # Groups of data in the NetCDF file (low and high resolution) | |
| 239 groups = ['lowres_products', 'highres_products'] | |
| 240 groups_names = ['Low resolution', 'High resolution'] | |
| 241 | |
| 242 # Variables that we want to plot | |
| 243 variable = ['backscatter', 'extinction', 'lidarratio', 'volumedepolarization', 'particledepolarization'] # ToDo Test particledepolarization when available | |
| 244 # Errors associated to the variables | |
| 245 variable_error = ['error_backscatter', 'error_extinction', 'error_lidarratio', | |
| 246 'positive_systematic_error_volumedepolarization', 'error_particledepolarization'] | |
| 247 variable_negative_error_volumedepol = 'negative_systematic_error_volumedepolarization' | |
| 248 # Other variables: | |
| 249 # cloud_mask (time, level) | |
| 250 | |
| 251 | |
| 252 for g in range(len(groups)): | |
| 253 # Read the altitude variable | |
| 254 altitude = read_variable_in_group(netcdf_file_path, groups[g], 'altitude') | |
| 255 | |
| 256 if altitude['metadata']['units'] == 'm': | |
| 257 altitude['data'] = altitude['data'] / 1000.0 | |
| 258 altitude['metadata']['units'] = 'km' | |
| 259 | |
| 260 # Initialize variables | |
| 261 data = None | |
| 262 variables_with_data = None | |
| 263 variables_axis = None | |
| 264 | |
| 265 # Read the data | |
| 266 for v in range(len(variable)): | |
| 267 # Initialize var | |
| 268 var = None | |
| 269 | |
| 270 try: | |
| 271 var, dim, met = read_variable_in_group2(netcdf_file_path, groups[g], variable[v]) | |
| 272 except TypeError: | |
| 273 print(f"The variable {variable[v]} is not present. Cannot unpack non-iterable NoneType object") | |
| 274 | |
| 275 if var is not None: # If the variable that we're trying to read exists | |
| 276 if var.max() != var.min(): # If there is data in the array | |
| 277 if variables_with_data is None: | |
| 278 variables_with_data = [variable[v]] | |
| 279 variables_axis = [met['long_name'] + ' [ ' + met['units'] + ' ]'] | |
| 280 else: | |
| 281 variables_with_data += [variable[v]] | |
| 282 variables_axis += [met['long_name'] + ' [ ' + met['units'] + ' ]'] | |
| 283 | |
| 284 var_error, dim_error, met_error = read_variable_in_group2(netcdf_file_path, groups[g], | |
| 285 variable_error[v]) | |
| 286 if variable[v] == 'volumedepolarization': | |
| 287 var_error_negative, dim, met = read_variable_in_group2(netcdf_file_path, groups[g], | |
| 288 variable_negative_error_volumedepol) | |
| 289 else: | |
| 290 var_error_negative = var_error | |
| 291 | |
| 292 # Add read data to the data variable to have it all grouped | |
| 293 if data is None: | |
| 294 # Start populating data and data_error | |
| 295 data = {variable[v]: var} | |
| 296 data_error_positive = {variable[v]: var_error} | |
| 297 data_error_negative = {variable[v]: var_error_negative} | |
| 298 else: | |
| 299 # Add more data to data and data_error | |
| 300 data[variable[v]] = var | |
| 301 data_error_positive[variable[v]] = var_error | |
| 302 data_error_negative[variable[v]] = var_error_negative | |
| 303 else: | |
| 304 print(f'There is no data for {variable[v]} in the {groups_names[g]} group.') | |
| 305 | |
| 306 dimensions = dim | |
| 307 | |
| 308 # Ask user if he/she wants to plot different temporal profiles in different rows | |
| 309 if g == 0: | |
| 310 # Initialize | |
| 311 timeinplots = False | |
| 312 timeinrows = False | |
| 313 # Ask user how he/she wants the plot | |
| 314 if len(dimensions['time']) > 1: | |
| 315 answer = input( | |
| 316 f"There are {len(dimensions['time'])} temporal profiles, do you want to plot them on different plots [y/n]? ") | |
| 317 if answer[0:1].lower() == 'y': | |
| 318 timeinplots = True | |
| 319 if timeinplots == False: | |
| 320 answer = input( | |
| 321 f"There are {len(dimensions['time'])} temporal profiles, do you want to plot them on different rows [y/n]? ") | |
| 322 if answer[0:1].lower() == 'y': | |
| 323 timeinrows = True | |
| 324 | |
| 325 if variables_with_data is not None: | |
| 326 timeiter = dimensions['time'] | |
| 327 if timeinplots == True: | |
| 328 dimensions_t = dimensions | |
| 329 for t in range(len(dimensions['time'])): | |
| 330 dimensions_t['time'] = [timeiter[t]] | |
| 331 for vr in range(len(variables_with_data)): | |
| 332 if vr == 0: | |
| 333 data_t = {variables_with_data[vr]: data[variables_with_data[vr]][:, t:t + 1, :]} | |
| 334 else: | |
| 335 data_t[variables_with_data[vr]] = data[variables_with_data[vr]][:, t:t + 1, :] | |
| 336 | |
| 337 fig = plot_vertical_profiles_plotly( | |
| 338 variables=variables_with_data, | |
| 339 data=data_t, | |
| 340 dimensions=dimensions_t, | |
| 341 y_array=altitude, | |
| 342 positive_error_array=data_error_positive, | |
| 343 negative_error_array=data_error_negative, | |
| 344 xaxis_label=variables_axis, | |
| 345 yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']', | |
| 346 title=netcdf_file[0:3].upper() + " - " + str( | |
| 347 datetime.fromtimestamp(dimensions_t['time'][0]).strftime( | |
| 348 '%d/%m/%Y %H:%M:%S')) + " - ELDA MWL " + | |
| 349 groups_names[g], | |
| 350 timeinrows=True | |
| 351 ) | |
| 352 | |
| 353 # Show the plot (in Jupyter Notebook or Jupyter Lab) | |
| 354 fig.show() | |
| 355 | |
| 356 # Save the plot | |
| 357 filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( | |
| 358 dimensions['time'][0]).strftime('%Y%m%d-%H%M%S') + "_" + groups[g] + ".html" | |
| 359 fig.write_html(html_path + filename) | |
| 360 | |
| 361 else: | |
| 362 # Plotting in one plot (even if there are different time profiles) | |
| 363 fig = plot_vertical_profiles_plotly( | |
| 364 variables_with_data, | |
| 365 data, | |
| 366 dimensions, | |
| 367 altitude, | |
| 368 data_error_positive, | |
| 369 data_error_negative, | |
| 370 xaxis_label=variables_axis, | |
| 371 yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']', | |
| 372 title=netcdf_file[0:3].upper() + " - " + str( | |
| 373 datetime.fromtimestamp(dimensions['time'][0]).strftime('%d/%m/%Y')) + " - ELDA MWL " + groups_names[ | |
| 374 g], | |
| 375 timeinrows=timeinrows | |
| 376 ) | |
| 377 | |
| 378 # Show the plot (in Jupyter Notebook or Jupyter Lab) | |
| 379 fig.show() | |
| 380 | |
| 381 # Save the plot | |
| 382 if timeinrows == True: | |
| 383 filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( | |
| 384 dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp( | |
| 385 dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + "_timeinrows.html" | |
| 386 else: | |
| 387 filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( | |
| 388 dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp( | |
| 389 dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + ".html" | |
| 390 | |
| 391 fig.write_html(html_path + filename) | |
| 392 | |
| 393 else: | |
| 394 print(f'There is no data to be plotted in the {groups_names[g]} group.') |
