• file
• latest
• revisions
• annotate
• diff
• comparison
• raw

eldamwl_plot.py@763a7ac22031 (annotated)

eldamwl_plot.py

Tue, 02 Jul 2024 08:37:49 +0000

author

pilar <pilar.gumaclaramunt@cnr.it>

date

Tue, 02 Jul 2024 08:37:49 +0000

changeset 2

763a7ac22031

parent 0

fce4cae19357

child 4

75213aef70b8

permissions

-rwxr-xr-x

Added the plotting of high and low resolution together in the same plot.

pilar@2	1	# import nc as nc
pilar@2	2	# import xarray as xr
pilar@0	3	from config import *
pilar@0	4	import plotly.graph_objects as go
pilar@0	5	from plotly.subplots import make_subplots
pilar@0	6	import netCDF4 as nc
pilar@0	7	from datetime import datetime
pilar@0	8	import numpy as np
pilar@0	9	import sys
pilar@0	10	from os.path import exists
pilar@0	11
pilar@0	12
pilar@0	13	def read_variable_in_group(nc_file_path, group_name, variable_name):
pilar@0	14	try:
pilar@0	15	with nc.Dataset(nc_file_path, 'r') as dataset:
pilar@0	16	# Find the group
pilar@0	17	group = dataset.groups[group_name]
pilar@0	18
pilar@0	19	# Access the variable
pilar@0	20	try:
pilar@0	21	variable = group.variables[variable_name]
pilar@0	22	dimensions = variable.dimensions
pilar@0	23	# print("Dimensions:", dimensions)
pilar@0	24	# variable_metadata = { variable.name, variable.long_name, variable.units, variable.ndim}
pilar@0	25
pilar@0	26	# variable_metadata.add = variable.long_name
pilar@0	27	variable_metadata = \
pilar@0	28	{'name': variable.name.capitalize(),
pilar@0	29	'long_name': variable.long_name.capitalize(),
pilar@0	30	'units': variable.units,
pilar@0	31	'ndim': variable.ndim}
pilar@0	32
pilar@0	33	# Create an index array for slicing the variable
pilar@0	34	index_array = [slice(None)] * variable.ndim
pilar@0	35	# Access the variable data using the index array
pilar@0	36	variable_data = variable[index_array]
pilar@0	37
pilar@0	38	# Fetch the dimension values
pilar@0	39	dimension_values = {}
pilar@0	40	for dim_name in dimensions:
pilar@0	41	dim_values = group.variables[dim_name][:]
pilar@0	42	dimension_values[dim_name] = dim_values
pilar@0	43	# Adjust the index array for slicing
pilar@0	44	index_array[variable.dimensions.index(dim_name)] = 0
pilar@0	45
pilar@0	46	# Now you have the variable data and its associated dimensions in a dictionary
pilar@0	47	result = {'data': variable_data, 'dimensions': dimension_values, 'metadata': variable_metadata}
pilar@0	48
pilar@0	49	# Print or use the result as needed
pilar@0	50	# print(result)
pilar@0	51	return result
pilar@0	52
pilar@0	53	except KeyError:
pilar@0	54	print("The specified variable does not exist in the group.")
pilar@0	55
pilar@0	56	except FileNotFoundError:
pilar@0	57	print(f"The file '{nc_file_path}' does not exist.")
pilar@0	58	except Exception as e:
pilar@0	59	print(f"The group {e} does not exist.")
pilar@0	60
pilar@0	61
pilar@2	62	def read_global_attribute(nc_file_path, attribute_name):
pilar@2	63	# Open the NetCDF file in read mode
pilar@2	64	with nc.Dataset(nc_file_path, 'r') as dataset:
pilar@2	65	# with netCDF4.Dataset(nc_file_path, 'r') as nc:
pilar@2	66	# Check if the attribute exists
pilar@2	67	if attribute_name in dataset.ncattrs():
pilar@2	68	# Read the attribute value
pilar@2	69	attribute_value = dataset.getncattr(attribute_name)
pilar@2	70	return attribute_value
pilar@2	71	else:
pilar@2	72	print(f"Attribute '{attribute_name}' not found in the NetCDF file.")
pilar@2	73
pilar@2	74
pilar@0	75	def read_variable_in_group2(nc_file_path, group_name, variable_name):
pilar@0	76	try:
pilar@0	77	with nc.Dataset(nc_file_path, 'r') as dataset:
pilar@0	78	# Find the group
pilar@0	79	group = dataset.groups[group_name]
pilar@0	80
pilar@0	81	# Access the variable
pilar@0	82	try:
pilar@0	83	variable = group.variables[variable_name]
pilar@0	84	dimensions = variable.dimensions
pilar@0	85	# print("Dimensions:", dimensions)
pilar@0	86	# variable_metadata = { variable.name, variable.long_name, variable.units, variable.ndim}
pilar@0	87
pilar@0	88	# variable_metadata.add = variable.long_name
pilar@0	89	variable_metadata = \
pilar@0	90	{'name': variable.name.capitalize(),
pilar@0	91	'long_name': variable.long_name.capitalize(),
pilar@0	92	'units': variable.units,
pilar@0	93	'ndim': variable.ndim}
pilar@0	94
pilar@0	95	# Create an index array for slicing the variable
pilar@0	96	index_array = [slice(None)] * variable.ndim
pilar@0	97	# Access the variable data using the index array
pilar@0	98	variable_data = variable[index_array]
pilar@0	99
pilar@0	100	# Fetch the dimension values
pilar@0	101	dimension_values = {}
pilar@0	102	for dim_name in dimensions:
pilar@0	103	dim_values = group.variables[dim_name][:]
pilar@2	104	dimension_values[dim_name] = dim_values.tolist()
pilar@0	105	# Adjust the index array for slicing
pilar@0	106	index_array[variable.dimensions.index(dim_name)] = 0
pilar@0	107
pilar@0	108	# Now you have the variable data and its associated dimensions in a dictionary
pilar@0	109	# result = {'data': variable_data, 'dimensions': dimension_values, 'metadata': variable_metadata}
pilar@0	110
pilar@0	111	# Print or use the result as needed
pilar@0	112	# print(result)
pilar@0	113	return variable_data, dimension_values, variable_metadata
pilar@0	114
pilar@0	115	except KeyError:
pilar@0	116	print(f"The variable {variable_name} does not exist in the {group_name} group.")
pilar@0	117	except FileNotFoundError:
pilar@0	118	print(f"The file {nc_file_path} does not exist.")
pilar@0	119	except Exception as e:
pilar@0	120	print(f"The group {e} does not exist.")
pilar@0	121
pilar@0	122
pilar@0	123	def plot_vertical_profiles_plotly(variables, data, dimensions, y_array, positive_error_array, negative_error_array,
pilar@0	124	xaxis_label, yaxis_label, title, timeinrows):
pilar@0	125	"""
pilar@0	126	ToDo Complete
pilar@0	127	Plot vertical profiles using plotly
pilar@0	128
pilar@0	129	Parameters:
pilar@0	130	-
pilar@0	131
pilar@0	132	Returns:
pilar@0	133	- fig: html
pilar@0	134	"""
pilar@0	135
pilar@0	136	# Set number of columns, number of rows, and subtitles if needed
pilar@0	137	ncols = len(variables)
pilar@0	138	subtitles = None
pilar@0	139	if timeinrows is True:
pilar@0	140	nrows = len(dimensions['time'])
pilar@0	141	if nrows > 1:
pilar@0	142	for t in range(nrows):
pilar@0	143	for v in range(ncols):
pilar@0	144	if subtitles is None:
pilar@0	145	subtitles = [str(datetime.fromtimestamp(dimensions['time'][t]))]
pilar@0	146	else:
pilar@0	147	subtitles += [str(datetime.fromtimestamp(dimensions['time'][t]))]
pilar@0	148	else:
pilar@0	149	nrows = 1
pilar@0	150	subtitles = ''
pilar@0	151
pilar@0	152	# Create figure
pilar@0	153	fig = make_subplots(rows=nrows, cols=ncols,
pilar@0	154	subplot_titles=subtitles,
pilar@0	155	shared_yaxes=True)
pilar@0	156
pilar@0	157	# Define colours
pilar@0	158	clrs = [['rgb(0,35,255)', 'rgb(0,20,150)', 'rgb(0,170,250)'], # UV
pilar@0	159	['rgb(90,210,50)', 'rgb(40,90,25)', 'rgb(120,250,80)'], # VIS
pilar@0	160	['rgb(250,30,30)', 'rgb(200,25,30)', 'rgb(250,125,0)']] # IR
pilar@0	161	clrs = np.array(clrs)
pilar@0	162
pilar@0	163	# Plot
pilar@0	164	for v in range(len(variables)):
pilar@0	165	for t in range(len(dimensions['time'])):
pilar@0	166	for w in range(len(dimensions['wavelength'])):
pilar@0	167	tn = str(datetime.fromtimestamp(dimensions['time'][t]))
pilar@0	168	wn = str(dimensions['wavelength'][w])
pilar@2	169	if timeinrows:
pilar@0	170	lgnd = wn + ' nm'
pilar@0	171	rown = t + 1
pilar@0	172	colorind = 0
pilar@2	173	# subtitle = str(datetime.fromtimestamp(dimensions['time'][t]))
pilar@0	174	else:
pilar@0	175	lgnd = wn + ' nm - ' + tn
pilar@0	176	rown = 1
pilar@0	177	colorind = t
pilar@0	178
pilar@0	179	# colour definition
pilar@0	180	if abs(dimensions['wavelength'][w] - 1064) < 1:
pilar@0	181	clr = clrs[2, colorind]
pilar@0	182	# clr = 'rgb(235,70,20)' # 'rgba(215,27,96,0.4)'
pilar@0	183	elif abs(dimensions['wavelength'][w] - 532) < 1:
pilar@0	184	clr = clrs[1, colorind]
pilar@0	185	elif abs(dimensions['wavelength'][w] - 355) < 1:
pilar@0	186	clr = clrs[0, colorind]
pilar@0	187	else:
pilar@0	188	clr = 'rgb(0,0,0)'
pilar@0	189
pilar@0	190	data_values = data[variables[v]][w, t, :]
pilar@2	191	if np.average(np.ma.masked_array(data_values, np.isnan(data_values))) is not np.ma.masked:
pilar@2	192	vertical_levels = y_array["data"][t, :]
pilar@2	193	error_values = positive_error_array[variables[v]][w, t, :]
pilar@2	194
pilar@2	195	# Add trace
pilar@2	196	fig.add_trace(go.Scatter(
pilar@2	197	x=data_values,
pilar@2	198	y=vertical_levels,
pilar@2	199	error_x=dict(array=error_values, width=0, thickness=0.1),
pilar@2	200	# ToDo Include distinction between positive and negative errors
pilar@2	201	mode='lines',
pilar@2	202	line=dict(width=2, color=clr),
pilar@2	203	name=lgnd),
pilar@2	204	row=rown, col=v + 1,
pilar@2	205	)
pilar@2	206
pilar@2	207	# ToDo Check formatting options (showlegend, mapbox)
pilar@2	208	fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1)
pilar@2	209	fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown,
pilar@2	210	col=v + 1)
pilar@2	211
pilar@2	212	fig.update_xaxes(ticks="outside")
pilar@2	213	fig.update_yaxes(ticks="outside", row=rown, col=v + 1)
pilar@2	214	# ToDo check tickformat (i.e. tickformat=":.2e")
pilar@2	215	if v == 0:
pilar@2	216	fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1)
pilar@2	217
pilar@2	218	# ToDo add boxes
pilar@2	219	a = fig.layout['xaxis']['domain']
pilar@2	220	# print(a)
pilar@2	221	b = fig.layout['yaxis']['domain']
pilar@2	222	# print(b)
pilar@2	223	# fig.add_hline(y=b[0], line=dict(color="yellow", width=2), row=rown, col=v+1)
pilar@2	224	# fig.add_hline(y=b[1], line=dict(color="purple", width=2), row=rown, col=v+1)
pilar@2	225	# fig.add_vline(x=a[0], line=dict(color="orange", width=2), row=rown, col=v+1)
pilar@2	226	# fig.add_vline(x=a[1], line=dict(color="pink", width=2), row=rown, col=v+1)
pilar@2	227	# fig.add_shape(type="line", xref="paper", yref="paper", x0=a[0], y0=b[0], x1=a[1], y1=b[1],
pilar@2	228	# line=dict(color="grey", width=2, ), row=rown, col=v+1)
pilar@2	229	# fig.add_shape(type="line", xref="paper", yref="paper", x0=1, y0=1, x1=1, y1=0,
pilar@2	230	# line=dict(color="blue", width=2, ), row=rown, col=v+1)
pilar@2	231	# fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=1, y1=0,
pilar@2	232	# line=dict(color="red", width=2, ), row=rown, col=v+1)
pilar@2	233	# fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=0, y1=1,
pilar@2	234	# line=dict(color="green", width=2, ), row=rown, col=v+1)
pilar@2	235	# fig.add_hline(y=0, line=dict(color="yellow", width=2))
pilar@2	236	# print(fig.data[0]['x'])
pilar@2	237
pilar@2	238	# Set axis labels and title
pilar@2	239	fig.update_layout(
pilar@2	240	yaxis_title=yaxis_label,
pilar@2	241	title=title,
pilar@2	242	template="seaborn", # You can change the template as per your preference
pilar@2	243	# ['ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation',
pilar@2	244	# 'xgridoff', 'ygridoff', 'gridon', 'none']
pilar@2	245	height=max(600, 450 * rown), width=min(max(400 * len(variables), 750), 1400),
pilar@2	246	)
pilar@2	247
pilar@2	248	return fig
pilar@2	249
pilar@2	250
pilar@2	251	def plot_vertical_profiles_plotly_allres(variables, data, dimensions, y_array, positive_error_array,
pilar@2	252	negative_error_array, xaxis_label, yaxis_label, title):
pilar@2	253	"""
pilar@2	254	ToDo Complete
pilar@2	255	Plot vertical profiles using plotly
pilar@2	256
pilar@2	257	Parameters:
pilar@2	258	-
pilar@2	259
pilar@2	260	Returns:
pilar@2	261	- fig: html
pilar@2	262	"""
pilar@2	263
pilar@2	264	# Set number of columns, number of rows, and subtitles if needed
pilar@2	265	ncols = len(variables)
pilar@2	266	nrows = 1
pilar@2	267	subtitles = ''
pilar@0	268
pilar@2	269	# Create figure
pilar@2	270	fig = make_subplots(rows=nrows, cols=ncols,
pilar@2	271	subplot_titles=subtitles,
pilar@2	272	shared_yaxes=True)
pilar@2	273
pilar@2	274	# Define colours
pilar@2	275	clrs = [['rgb(0,35,255)', 'rgb(0,20,150)'], # UV
pilar@2	276	['rgb(90,210,50)', 'rgb(40,90,25)'], # VIS
pilar@2	277	['rgb(250,30,30)', 'rgb(130,30,10)']] # IR
pilar@2	278	clrs = np.array(clrs)
pilar@2	279
pilar@2	280	# Plot
pilar@2	281	for v in range(len(variables)):
pilar@2	282	for t in range(len(dimensions['time'])):
pilar@2	283	for w in range(len(dimensions['wavelength_res'])):
pilar@2	284	wl = float(dimensions['wavelength_res'][w].split(' ')[0])
pilar@2	285	wn = str(dimensions['wavelength_res'][w])
pilar@2	286	lgnd = wn.replace(' ', ' nm ')
pilar@2	287	rown = 1
pilar@2	288	if wn.find('LR') > 0:
pilar@2	289	colorind = 0
pilar@2	290	else:
pilar@2	291	colorind = 1
pilar@2	292
pilar@2	293	# colour definition
pilar@2	294	if abs(wl - 1064) < 1:
pilar@2	295	clr = clrs[2, colorind]
pilar@2	296	elif abs(wl - 532) < 1:
pilar@2	297	clr = clrs[1, colorind]
pilar@2	298	elif abs(wl - 355) < 1:
pilar@2	299	clr = clrs[0, colorind]
pilar@2	300	else:
pilar@2	301	clr = 'rgb(0,0,0)'
pilar@0	302
pilar@2	303	data_values = data[variables[v]][w, t, :]
pilar@2	304	if np.average(np.ma.masked_array(data_values, np.isnan(data_values))) is not np.ma.masked:
pilar@2	305	vertical_levels = y_array["data"][t, :]
pilar@2	306	error_values = positive_error_array[variables[v]][w, t, :]
pilar@2	307	# Add trace
pilar@2	308	fig.add_trace(go.Scatter(
pilar@2	309	x=data_values,
pilar@2	310	y=vertical_levels,
pilar@2	311	error_x=dict(array=error_values, width=0, thickness=0.1),
pilar@2	312	# ToDo Include distinction between positive and negative errors
pilar@2	313	mode='lines',
pilar@2	314	line=dict(width=2, color=clr),
pilar@2	315	name=lgnd),
pilar@2	316	row=rown, col=v + 1,
pilar@2	317	)
pilar@2	318
pilar@2	319	# ToDo Check formatting options (showlegend, mapbox)
pilar@2	320	fig.update_xaxes({"title": xaxis_label[v]}, row=rown, col=v + 1)
pilar@2	321	fig.update_xaxes({"mirror": True, "ticks": 'outside', "showline": True, "color": "black"}, row=rown,
pilar@2	322	col=v + 1)
pilar@0	323
pilar@2	324	fig.update_xaxes(ticks="outside")
pilar@2	325	fig.update_yaxes(ticks="outside", row=rown, col=v + 1)
pilar@2	326	# ToDo check tickformat (i.e. tickformat=":.2e")
pilar@2	327	if v == 0:
pilar@2	328	fig.update_yaxes({"title": yaxis_label}, row=rown, col=v + 1)
pilar@2	329
pilar@2	330	# ToDo add boxes
pilar@2	331	a = fig.layout['xaxis']['domain']
pilar@2	332	# print(a)
pilar@2	333	b = fig.layout['yaxis']['domain']
pilar@2	334	# print(b)
pilar@2	335	# fig.add_hline(y=b[0], line=dict(color="yellow", width=2), row=rown, col=v+1)
pilar@2	336	# fig.add_hline(y=b[1], line=dict(color="purple", width=2), row=rown, col=v+1)
pilar@2	337	# fig.add_vline(x=a[0], line=dict(color="orange", width=2), row=rown, col=v+1)
pilar@2	338	# fig.add_vline(x=a[1], line=dict(color="pink", width=2), row=rown, col=v+1)
pilar@2	339	# fig.add_shape(type="line", xref="paper", yref="paper", x0=a[0], y0=b[0], x1=a[1], y1=b[1],
pilar@2	340	# line=dict(color="grey", width=2, ), row=rown, col=v+1)
pilar@2	341	# fig.add_shape(type="line", xref="paper", yref="paper", x0=1, y0=1, x1=1, y1=0,
pilar@2	342	# line=dict(color="blue", width=2, ), row=rown, col=v+1)
pilar@2	343	# fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=1, y1=0,
pilar@2	344	# line=dict(color="red", width=2, ), row=rown, col=v+1)
pilar@2	345	# fig.add_shape(type="line", xref="paper", yref="paper", x0=0, y0=0, x1=0, y1=1,
pilar@2	346	# line=dict(color="green", width=2, ), row=rown, col=v+1)
pilar@2	347	# fig.add_hline(y=0, line=dict(color="yellow", width=2))
pilar@2	348	# print(fig.data[0]['x'])
pilar@0	349
pilar@0	350	# Set axis labels and title
pilar@0	351	fig.update_layout(
pilar@0	352	yaxis_title=yaxis_label,
pilar@0	353	title=title,
pilar@0	354	template="seaborn", # You can change the template as per your preference
pilar@0	355	# ['ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation',
pilar@0	356	# 'xgridoff', 'ygridoff', 'gridon', 'none']
pilar@0	357	height=max(600, 450 * rown), width=min(max(400 * len(variables), 750), 1400),
pilar@0	358	)
pilar@0	359
pilar@0	360	return fig
pilar@0	361
pilar@0	362
pilar@0	363	# Files for testing
pilar@2	364	# Backscatter and extinction, no error:
pilar@2	365	# netcdf_file = 'hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc'
pilar@2	366	# Extinction with error:
pilar@2	367	# netcdf_file = 'hpb_012_0000627_202308240200_202308240400_20230824hpb0200_ELDAmwl_v0.0.1.nc'
pilar@0	368
pilar@0	369
pilar@0	370	try:
pilar@0	371	netcdf_file = sys.argv[1]
pilar@0	372	except:
pilar@0	373	print('Syntax is incorrect, you need to pass the filename of the file to be plotted.\n'
pilar@0	374	'Example: ./eldamwl_plot.py hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc')
pilar@0	375	exit()
pilar@0	376
pilar@2	377	# Print help
pilar@2	378	if netcdf_file == '-h' or netcdf_file == '--h':
pilar@2	379	print('Syntax: ./eldamwl_plot.py hpb_012_0000598_201810172100_201810172300_20181017oh00_ELDAmwl_v0.0.1.nc')
pilar@2	380	exit()
pilar@2	381
pilar@0	382	netcdf_file_path = netcdf_path + netcdf_file
pilar@0	383
pilar@0	384	file_exists = exists(netcdf_file_path)
pilar@0	385	if not file_exists:
pilar@0	386	print('The file that you provided does not exist.')
pilar@0	387	exit()
pilar@0	388
pilar@0	389	# Groups of data in the NetCDF file (low and high resolution)
pilar@0	390	groups = ['lowres_products', 'highres_products']
pilar@0	391	groups_names = ['Low resolution', 'High resolution']
pilar@0	392
pilar@0	393	# Variables that we want to plot
pilar@2	394	variable = ['backscatter', 'extinction', 'lidarratio', 'volumedepolarization',
pilar@2	395	'particledepolarization'] # ToDo Test particledepolarization when available
pilar@0	396	# Errors associated to the variables
pilar@0	397	variable_error = ['error_backscatter', 'error_extinction', 'error_lidarratio',
pilar@0	398	'positive_systematic_error_volumedepolarization', 'error_particledepolarization']
pilar@0	399	variable_negative_error_volumedepol = 'negative_systematic_error_volumedepolarization'
pilar@2	400	# ToDo Consider adding other variables:
pilar@0	401	# cloud_mask (time, level)
pilar@0	402
pilar@2	403	# Read global attributes
pilar@2	404	station_id = read_global_attribute(netcdf_file_path, 'station_ID')
pilar@2	405	station_location = read_global_attribute(netcdf_file_path, 'location')
pilar@0	406
pilar@2	407	# Read data
pilar@0	408	for g in range(len(groups)):
pilar@0	409	# Read the altitude variable
pilar@0	410	altitude = read_variable_in_group(netcdf_file_path, groups[g], 'altitude')
pilar@0	411
pilar@0	412	if altitude['metadata']['units'] == 'm':
pilar@0	413	altitude['data'] = altitude['data'] / 1000.0
pilar@0	414	altitude['metadata']['units'] = 'km'
pilar@0	415
pilar@0	416	# Initialize variables
pilar@0	417	data = None
pilar@0	418	variables_with_data = None
pilar@0	419	variables_axis = None
pilar@0	420
pilar@0	421	# Read the data
pilar@0	422	for v in range(len(variable)):
pilar@0	423	# Initialize var
pilar@0	424	var = None
pilar@0	425
pilar@0	426	try:
pilar@0	427	var, dim, met = read_variable_in_group2(netcdf_file_path, groups[g], variable[v])
pilar@0	428	except TypeError:
pilar@0	429	print(f"The variable {variable[v]} is not present. Cannot unpack non-iterable NoneType object")
pilar@0	430
pilar@0	431	if var is not None: # If the variable that we're trying to read exists
pilar@0	432	if var.max() != var.min(): # If there is data in the array
pilar@0	433	if variables_with_data is None:
pilar@0	434	variables_with_data = [variable[v]]
pilar@0	435	variables_axis = [met['long_name'] + ' [ ' + met['units'] + ' ]']
pilar@0	436	else:
pilar@0	437	variables_with_data += [variable[v]]
pilar@0	438	variables_axis += [met['long_name'] + ' [ ' + met['units'] + ' ]']
pilar@0	439
pilar@0	440	var_error, dim_error, met_error = read_variable_in_group2(netcdf_file_path, groups[g],
pilar@0	441	variable_error[v])
pilar@0	442	if variable[v] == 'volumedepolarization':
pilar@0	443	var_error_negative, dim, met = read_variable_in_group2(netcdf_file_path, groups[g],
pilar@0	444	variable_negative_error_volumedepol)
pilar@0	445	else:
pilar@0	446	var_error_negative = var_error
pilar@0	447
pilar@0	448	# Add read data to the data variable to have it all grouped
pilar@0	449	if data is None:
pilar@0	450	# Start populating data and data_error
pilar@0	451	data = {variable[v]: var}
pilar@0	452	data_error_positive = {variable[v]: var_error}
pilar@0	453	data_error_negative = {variable[v]: var_error_negative}
pilar@0	454	else:
pilar@0	455	# Add more data to data and data_error
pilar@0	456	data[variable[v]] = var
pilar@0	457	data_error_positive[variable[v]] = var_error
pilar@0	458	data_error_negative[variable[v]] = var_error_negative
pilar@0	459	else:
pilar@0	460	print(f'There is no data for {variable[v]} in the {groups_names[g]} group.')
pilar@0	461
pilar@0	462	dimensions = dim
pilar@0	463
pilar@2	464	# Save lowres and highres in variables to plot them together afterwards
pilar@2	465	if groups[g] == 'lowres_products':
pilar@2	466	altitude_lr = altitude
pilar@2	467	dimensions_lr = dim
pilar@2	468	data_lr = data
pilar@2	469	data_error_positive_lr = data_error_positive
pilar@2	470	data_error_negative_lr = data_error_negative
pilar@2	471	variables_with_data_lr = variables_with_data
pilar@2	472	variables_axis_lr = variables_axis
pilar@2	473	if groups[g] == 'highres_products':
pilar@2	474	altitude_hr = altitude
pilar@2	475	dimensions_hr = dim
pilar@2	476	data_hr = data
pilar@2	477	data_error_positive_hr = data_error_positive
pilar@2	478	data_error_negative_hr = data_error_negative
pilar@2	479	variables_with_data_hr = variables_with_data
pilar@2	480	variables_axis_hr = variables_axis
pilar@2	481
pilar@2	482	# If interactive is set to True, ask user if he/she wants to plot different temporal profiles in different rows
pilar@2	483	# If interactive is set to False, proceed with default values
pilar@0	484	if g == 0:
pilar@0	485	# Initialize
pilar@2	486	timeinplots = True
pilar@0	487	timeinrows = False
pilar@2	488	# If interactive (config.py) = True, ask user how he/she wants the plot
pilar@2	489	if interactive and len(dimensions['time']) > 1:
pilar@0	490	answer = input(
pilar@2	491	f"There are {len(dimensions['time'])} temporal profiles, "
pilar@2	492	f"do you want to plot them on different plots [y/n]? ")
pilar@2	493	if answer[0:1].lower() == 'n':
pilar@2	494	timeinplots = False
pilar@2	495	if not timeinplots:
pilar@0	496	answer = input(
pilar@2	497	f"There are {len(dimensions['time'])} temporal profiles, "
pilar@2	498	f"do you want to plot them on different rows [y/n]? ")
pilar@0	499	if answer[0:1].lower() == 'y':
pilar@0	500	timeinrows = True
pilar@0	501
pilar@0	502	if variables_with_data is not None:
pilar@0	503	timeiter = dimensions['time']
pilar@2	504	if timeinplots:
pilar@2	505	# Default plotting
pilar@0	506	dimensions_t = dimensions
pilar@0	507	for t in range(len(dimensions['time'])):
pilar@0	508	dimensions_t['time'] = [timeiter[t]]
pilar@0	509	for vr in range(len(variables_with_data)):
pilar@0	510	if vr == 0:
pilar@0	511	data_t = {variables_with_data[vr]: data[variables_with_data[vr]][:, t:t + 1, :]}
pilar@0	512	else:
pilar@0	513	data_t[variables_with_data[vr]] = data[variables_with_data[vr]][:, t:t + 1, :]
pilar@0	514
pilar@0	515	fig = plot_vertical_profiles_plotly(
pilar@0	516	variables=variables_with_data,
pilar@0	517	data=data_t,
pilar@0	518	dimensions=dimensions_t,
pilar@0	519	y_array=altitude,
pilar@0	520	positive_error_array=data_error_positive,
pilar@0	521	negative_error_array=data_error_negative,
pilar@0	522	xaxis_label=variables_axis,
pilar@0	523	yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']',
pilar@2	524	title=station_location + " (" + station_id.upper() + ") - " + str(
pilar@2	525	datetime.fromtimestamp(dimensions_t['time'][0]).strftime('%d/%m/%Y %H:%M:%S')) +
pilar@2	526	" - ELDA MWL " + groups_names[g],
pilar@0	527	timeinrows=True
pilar@0	528	)
pilar@0	529
pilar@0	530	# Show the plot (in Jupyter Notebook or Jupyter Lab)
pilar@2	531	# fig.show()
pilar@0	532
pilar@0	533	# Save the plot
pilar@0	534	filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp(
pilar@0	535	dimensions['time'][0]).strftime('%Y%m%d-%H%M%S') + "_" + groups[g] + ".html"
pilar@0	536	fig.write_html(html_path + filename)
pilar@0	537
pilar@0	538	else:
pilar@0	539	# Plotting in one plot (even if there are different time profiles)
pilar@0	540	fig = plot_vertical_profiles_plotly(
pilar@0	541	variables_with_data,
pilar@0	542	data,
pilar@0	543	dimensions,
pilar@0	544	altitude,
pilar@0	545	data_error_positive,
pilar@0	546	data_error_negative,
pilar@0	547	xaxis_label=variables_axis,
pilar@0	548	yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']',
pilar@2	549	title=station_location + " (" + station_id.upper() + ") - " + str(
pilar@0	550	datetime.fromtimestamp(dimensions['time'][0]).strftime('%d/%m/%Y')) + " - ELDA MWL " + groups_names[
pilar@0	551	g],
pilar@0	552	timeinrows=timeinrows
pilar@0	553	)
pilar@0	554
pilar@0	555	# Show the plot (in Jupyter Notebook or Jupyter Lab)
pilar@2	556	# fig.show()
pilar@0	557
pilar@0	558	# Save the plot
pilar@2	559	if timeinrows:
pilar@2	560	filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp(
pilar@0	561	dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp(
pilar@0	562	dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + "_timeinrows.html"
pilar@0	563	else:
pilar@2	564	filename = "ELDAmwl_" + station_id.upper() + "_" + datetime.fromtimestamp(
pilar@0	565	dimensions['time'][0]).strftime('%Y%m%d%H%M%S') + '-' + datetime.fromtimestamp(
pilar@0	566	dimensions['time'][-1]).strftime('%Y%m%d%H%M%S') + "_" + groups[g] + ".html"
pilar@0	567
pilar@0	568	fig.write_html(html_path + filename)
pilar@0	569
pilar@0	570	else:
pilar@0	571	print(f'There is no data to be plotted in the {groups_names[g]} group.')
pilar@2	572
pilar@2	573	# Plotting high and low resolution together
pilar@2	574	if len(variables_with_data_lr) > 0 and len(variables_with_data_hr) > 0:
pilar@2	575	print(f'We have low and high resolution data, let''s combine it!')
pilar@2	576
pilar@2	577	# Variables with data
pilar@2	578	variables_with_data = list(variables_with_data_lr)
pilar@2	579	variables_with_data.extend(x for x in variables_with_data_hr if x not in variables_with_data)
pilar@2	580	print(f'Variables with data: {variables_with_data}')
pilar@2	581
pilar@2	582	# Variables axis
pilar@2	583	variables_axis = list(variables_axis_lr)
pilar@2	584	variables_axis.extend(x for x in variables_axis_hr if x not in variables_axis)
pilar@2	585
pilar@2	586	# Joining dimensions for both resolutions
pilar@2	587	# Wavelength
pilar@2	588	wav = list(dimensions_lr['wavelength'])
pilar@2	589	wav.extend(x for x in dimensions_hr['wavelength'] if x not in wav)
pilar@2	590	dimensions['wavelength'] = wav
pilar@2	591	# Add LR/HR to the wavelength
pilar@2	592	wlr = list(dimensions_lr['wavelength'])
pilar@2	593	for r in range(len(wlr)):
pilar@2	594	wlr[r] = str(wlr[r]) + ' (LR)'
pilar@2	595	whr = list(dimensions_hr['wavelength'])
pilar@2	596	for r in range(len(whr)):
pilar@2	597	whr[r] = str(whr[r]) + ' (HR)'
pilar@2	598	dimensions['wavelength_res'] = wlr + whr
pilar@2	599	# Time
pilar@2	600	tim = list(dimensions_lr['time'])
pilar@2	601	tim.extend(x for x in dimensions_hr['time'] if x not in tim)
pilar@2	602	dimensions['time'] = tim
pilar@2	603	# Level
pilar@2	604	lev = list(dimensions_lr['level'])
pilar@2	605	lev.extend(x for x in dimensions_hr['level'] if x not in lev)
pilar@2	606	dimensions['level'] = lev
pilar@2	607
pilar@2	608	# Populate the arrays with data and errors
pilar@2	609	for v in variables_with_data:
pilar@2	610	# Create array to store the joint data
pilar@2	611	# data_var = np.ma.masked_all((len(dimensions['wavelength_res']), len(tim), len(lev)))
pilar@2	612	data_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev)))
pilar@2	613	data_error_positive_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev)))
pilar@2	614	data_error_negative_var = np.ma.zeros((len(dimensions['wavelength_res']), len(tim), len(lev)))
pilar@2	615	# ToDo check if needed
pilar@2	616	data_var = np.ma.masked_values(data_var, 0)
pilar@2	617	data_error_positive_var = np.ma.masked_values(data_error_positive_var, 0)
pilar@2	618	data_error_negative_var = np.ma.masked_values(data_error_negative_var, 0)
pilar@2	619	for w in range(len(dimensions['wavelength_res'])):
pilar@2	620	for t in range(len(dimensions['time'])):
pilar@2	621	# Incorporate data into the main array
pilar@2	622	print(f"{v} - Wavelength: {dimensions['wavelength_res'][w]}. Time: {dimensions['time'][t]}.")
pilar@2	623
pilar@2	624	# High resolution
pilar@2	625	if dimensions['wavelength_res'][w].find('HR') >= 0:
pilar@2	626	# Is there HR data?
pilar@2	627	if v in data_hr:
pilar@2	628	try:
pilar@2	629	# Find position of w and t
pilar@2	630	wl = float(dimensions['wavelength_res'][w].split(' ')[0])
pilar@2	631	wp = dimensions_hr['wavelength'].index(wl)
pilar@2	632	tp = dimensions_hr['time'].index(dimensions['time'][t])
pilar@2	633	# Check if there is data
pilar@2	634	a = data_hr[v][wp, tp, :]
pilar@2	635	d = np.average(np.ma.masked_array(a, np.isnan(a)))
pilar@2	636	if np.ma.min(d) is not np.ma.masked:
pilar@2	637	print(f'\tWe have HR data for {v}')
pilar@2	638	# Save data into common structure
pilar@2	639	data_var[w, t, :] = a
pilar@2	640	if np.average(np.ma.masked_array(data_error_positive_hr[v][wp, tp, :], np.isnan(
pilar@2	641	data_error_positive_hr[v][wp, tp, :]))) is not np.ma.masked:
pilar@2	642	data_error_positive_var[w, t, :] = data_error_positive_hr[v][wp, tp, :]
pilar@2	643	if np.average(np.ma.masked_array(data_error_negative_hr[v][wp, tp, :], np.isnan(
pilar@2	644	data_error_negative_hr[v][wp, tp, :]))) is not np.ma.masked:
pilar@2	645	data_error_negative_var[w, t, :] = data_error_negative_hr[v][wp, tp, :]
pilar@2	646	except ValueError:
pilar@2	647	pass
pilar@2	648
pilar@2	649	# Low resolution
pilar@2	650	if dimensions['wavelength_res'][w].find('LR') >= 0:
pilar@2	651	# Is there LR data?
pilar@2	652	if v in data_lr:
pilar@2	653	# Find position of w and t
pilar@2	654	try:
pilar@2	655	wl = float(dimensions['wavelength_res'][w].split(' ')[0])
pilar@2	656	wp = dimensions_hr['wavelength'].index(wl)
pilar@2	657	if wl in dimensions_lr['wavelength']:
pilar@2	658	wp = dimensions_lr['wavelength'].index(dimensions['wavelength'][w])
pilar@2	659	else:
pilar@2	660	wp = -1
pilar@2	661	# if dimensions['time'][t] in dimensions_lr['time']:
pilar@2	662	tp = dimensions_lr['time'].index(dimensions['time'][t])
pilar@2	663	# else:
pilar@2	664	# tp = -1
pilar@2	665	if wp > -1 and tp > -1:
pilar@2	666	# Check if there is data
pilar@2	667	a = data_lr[v][wp, tp, :]
pilar@2	668	d = np.average(np.ma.masked_array(a, np.isnan(a)))
pilar@2	669	if np.ma.min(d) is not np.ma.masked:
pilar@2	670	print(f'\tWe have LR data for {v}')
pilar@2	671	# Save data into common structure
pilar@2	672	data_var[w, t, :] = a
pilar@2	673	if np.average(np.ma.masked_array(data_error_positive_lr[v][wp, tp, :], np.isnan(
pilar@2	674	data_error_positive_lr[v][wp, tp, :]))) is not np.ma.masked:
pilar@2	675	data_error_positive_var[w, t, :] = data_error_positive_lr[v][wp, tp, :]
pilar@2	676	if np.average(np.ma.masked_array(data_error_negative_lr[v][wp, tp, :], np.isnan(
pilar@2	677	data_error_negative_lr[v][wp, tp, :]))) is not np.ma.masked:
pilar@2	678	data_error_negative_var[w, t, :] = data_error_negative_lr[v][wp, tp, :]
pilar@2	679	except ValueError:
pilar@2	680	# i.e. no 1064 wavelength in dimensions_lr
pilar@2	681	pass
pilar@2	682	# except IndexError:
pilar@2	683	# pass
pilar@2	684
pilar@2	685	# Store data in the global matrices
pilar@2	686	if variables_with_data.index(v) == 0:
pilar@2	687	data_all = {v: data_var}
pilar@2	688	data_error_positive = {v: data_error_positive_var}
pilar@2	689	data_error_negative = {v: data_error_negative_var}
pilar@2	690	else:
pilar@2	691	data_all[v] = data_var
pilar@2	692	data_error_positive[v] = data_error_positive_var
pilar@2	693	data_error_negative[v] = data_error_negative_var
pilar@2	694
pilar@2	695	# Plot
pilar@2	696	timeiter = dimensions['time']
pilar@2	697	# Default plotting
pilar@2	698	dimensions_t = dimensions
pilar@2	699	for t in range(len(dimensions['time'])):
pilar@2	700	print(f'plotting hr & lr for {timeiter[t]}')
pilar@2	701	dimensions_t['time'] = [timeiter[t]]
pilar@2	702	for vr in range(len(variables_with_data)):
pilar@2	703	if vr == 0:
pilar@2	704	data_t = {variables_with_data[vr]: data_all[variables_with_data[vr]][:, t:t + 1, :]}
pilar@2	705	else:
pilar@2	706	data_t[variables_with_data[vr]] = data_all[variables_with_data[vr]][:, t:t + 1, :]
pilar@2	707
pilar@2	708	fig = plot_vertical_profiles_plotly_allres(
pilar@2	709	variables=variables_with_data,
pilar@2	710	data=data_t,
pilar@2	711	dimensions=dimensions_t,
pilar@2	712	y_array=altitude,
pilar@2	713	positive_error_array=data_error_positive,
pilar@2	714	negative_error_array=data_error_negative,
pilar@2	715	xaxis_label=variables_axis,
pilar@2	716	yaxis_label=altitude['metadata']['long_name'] + ' [' + altitude['metadata']['units'] + ']',
pilar@2	717	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")
pilar@2	718
pilar@2	719	# Show the plot (in Jupyter Notebook or Jupyter Lab)
pilar@2	720	fig.show()
pilar@2	721
pilar@2	722	# Save the plot
pilar@2	723	filename = "ELDAmwl_" + netcdf_file[0:3].upper() + "_" + datetime.fromtimestamp(
pilar@2	724	dimensions['time'][0]).strftime('%Y%m%d-%H%M%S') + "_all.html"
pilar@2	725	fig.write_html(html_path + filename)
pilar@2	726