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

eldamwl_plot.py@fce4cae19357 (annotated)

eldamwl_plot.py

Mon, 29 Jan 2024 11:07:31 +0000

author

pilar <pilar.gumaclaramunt@cnr.it>

date

Mon, 29 Jan 2024 11:07:31 +0000

changeset 0

fce4cae19357

child 2

763a7ac22031

permissions

-rwxr-xr-x

First version

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