diff -r f0db47eb70b1 -r 763a7ac22031 eldamwl_plot.py --- a/eldamwl_plot.py Mon Jan 29 11:16:12 2024 +0000 +++ b/eldamwl_plot.py Tue Jul 02 08:37:49 2024 +0000 @@ -1,6 +1,5 @@ -from typing import List, Any - -import xarray as xr +# import nc as nc +# import xarray as xr from config import * import plotly.graph_objects as go from plotly.subplots import make_subplots @@ -60,6 +59,19 @@ print(f"The group {e} does not exist.") +def read_global_attribute(nc_file_path, attribute_name): + # Open the NetCDF file in read mode + with nc.Dataset(nc_file_path, 'r') as dataset: + # with netCDF4.Dataset(nc_file_path, 'r') as nc: + # Check if the attribute exists + if attribute_name in dataset.ncattrs(): + # Read the attribute value + attribute_value = dataset.getncattr(attribute_name) + return attribute_value + else: + print(f"Attribute '{attribute_name}' not found in the NetCDF file.") + + def read_variable_in_group2(nc_file_path, group_name, variable_name): try: with nc.Dataset(nc_file_path, 'r') as dataset: @@ -89,7 +101,7 @@ dimension_values = {} for dim_name in dimensions: dim_values = group.variables[dim_name][:] - dimension_values[dim_name] = dim_values + dimension_values[dim_name] = dim_values.tolist() # Adjust the index array for slicing index_array[variable.dimensions.index(dim_name)] = 0 @@ -154,17 +166,15 @@ for w in range(len(dimensions['wavelength'])): tn = str(datetime.fromtimestamp(dimensions['time'][t])) wn = str(dimensions['wavelength'][w]) - if timeinrows == True: + if timeinrows: lgnd = wn + ' nm' rown = t + 1 - colorint = 1 colorind = 0 - subtitle = str(datetime.fromtimestamp(dimensions['time'][t])) + # subtitle = str(datetime.fromtimestamp(dimensions['time'][t])) else: lgnd = wn + ' nm - ' + tn rown = 1 colorind = t - colorint = 1 - (0.4 * t) # colour definition if abs(dimensions['wavelength'][w] - 1064) < 1: @@ -178,30 +188,164 @@ clr = 'rgb(0,0,0)' data_values = data[variables[v]][w, t, :] - error_values = positive_error_array[variables[v]][w, t, :] - vertical_levels = y_array["data"][t, :] + if np.average(np.ma.masked_array(data_values, np.isnan(data_values))) is not np.ma.masked: + vertical_levels = y_array["data"][t, :] + error_values = positive_error_array[variables[v]][w, t, :] + + # Add trace + fig.add_trace(go.Scatter( + x=data_values, + y=vertical_levels, + error_x=dict(array=error_values, width=0, thickness=0.1), + # ToDo Include distinction between positive and negative errors + mode='lines', + line=dict(width=2, color=clr), + name=lgnd), + row=rown, col=v + 1, + ) + + # ToDo Check formatting options (showlegend, mapbox) + fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1) + fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown, + col=v + 1) + + fig.update_xaxes(ticks="outside") + fig.update_yaxes(ticks="outside", row=rown, col=v + 1) + # ToDo check tickformat (i.e. tickformat=":.2e") + if v == 0: + fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1) + + # ToDo add boxes + a = fig.layout['xaxis']['domain'] + # print(a) + b = fig.layout['yaxis']['domain'] + # print(b) + # fig.add_hline(y=b[0], line=dict(color="yellow", width=2), row=rown, col=v+1) + # fig.add_hline(y=b[1], line=dict(color="purple", width=2), row=rown, col=v+1) + # fig.add_vline(x=a[0], line=dict(color="orange", width=2), row=rown, col=v+1) + # fig.add_vline(x=a[1], line=dict(color="pink", width=2), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=a[0], y0=b[0], x1=a[1], y1=b[1], + # line=dict(color="grey", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=1, y0=1, x1=1, y1=0, + # line=dict(color="blue", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=1, y1=0, + # line=dict(color="red", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=0, y1=1, + # line=dict(color="green", width=2, ), row=rown, col=v+1) + # fig.add_hline(y=0, line=dict(color="yellow", width=2)) + # print(fig.data[0]['x']) + + # Set axis labels and title + fig.update_layout( + yaxis_title=yaxis_label, + title=title, + template="seaborn", # You can change the template as per your preference + # ['ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation', + # 'xgridoff', 'ygridoff', 'gridon', 'none'] + height=max(600, 450 * rown), width=min(max(400 * len(variables), 750), 1400), + ) + + return fig + + +def plot_vertical_profiles_plotly_allres(variables, data, dimensions, y_array, positive_error_array, + negative_error_array, xaxis_label, yaxis_label, title): + """ + ToDo Complete + Plot vertical profiles using plotly + + Parameters: + - + + Returns: + - fig: html + """ + + # Set number of columns, number of rows, and subtitles if needed + ncols = len(variables) + nrows = 1 + subtitles = '' - # Add trace - fig.add_trace(go.Scatter( - x=data_values, - y=vertical_levels, - error_x=dict(array=error_values, width=0, thickness=0.1), - # ToDo Include distinction between positive and negative errors - mode='lines', - line=dict(width=2, color=clr), - name=lgnd), - row=rown, col=v + 1, - ) + # Create figure + fig = make_subplots(rows=nrows, cols=ncols, + subplot_titles=subtitles, + shared_yaxes=True) + + # Define colours + clrs = [['rgb(0,35,255)', 'rgb(0,20,150)'], # UV + ['rgb(90,210,50)', 'rgb(40,90,25)'], # VIS + ['rgb(250,30,30)', 'rgb(130,30,10)']] # IR + clrs = np.array(clrs) + + # Plot + for v in range(len(variables)): + for t in range(len(dimensions['time'])): + for w in range(len(dimensions['wavelength_res'])): + wl = float(dimensions['wavelength_res'][w].split(' ')[0]) + wn = str(dimensions['wavelength_res'][w]) + lgnd = wn.replace(' ', ' nm ') + rown = 1 + if wn.find('LR') > 0: + colorind = 0 + else: + colorind = 1 + + # colour definition + if abs(wl - 1064) < 1: + clr = clrs[2, colorind] + elif abs(wl - 532) < 1: + clr = clrs[1, colorind] + elif abs(wl - 355) < 1: + clr = clrs[0, colorind] + else: + clr = 'rgb(0,0,0)' - # ToDo Check formatting options (showlegend, mapbox) - fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1) - fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown, - col=v + 1) + data_values = data[variables[v]][w, t, :] + if np.average(np.ma.masked_array(data_values, np.isnan(data_values))) is not np.ma.masked: + vertical_levels = y_array["data"][t, :] + error_values = positive_error_array[variables[v]][w, t, :] + # Add trace + fig.add_trace(go.Scatter( + x=data_values, + y=vertical_levels, + error_x=dict(array=error_values, width=0, thickness=0.1), + # ToDo Include distinction between positive and negative errors + mode='lines', + line=dict(width=2, color=clr), + name=lgnd), + row=rown, col=v + 1, + ) + + # ToDo Check formatting options (showlegend, mapbox) + fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1) + fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown, + col=v + 1) - fig.update_xaxes(ticks="outside") - fig.update_yaxes(ticks="outside", row=rown, col=v + 1) # ToDo check tickformat (i.e. tickformat=":.2e") - if v == 0: - fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1) + fig.update_xaxes(ticks="outside") + fig.update_yaxes(ticks="outside", row=rown, col=v + 1) + # ToDo check tickformat (i.e. tickformat=":.2e") + if v == 0: + fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1) + + # ToDo add boxes + a = fig.layout['xaxis']['domain'] + # print(a) + b = fig.layout['yaxis']['domain'] + # print(b) + # fig.add_hline(y=b[0], line=dict(color="yellow", width=2), row=rown, col=v+1) + # fig.add_hline(y=b[1], line=dict(color="purple", width=2), row=rown, col=v+1) + # fig.add_vline(x=a[0], line=dict(color="orange", width=2), row=rown, col=v+1) + # fig.add_vline(x=a[1], line=dict(color="pink", width=2), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=a[0], y0=b[0], x1=a[1], y1=b[1], + # line=dict(color="grey", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=1, y0=1, x1=1, y1=0, + # line=dict(color="blue", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=1, y1=0, + # line=dict(color="red", width=2, ), row=rown, col=v+1) + # fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=0, y1=1, + # line=dict(color="green", width=2, ), row=rown, col=v+1) + # fig.add_hline(y=0, line=dict(color="yellow", width=2)) + # print(fig.data[0]['x']) # Set axis labels and title fig.update_layout( @@ -217,8 +361,10 @@ # Files for testing -# netcdf_file = 'hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc' # Backscatter and extinction, no error -# netcdf_file = 'hpb_012_0000627_202308240200_202308240400_20230824hpb0200_ELDAmwl_v0.0.1.nc' # Extinction with error +# Backscatter and extinction, no error: +# netcdf_file = 'hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc' +# Extinction with error: +# netcdf_file = 'hpb_012_0000627_202308240200_202308240400_20230824hpb0200_ELDAmwl_v0.0.1.nc' try: @@ -228,6 +374,11 @@ 'Example: ./eldamwl_plot.py hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc') exit() +# Print help +if netcdf_file == '-h' or netcdf_file == '--h': + print('Syntax: ./eldamwl_plot.py hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc') + exit() + netcdf_file_path = netcdf_path + netcdf_file file_exists = exists(netcdf_file_path) @@ -240,15 +391,20 @@ groups_names = ['Low resolution', 'High resolution'] # Variables that we want to plot -variable = ['backscatter', 'extinction', 'lidarratio', 'volumedepolarization', 'particledepolarization'] # ToDo Test particledepolarization when available +variable = ['backscatter', 'extinction', 'lidarratio', 'volumedepolarization', + 'particledepolarization'] # ToDo Test particledepolarization when available # Errors associated to the variables variable_error = ['error_backscatter', 'error_extinction', 'error_lidarratio', 'positive_systematic_error_volumedepolarization', 'error_particledepolarization'] variable_negative_error_volumedepol = 'negative_systematic_error_volumedepolarization' -# Other variables: +# ToDo Consider adding other variables: # cloud_mask (time, level) +# Read global attributes +station_id = read_global_attribute(netcdf_file_path, 'station_ID') +station_location = read_global_attribute(netcdf_file_path, 'location') +# Read data for g in range(len(groups)): # Read the altitude variable altitude = read_variable_in_group(netcdf_file_path, groups[g], 'altitude') @@ -305,26 +461,48 @@ dimensions = dim - # Ask user if he/she wants to plot different temporal profiles in different rows + # Save lowres and highres in variables to plot them together afterwards + if groups[g] == 'lowres_products': + altitude_lr = altitude + dimensions_lr = dim + data_lr = data + data_error_positive_lr = data_error_positive + data_error_negative_lr = data_error_negative + variables_with_data_lr = variables_with_data + variables_axis_lr = variables_axis + if groups[g] == 'highres_products': + altitude_hr = altitude + dimensions_hr = dim + data_hr = data + data_error_positive_hr = data_error_positive + data_error_negative_hr = data_error_negative + variables_with_data_hr = variables_with_data + variables_axis_hr = variables_axis + + # If interactive is set to True, ask user if he/she wants to plot different temporal profiles in different rows + # If interactive is set to False, proceed with default values if g == 0: # Initialize - timeinplots = False + timeinplots = True timeinrows = False - # Ask user how he/she wants the plot - if len(dimensions['time']) > 1: + # If interactive (config.py) = True, ask user how he/she wants the plot + if interactive and len(dimensions['time']) > 1: answer = input( - f"There are {len(dimensions['time'])} temporal profiles, do you want to plot them on different plots [y/n]? ") - if answer[0:1].lower() == 'y': - timeinplots = True - if timeinplots == False: + f"There are {len(dimensions['time'])} temporal profiles, " + f"do you want to plot them on different plots [y/n]? ") + if answer[0:1].lower() == 'n': + timeinplots = False + if not timeinplots: answer = input( - f"There are {len(dimensions['time'])} temporal profiles, do you want to plot them on different rows [y/n]? ") + f"There are {len(dimensions['time'])} temporal profiles, " + f"do you want to plot them on different rows [y/n]? ") if answer[0:1].lower() == 'y': timeinrows = True if variables_with_data is not None: timeiter = dimensions['time'] - if timeinplots == True: + if timeinplots: + # Default plotting dimensions_t = dimensions for t in range(len(dimensions['time'])): dimensions_t['time'] = [timeiter[t]] @@ -343,15 +521,14 @@ negative_error_array=data_error_negative, xaxis_label=variables_axis, yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']', - title=netcdf_file[0:3].upper() + " - " + str( - datetime.fromtimestamp(dimensions_t['time'][0]).strftime( - '%d/%m/%Y %H:%M:%S')) + " - ELDA MWL " + - groups_names[g], + title=station_location + " (" + station_id.upper() + ") - " + str( + datetime.fromtimestamp(dimensions_t['time'][0]).strftime('%d/%m/%Y %H:%M:%S')) + + " - ELDA MWL " + groups_names[g], timeinrows=True ) # Show the plot (in Jupyter Notebook or Jupyter Lab) - fig.show() + # fig.show() # Save the plot filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( @@ -369,22 +546,22 @@ data_error_negative, xaxis_label=variables_axis, yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']', - title=netcdf_file[0:3].upper() + " - " + str( + title=station_location + " (" + station_id.upper() + ") - " + str( datetime.fromtimestamp(dimensions['time'][0]).strftime('%d/%m/%Y')) + " - ELDA MWL " + groups_names[ g], timeinrows=timeinrows ) # Show the plot (in Jupyter Notebook or Jupyter Lab) - fig.show() + # fig.show() # Save the plot - if timeinrows == True: - filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( + if timeinrows: + filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp( dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp( dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + "_timeinrows.html" else: - filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( + filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp( dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp( dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + ".html" @@ -392,3 +569,158 @@ else: print(f'There is no data to be plotted in the {groups_names[g]} group.') + +# Plotting high and low resolution together +if len(variables_with_data_lr) > 0 and len(variables_with_data_hr) > 0: + print(f'We have low and high resolution data, let''s combine it!') + + # Variables with data + variables_with_data = list(variables_with_data_lr) + variables_with_data.extend(x for x in variables_with_data_hr if x not in variables_with_data) + print(f'Variables with data: {variables_with_data}') + + # Variables axis + variables_axis = list(variables_axis_lr) + variables_axis.extend(x for x in variables_axis_hr if x not in variables_axis) + + # Joining dimensions for both resolutions + # Wavelength + wav = list(dimensions_lr['wavelength']) + wav.extend(x for x in dimensions_hr['wavelength'] if x not in wav) + dimensions['wavelength'] = wav + # Add LR/HR to the wavelength + wlr = list(dimensions_lr['wavelength']) + for r in range(len(wlr)): + wlr[r] = str(wlr[r]) + ' (LR)' + whr = list(dimensions_hr['wavelength']) + for r in range(len(whr)): + whr[r] = str(whr[r]) + ' (HR)' + dimensions['wavelength_res'] = wlr + whr + # Time + tim = list(dimensions_lr['time']) + tim.extend(x for x in dimensions_hr['time'] if x not in tim) + dimensions['time'] = tim + # Level + lev = list(dimensions_lr['level']) + lev.extend(x for x in dimensions_hr['level'] if x not in lev) + dimensions['level'] = lev + + # Populate the arrays with data and errors + for v in variables_with_data: + # Create array to store the joint data + # data_var = np.ma.masked_all((len(dimensions['wavelength_res']), len(tim), len(lev))) + data_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev))) + data_error_positive_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev))) + data_error_negative_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev))) + # ToDo check if needed + data_var = np.ma.masked_values(data_var, 0) + data_error_positive_var = np.ma.masked_values(data_error_positive_var, 0) + data_error_negative_var = np.ma.masked_values(data_error_negative_var, 0) + for w in range(len(dimensions['wavelength_res'])): + for t in range(len(dimensions['time'])): + # Incorporate data into the main array + print(f"{v} - Wavelength: {dimensions['wavelength_res'][w]}. Time: {dimensions['time'][t]}.") + + # High resolution + if dimensions['wavelength_res'][w].find('HR') >= 0: + # Is there HR data? + if v in data_hr: + try: + # Find position of w and t + wl = float(dimensions['wavelength_res'][w].split(' ')[0]) + wp = dimensions_hr['wavelength'].index(wl) + tp = dimensions_hr['time'].index(dimensions['time'][t]) + # Check if there is data + a = data_hr[v][wp, tp, :] + d = np.average(np.ma.masked_array(a, np.isnan(a))) + if np.ma.min(d) is not np.ma.masked: + print(f'\tWe have HR data for {v}') + # Save data into common structure + data_var[w, t, :] = a + if np.average(np.ma.masked_array(data_error_positive_hr[v][wp, tp, :], np.isnan( + data_error_positive_hr[v][wp, tp, :]))) is not np.ma.masked: + data_error_positive_var[w, t, :] = data_error_positive_hr[v][wp, tp, :] + if np.average(np.ma.masked_array(data_error_negative_hr[v][wp, tp, :], np.isnan( + data_error_negative_hr[v][wp, tp, :]))) is not np.ma.masked: + data_error_negative_var[w, t, :] = data_error_negative_hr[v][wp, tp, :] + except ValueError: + pass + + # Low resolution + if dimensions['wavelength_res'][w].find('LR') >= 0: + # Is there LR data? + if v in data_lr: + # Find position of w and t + try: + wl = float(dimensions['wavelength_res'][w].split(' ')[0]) + wp = dimensions_hr['wavelength'].index(wl) + if wl in dimensions_lr['wavelength']: + wp = dimensions_lr['wavelength'].index(dimensions['wavelength'][w]) + else: + wp = -1 + # if dimensions['time'][t] in dimensions_lr['time']: + tp = dimensions_lr['time'].index(dimensions['time'][t]) + # else: + # tp = -1 + if wp > -1 and tp > -1: + # Check if there is data + a = data_lr[v][wp, tp, :] + d = np.average(np.ma.masked_array(a, np.isnan(a))) + if np.ma.min(d) is not np.ma.masked: + print(f'\tWe have LR data for {v}') + # Save data into common structure + data_var[w, t, :] = a + if np.average(np.ma.masked_array(data_error_positive_lr[v][wp, tp, :], np.isnan( + data_error_positive_lr[v][wp, tp, :]))) is not np.ma.masked: + data_error_positive_var[w, t, :] = data_error_positive_lr[v][wp, tp, :] + if np.average(np.ma.masked_array(data_error_negative_lr[v][wp, tp, :], np.isnan( + data_error_negative_lr[v][wp, tp, :]))) is not np.ma.masked: + data_error_negative_var[w, t, :] = data_error_negative_lr[v][wp, tp, :] + except ValueError: + # i.e. no 1064 wavelength in dimensions_lr + pass + # except IndexError: + # pass + + # Store data in the global matrices + if variables_with_data.index(v) == 0: + data_all = {v: data_var} + data_error_positive = {v: data_error_positive_var} + data_error_negative = {v: data_error_negative_var} + else: + data_all[v] = data_var + data_error_positive[v] = data_error_positive_var + data_error_negative[v] = data_error_negative_var + + # Plot + timeiter = dimensions['time'] + # Default plotting + dimensions_t = dimensions + for t in range(len(dimensions['time'])): + print(f'plotting hr & lr for {timeiter[t]}') + dimensions_t['time'] = [timeiter[t]] + for vr in range(len(variables_with_data)): + if vr == 0: + data_t = {variables_with_data[vr]: data_all[variables_with_data[vr]][:, t:t + 1, :]} + else: + data_t[variables_with_data[vr]] = data_all[variables_with_data[vr]][:, t:t + 1, :] + + fig = plot_vertical_profiles_plotly_allres( + variables=variables_with_data, + data=data_t, + dimensions=dimensions_t, + y_array=altitude, + positive_error_array=data_error_positive, + negative_error_array=data_error_negative, + xaxis_label=variables_axis, + yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']', + title=station_location + " (" + station_id.upper() + ") - " + str(datetime.fromtimestamp(dimensions_t['time'][0]).strftime('%d/%m/%Y %H:%M:%S')) + " - ELDA MWL - High and low resolution") + + # Show the plot (in Jupyter Notebook or Jupyter Lab) + fig.show() + + # Save the plot + filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp( + dimensions['time'][0]).strftime('%Y%m%d-%H%M%S') + "_all.html" + fig.write_html(html_path + filename) +