Changed the timestamps to UTC and the filenames of the plots default tip

Wed, 28 Aug 2024 11:28:24 +0000

author
pilar <pilar.gumaclaramunt@cnr.it>
date
Wed, 28 Aug 2024 11:28:24 +0000
changeset 4
75213aef70b8
parent 3
bd2181f8295e

Changed the timestamps to UTC and the filenames of the plots

eldamwl_plot.py file | annotate | diff | comparison | revisions
--- a/eldamwl_plot.py	Tue Jul 02 10:29:59 2024 +0000
+++ b/eldamwl_plot.py	Wed Aug 28 11:28:24 2024 +0000
@@ -4,12 +4,22 @@
 import plotly.graph_objects as go
 from plotly.subplots import make_subplots
 import netCDF4 as nc
-from datetime import datetime
+from datetime import datetime, timezone
 import numpy as np
 import sys
 from os.path import exists
 
 
+def check_group_presence(nc_file_path, group_name):
+    with nc.Dataset(nc_file_path, 'r') as dataset:
+        if group_name in dataset.groups:
+            print(f'{group_name} group exists')
+            return True
+        else:
+            print(f'{group_name} group does not exist')
+            return False
+
+
 def read_variable_in_group(nc_file_path, group_name, variable_name):
     try:
         with nc.Dataset(nc_file_path, 'r') as dataset:
@@ -365,6 +375,8 @@
 # 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'
+# v0.0.3
+# pot_012_0000721_202406081503_202406081601_20240608pot150N_ELDAmwl_v0.0.3.nc
 
 
 try:
@@ -405,322 +417,355 @@
 station_location = read_global_attribute(netcdf_file_path, 'location')
 
 # Read data
+h = 0
 for g in range(len(groups)):
-    # Read the altitude variable
-    altitude = read_variable_in_group(netcdf_file_path, groups[g], 'altitude')
+    if check_group_presence(netcdf_file_path, groups[g]):
+        # Read the altitude variable
+        altitude = read_variable_in_group(netcdf_file_path, groups[g], 'altitude')
 
-    if altitude['metadata']['units'] == 'm':
-        altitude['data'] = altitude['data'] / 1000.0
-        altitude['metadata']['units'] = 'km'
+        if altitude['metadata']['units'] == 'm':
+            altitude['data'] = altitude['data'] / 1000.0
+            altitude['metadata']['units'] = 'km'
 
-    # Initialize variables
-    data = None
-    variables_with_data = None
-    variables_axis = None
+        # Initialize variables
+        data = None
+        variables_with_data = None
+        variables_axis = None
 
-    # Read the data
-    for v in range(len(variable)):
-        # Initialize var
-        var = None
+        # Read the data
+        for v in range(len(variable)):
+            # Initialize var
+            var = None
 
-        try:
-            var, dim, met = read_variable_in_group2(netcdf_file_path, groups[g], variable[v])
-        except TypeError:
-            print(f"The variable {variable[v]} is not present. Cannot unpack non-iterable NoneType object")
+            try:
+                var, dim, met = read_variable_in_group2(netcdf_file_path, groups[g], variable[v])
+            except TypeError:
+                print(f"The variable {variable[v]} is not present. Cannot unpack non-iterable NoneType object")
 
-        if var is not None:  # If the variable that we're trying to read exists
-            if var.max() != var.min():  # If there is data in the array
-                if variables_with_data is None:
-                    variables_with_data = [variable[v]]
-                    variables_axis = [met['long_name'] + ' [ ' + met['units'] + ' ]']
-                else:
-                    variables_with_data += [variable[v]]
-                    variables_axis += [met['long_name'] + ' [ ' + met['units'] + ' ]']
+            if var is not None:  # If the variable that we're trying to read exists
+                if var.max() != var.min():  # If there is data in the array
+                    if variables_with_data is None:
+                        variables_with_data = [variable[v]]
+                        variables_axis = [met['long_name'] + ' [ ' + met['units'] + ' ]']
+                    else:
+                        variables_with_data += [variable[v]]
+                        variables_axis += [met['long_name'] + ' [ ' + met['units'] + ' ]']
 
-                var_error, dim_error, met_error = read_variable_in_group2(netcdf_file_path, groups[g],
-                                                                          variable_error[v])
-                if variable[v] == 'volumedepolarization':
-                    var_error_negative, dim, met = read_variable_in_group2(netcdf_file_path, groups[g],
-                                                                           variable_negative_error_volumedepol)
-                else:
-                    var_error_negative = var_error
+                    var_error, dim_error, met_error = read_variable_in_group2(netcdf_file_path, groups[g],
+                                                                              variable_error[v])
+                    if variable[v] == 'volumedepolarization':
+                        var_error_negative, dim, met = read_variable_in_group2(netcdf_file_path, groups[g],
+                                                                               variable_negative_error_volumedepol)
+                    else:
+                        var_error_negative = var_error
 
-                # Add read data to the data variable to have it all grouped
-                if data is None:
-                    # Start populating data and data_error
-                    data = {variable[v]: var}
-                    data_error_positive = {variable[v]: var_error}
-                    data_error_negative = {variable[v]: var_error_negative}
+                    # Add read data to the data variable to have it all grouped
+                    if data is None:
+                        # Start populating data and data_error
+                        data = {variable[v]: var}
+                        data_error_positive = {variable[v]: var_error}
+                        data_error_negative = {variable[v]: var_error_negative}
+                    else:
+                        # Add more data to data and data_error
+                        data[variable[v]] = var
+                        data_error_positive[variable[v]] = var_error
+                        data_error_negative[variable[v]] = var_error_negative
                 else:
-                    # Add more data to data and data_error
-                    data[variable[v]] = var
-                    data_error_positive[variable[v]] = var_error
-                    data_error_negative[variable[v]] = var_error_negative
-            else:
-                print(f'There is no data for {variable[v]} in the {groups_names[g]} group.')
+                    print(f'There is no data for {variable[v]} in the {groups_names[g]} group.')
 
-    dimensions = dim
+        dimensions = dim
 
-    # 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
+        # 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 = True
-        timeinrows = False
-        # 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, "
-                f"do you want to plot them on different plots [y/n]? ")
-            if answer[0:1].lower() == 'n':
-                timeinplots = False
-            if not timeinplots:
+        # 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 h == 0:
+            # Initialize
+            timeinplots = True
+            timeinrows = False
+            # 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, "
-                    f"do you want to plot them on different rows [y/n]? ")
-                if answer[0:1].lower() == 'y':
-                    timeinrows = True
+                    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, "
+                        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:
+                # Default plotting
+                dimensions_t = dimensions
+                for t in range(len(dimensions['time'])):
+                    dimensions_t['time'] = [timeiter[t]]
+                    for vr in range(len(variables_with_data)):
+                        if vr == 0:
+                            data_t = {variables_with_data[vr]: data[variables_with_data[vr]][:, t:t + 1, :]}
+                        else:
+                            data_t[variables_with_data[vr]] = data[variables_with_data[vr]][:, t:t + 1, :]
+
+                    title = station_location + " (" + station_id.upper() + ") - " + str(
+                        datetime.fromtimestamp(dimensions_t['time'][0], tz=timezone.utc).strftime(
+                            '%d/%m/%Y %H:%M:%S')) + " - ELDA MWL " + groups_names[g]
 
-    if variables_with_data is not None:
-        timeiter = dimensions['time']
-        if timeinplots:
-            # Default plotting
-            dimensions_t = dimensions
-            for t in range(len(dimensions['time'])):
-                dimensions_t['time'] = [timeiter[t]]
-                for vr in range(len(variables_with_data)):
-                    if vr == 0:
-                        data_t = {variables_with_data[vr]: data[variables_with_data[vr]][:, t:t + 1, :]}
+                    # Timestamps for the filename
+                    t0 = datetime.fromtimestamp(timeiter[t], tz=timezone.utc).strftime('%Y%m%d%H%M')
+                    if t + 1 < len(timeiter):
+                        t1 = datetime.fromtimestamp(
+                                timeiter[t+1], tz=timezone.utc).strftime('%Y%m%d%H%M')
                     else:
-                        data_t[variables_with_data[vr]] = data[variables_with_data[vr]][:, t:t + 1, :]
+                        t1 = netcdf_file[29:41]
+
+                    filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + t0 + "-" + t1 + "_" + \
+                                   groups[g] + ".html"
 
+                    fig = plot_vertical_profiles_plotly(
+                        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=title,
+                        timeinrows=True
+                    )
+
+                    # Show the plot (in Jupyter Notebook or Jupyter Lab)
+                    # fig.show()
+
+                    # Save the plot
+                    fig.write_html(html_path + filename)
+
+            else:
+                # Plotting in one plot (even if there are different time profiles)
                 fig = plot_vertical_profiles_plotly(
-                    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,
+                    variables_with_data,
+                    data,
+                    dimensions,
+                    altitude,
+                    data_error_positive,
+                    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 " + groups_names[g],
-                    timeinrows=True
+                        datetime.fromtimestamp(dimensions['time'][0], tz=timezone.utc).strftime(
+                            '%d/%m/%Y')) + " - ELDA MWL " + groups_names[
+                              g],
+                    timeinrows=timeinrows
                 )
 
                 # 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') + "_" + groups[g] + ".html"
+                if timeinrows:
+                    filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp(
+                        dimensions['time'][0], tz=timezone.utc).strftime('%Y%m%d%H%M') + '-' + datetime.fromtimestamp(
+                        dimensions['time'][-1], tz=timezone.utc).strftime('%Y%m%d%H%M') + "_" + groups[
+                                   g] + "_timeinrows.html"
+                else:
+                    filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp(
+                        dimensions['time'][0], tz=timezone.utc).strftime('%Y%m%d%H%M') + '-' + datetime.fromtimestamp(
+                        dimensions['time'][-1], tz=timezone.utc).strftime('%Y%m%d%H%M') + "_" + groups[g] + ".html"
+
                 fig.write_html(html_path + filename)
 
+            dimensions['time'] = timeiter
+
         else:
-            # Plotting in one plot (even if there are different time profiles)
-            fig = plot_vertical_profiles_plotly(
-                variables_with_data,
-                data,
-                dimensions,
-                altitude,
-                data_error_positive,
-                data_error_negative,
+            print(f'There is no data to be plotted in the {groups_names[g]} group.')
+
+        h = h + 1
+
+# Plotting high and low resolution together
+if 'variables_with_data_lr' in locals() and 'variables_with_data_hr' in locals():
+    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'])
+        dhrw = list(dimensions_hr['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(dhrw)
+        for r in range(len(whr)):
+            whr[r] = str(whr[r]) + ' (HR)'
+        dimensions['wavelength_res'] = wlr + whr
+        dimensions_hr['wavelength'] = dhrw
+        dimensions['wavelength'] = wav
+        # 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)
+                                wp = dhrw.index(wl)  # ToDo test
+                                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)    # ToDo check
+                                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['time'][0]).strftime('%d/%m/%Y')) + " - ELDA MWL " + groups_names[
-                          g],
-                timeinrows=timeinrows
-            )
+                    datetime.fromtimestamp(dimensions_t['time'][0], tz=timezone.utc).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
-            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_" + 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"
-
-            fig.write_html(html_path + filename)
-
-    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
+            # Timestamps for the filename
+            t0 = datetime.fromtimestamp(timeiter[t], tz=timezone.utc).strftime('%Y%m%d%H%M')
+            if t + 1 < len(timeiter):
+                t1 = datetime.fromtimestamp(
+                    timeiter[t + 1], tz=timezone.utc).strftime('%Y%m%d%H%M')
+            else:
+                t1 = netcdf_file[29:41]
 
-                # 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)
-
+            filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + t0 + "-" + t1 + "_all.html"
+            fig.write_html(html_path + filename)

mercurial