Commit 9987cda5 authored by Markus Gaug's avatar Markus Gaug
Browse files

new plot for volcano data public

parent d3c85327
Pipeline #7265 failed with stages
in 0 seconds
from matplotlib import pyplot as plt
import numpy as np
import json
import sys
from os.path import basename
from scipy.optimize import curve_fit
from astropy import units as u
import pandas as pd
import glob
import re
# LIDAR analysis modules
from licel_reader import LICEL
from lidar_analysis import *
from rayleigh import Rayleigh
from observatory import Observatory
from gdas_reader import GDASReader
from atmospheric_profiles import AtmosphericProfiles
import models
np.set_printoptions(threshold=sys.maxsize)
# MAGIC file
magicfile = "MAGIClidar_2021_09_22_21_09_03_fixed.data"
magic_tab = pd.read_table(magicfile,sep=",\s+", skip_blank_lines=True)
#print(magic_tab['x'])
#exit(0)
# required configuration files
gdas_config = 'config/entry_data_gdas.json'
obs_config = 'config/ObservatoryData.json'
# initialize the molecular profile from GDAS
# location of the corresponding GDAS file, needs to be adapted!!!
gdas_file = [ '../../Measurements/LaPalma/GDAS/gdas1.sep21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name) # get the observatory coordinates and altitude
gdas = GDASReader(gdas_file, obs) # read the GDAS file
ap = AtmosphericProfiles(gdas_config, gdas) # initialize the atmospheric profiles of the GDAS file
H0 = ap.H0_WS # altitude of the observatory
ap.calculate_results_for_all_profiles(plot=False, save_to_file=False, metrics=False) # calculuate the atmospheric profiles of the GDAS file
dates = list(ap.entry_data.keys())[0] # available dates in GDAS file
model = models.Models() # intialize the fitting models
# initialize the axis ranges
x_min = 0 # in km
x_max = 7.8 # in km
y_min = 11.8
y_max = 21.5
# molecular fit ranges
x_mol0 = 1.8 # in km
x_mol1 = 2.3 # in km
max_ranges_visualized = [ 8. , 0.5, # 355 nm in km
4. , 8., # 387 nm in km
8., 0.5, # 532 nm in km
0.75, 0.5 ] # near range in km
min_ranges_visualized = [ 0.45 , 0.45, # 355 nm in km
0.45 , 0.45, # 387 nm in km
0.45, 0.45, # 532 nm in km
0. , 0. ] # near range in km
min_magic = 0.7
min_ranges_beta = [ 0.5 ]
max_ranges_beta = [ 4. ]
# standard channel configuration
labels = [ '355 nm', '355 nm (photon counting)',
'387 nm (analog)', '387 nm',
'532 nm', '532 nm (photon counting)',
'near range (analog)', 'near range' ]
uv_id = 0
raman_id = 2
green_id = 4
nr_id = 6
offsets = [ 2 , -7 , 10, -5 , -7, -9, -3, -5 ] # internal Licel electronic time offsets
bg_range = 5000 # Background subtraction (both atmospheric and electronic)
used_elastic_channels = [0,3,4] # Used elastic channels for molecular profile fits and Klett inversion
dx_range = 2. # current uncertainty of range in m
# apply a constant re-binning (-1 means no re-binning)
rebin = -1
# apply a dynamic re-binning threshold in photo-electrons (-1 means no re-binning)
rebin_thres = -1 # 0.1 # 4 # in p.e.
rebin_seq = 50 # increase re-binning every x slices
rebin_fac = 2 # increase re-binning by x after every rebin_seq slices
y_raman = np.zeros(1)
y_green = np.zeros(1)
y_uv = np.zeros(1)
x_raman = np.zeros(1)
x_green = np.zeros(1)
x_uv = np.zeros(1)
y_nr = np.zeros(1)
b_ids = np.zeros(1)
alphaaer_filt = np.zeros(1)
pc_thres = 1
pc_thres_nr = 0.01
pc_sum_raman = 0.4 # photon counting sum in each (new) bin of the sample, after re-binning
pc_sum_green = 0.7 # photon counting sum in each (new) bin of the sample, after re-binning
pc_sum_uv = 0.04 # photon counting sum in each (new) bin of the sample, after re-binning
nr_fac = 1800
absolute_offset = 0 #100 # 80
nr_offset = 10 #50
y_filt_window = 61 # 47 #67
y_filt_poly = 4 # 5 #5
vertical_axis = True
dateRegex = re.compile(r'(\d\d\d\d\d\d).(\d\d\d-\d\d\d\d)') # 20210922.200202
# function to evaluate the molecular profile at altitude h
def eval_molecular_profile(h: np.array, param: float, hlo: float, costheta: float, beta_mol_0: float):
#print (h,hlo,h,beta_mol_0,H0)
return param - np.log(Rayleigh.ptMAGIC) + ap.eval_f(dates, hlo, h*u.km, costheta, beta_mol_0, H0)
def show_file(filename):
#plt.figure() # initialize a new figure for each data file
fig, ax = plt.subplots(constrained_layout=True,figsize=[5,7.5])
licel = LICEL(filename) # initialize Licel reader for filename
licel.print_header() # print header information
ddate = licel.start_datetime.date().isoformat()
dtime = licel.start_datetime.time().isoformat()
costheta = np.cos(licel.zenith_ang * np.pi / 180.) # needed for LIDAR inversion
integrate_mol_from = 0 * u.km # start point of extinction integration
x_max_used = x_max * costheta
is_glue = True
id_amp = 2
id_pc = 3
testplot = False # plot a test plot to control the gluing
rebin_glue = -1 # do not rebin for gluing
deltaR = 1000 # initial gluing window
deltae = 1.5
fmin = 0.01
fmax = 180.
id0 = 30
licel.channels[id_amp].offset_range = offsets[id_amp] # correct for the different zero-offset in channel_id 1
licel.channels[id_pc].offset_range = offsets[id_pc] # correct for the different zero-offset in channel_id 1
if (is_glue):
(a_Raman, b_Raman) = glue_lange(licel,id_amp,id_pc,rebin_glue,
bg_range=bg_range+(int)(licel.channels[id_amp].num_data/2),bg_range2=(int)(licel.channels[id_amp].num_data/2),
deltaR=deltaR,
deltaa=deltae,fmin=fmin,fmax=fmax,plot=testplot,id0=id0) #2000.,7500.,deadtime,testplot)
if (testplot):
plt.show()
for (channel_id,channel) in enumerate(licel.channels): # loop over all channels
channel.offset_range = offsets[channel_id] # correct for the different zero-offset
if (channel_id > 3): # work-around, until David has set correct bin widths for the new Licels
channel.bin_width = 3.75
channel.print_header() # print channel header
if ((channel_id == uv_id+1) or (channel_id == green_id+1)): # or channel_id == 6):
continue
if (channel_id < nr_id):
rebin = 0
else:
rebin = 3
# get the data to plot
if (channel_id == raman_id): # channel 2 has a different ground line offset at about half of the signal range, and need separate background estimation
(x_range,y_sig,y_err,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=rebin_thres,rebin_seq=rebin_seq,rebin_fac=rebin_fac,
bg_range=bg_range+(int)(channel.num_data/2),bg_range2=(int)(channel.num_data/2),
verbose=False)
else:
(x_range,y_sig,y_err,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=rebin_thres,rebin_seq=rebin_seq,rebin_fac=rebin_fac,
bg_range=bg_range,
verbose=True)
if (not channel.analog): # pc channels have one more bin
x_range = x_range[:-1]
y_sig = y_sig[:-1]
y_err = y_err[:-1]
if (channel_id == 0):
x_range = x_range - 2 # move x by 2 meters to the left for fine-adjustment
if (channel_id < nr_id):
x_range = x_range + absolute_offset
else:
x_range = x_range + absolute_offset + nr_offset
if (is_glue and channel.analog): # apply the gluing here!!
if (channel_id == nr_id or channel_id == green_id or channel_id == uv_id):
#print ('y_sig=', y_sig)
#print ('x_range=', x_range)
print ('apply glueing here: ', a_Raman)
y_sig = y_sig*a_Raman
else:
print ('apply glueing here: ', a_Raman, b_Raman)
y_sig = y_sig*a_Raman# + b_Raman
if (channel_id == raman_id): # Raman line, analog, save it to a separate variable
y_raman = y_sig
pc_idx_raman = (np.where((y_sig <= pc_thres) & (x_range > 1000.)))[0][0]
print ('pc_idx=', pc_idx_raman)
if (channel_id == (raman_id + 1)):
y_raman[pc_idx_raman:] = y_sig[pc_idx_raman:] # here we join the analog line in
(x_raman,y_raman,yerr_raman,b_ids,b_widths) = dynamic_rebinning(channel, x_range,y_raman,offset,offseterr,pc_thres=pc_sum_raman,rebin_start=10)
#print ('Y Raman: ',y_raman, ' x Raman: ', x_raman)
if (channel_id == green_id):
(x_green,y_green,yerr_green,b_ids,b_widths) = dynamic_rebinning(channel, x_range,y_sig,offset,offseterr,pc_thres=pc_sum_green,rebin_start=10, verbose=False)
#print ('Y green: ',y_green, ' x green: ', x_green)
if (channel_id == uv_id):
(x_uv,y_uv,yerr_uv,b_ids,b_widths) = dynamic_rebinning(channel, x_range,y_sig,offset,offseterr,pc_thres=pc_sum_uv,rebin_start=10)
#print ('Y UV: ',y_uv, ' x UV: ', x_uv)
if (channel_id == nr_id): # Raman line, analog, save it to a separate variable
y_nr = y_sig
#print (x_range, y_sig)
pc_idx_nr = (np.where((y_sig <= pc_thres_nr) & (x_range > 10.)))[0][0]
print ('pc_idx_nr=', pc_idx_nr)
y_sig = y_sig * nr_fac
#print ('y_sig=', y_sig)
if (channel_id == (nr_id + 1)):
y_nr[pc_idx_nr:] = y_sig[pc_idx_nr:] # here we join the analog line in
y_nr = y_nr * nr_fac
# discard any negative data for the logarithm in plot
# discard range with too high noise
if (channel_id == uv_id):
ids = np.where((y_uv>0) & (x_uv>10.) & (x_uv < max_ranges_visualized[channel_id]*1000) & (x_uv > min_ranges_visualized[channel_id]*1000))
elif (channel_id == raman_id +1):
ids = np.where((y_raman>0) & (x_raman>10.) & (x_raman < max_ranges_visualized[channel_id]*1000) & (x_raman > min_ranges_visualized[channel_id]*1000))
elif (channel_id == green_id):
ids = np.where((y_green>0) & (x_green>10.) & (x_green < max_ranges_visualized[channel_id]*1000) & (x_green > min_ranges_visualized[channel_id]*1000))
elif (channel_id < nr_id):
ids = np.where((y_sig>0) & (x_range>200.) & (x_range < max_ranges_visualized[channel_id]*1000) & (x_range > min_ranges_visualized[channel_id]*1000))
else:
ids = np.where((y_sig>0) & (x_range>10.) & (x_range < max_ranges_visualized[channel_id]*1000) & (x_range > min_ranges_visualized[channel_id]*1000))
x_range = x_range[ids]
x_height= x_range/1000. * costheta # convert x-axis from range to height for MAGIC analysis
y_sig = y_sig[ids]
y_err = y_err[ids]
if (channel_id == uv_id):
#print ('Y: ',y_raman[ids], ' x: ', x_raman[ids])
sr = np.log(y_uv[ids] * x_uv[ids]**2) # standard logarithm of range-corrected signal
x_height = x_uv[ids]/1000. * costheta
elif (channel_id == raman_id + 1):
#print ('Y: ',y_raman[ids], ' x: ', x_raman[ids])
sr = np.log(y_raman[ids] * x_raman[ids]**2) # standard logarithm of range-corrected signal
x_height = x_raman[ids]/1000. * costheta
elif (channel_id == green_id):
#print ('Y: ',y_raman[ids], ' x: ', x_raman[ids])
sr = np.log(y_green[ids] * x_green[ids]**2) # standard logarithm of range-corrected signal
x_height = x_green[ids]/1000. * costheta
elif (channel_id == nr_id):
print ('Y: ',len(y_sig), ' x: ', len(x_range))
sr = np.log(y_sig * x_range**2) # standard logarithm of range-corrected signal
elif (channel_id == nr_id + 1):
print ('Y: ',len(y_nr[ids]), ' x: ', len(x_range))
sr = np.log(y_nr[ids] * x_range**2) # standard logarithm of range-corrected signal
else:
sr = np.log(y_sig * x_range**2) # standard logarithm of range-corrected signal
#y_err = 1/sr * (y_err * x_range**2 + y_sig * 2 * x_range * dx_range )
if ((channel_id != raman_id) and (channel_id != nr_id)):
if (vertical_axis):
ax.plot(sr, x_height,label=labels[channel_id], alpha = 0.5)
else:
ax.plot(x_height, sr,label=labels[channel_id], alpha = 0.5)
if channel_id in used_elastic_channels: # molecular fit and Klett only for selected elastic channels
ray = Rayleigh(channel.wavelength_nm * u.nm) # initialize Rayleigh backscatter coefficient for used wavelength
beta_mol_0 = ray.dbeta_domega((180 * u.deg).to(u.rad)) # mol. backscatter coefficient at altitude 0 from LIDAR
ids2 = (x_height > x_mol0*costheta) & (x_height < x_mol1*costheta)
# fit molecular profile to data and obtain the LIDAR constant c0
ans = curve_fit(lambda xl, c0: eval_molecular_profile(xl, c0, integrate_mol_from, costheta, beta_mol_0), x_height[ids2] + H0.to(u.km).value, sr[ids2],p0=[ 15. ])
y_fit = eval_molecular_profile(x_height[ids2] + H0.to(u.km).value, *ans[0], integrate_mol_from, costheta, beta_mol_0)
#chi2_red = np.sum(((sr[ids2] - y_fit)/y_err[ids2])**2) / (y_fit.size - 1)
#print(*ans[0], chi2_red, y_err[ids2])
#plt.plot(x_height[ids2], y_fit, label = r"molecular fit {:.0f} nm, $\chi^{{2}}$/NDF={:.1f}".format(channel.wavelength_nm,chi2_red))ç
ax.plot(x_height[ids2], y_fit, label = r"molecular fit {:.0f} nm".format(channel.wavelength_nm))
if (channel_id == raman_id+1):
print ('b_widths: len: ', len(b_widths[ids]), ' cont: ', b_widths[ids], ' x: len: ',len(x_height), ' cont: ', x_height)
ax2 = ax.twiny()
ax2.set_xlim(-1e-4,1e-2)
ax2.set_xlabel(r"Extinction coefficient (m$^{-1}$)")
for i in np.arange(0,1):
y_poly = y_filt_poly
y_window = y_filt_window# + 4*i
alphaaer_filt = invert_raman_savgol_filt(x_height, y_raman[ids],
eval_molecular_profile(x_height + H0.to(u.km).value, *ans[0], integrate_mol_from, costheta, beta_mol_0),
b_widths[ids], y_window, y_poly )
#print ('alphaaer', x_height, alphaaer_filt)
#ax2.set_xticks(new_tick_locations)
#ax2.set_xticklabels(tick_function(new_tick_locations))
ids_a = np.where((x_height > min_ranges_beta) & (x_height < max_ranges_beta))
#ax2.plot(alphaaer_filt[ids_a],x_height[ids_a], label=r"$\beta$, w={:d}".format(y_window))
ax2.plot(alphaaer_filt[ids_a],x_height[ids_a], label=r"$\alpha_{aer}$")
if (vertical_axis):
ax.text(y_min+(y_max-y_min)*0.65,x_min+(x_max_used-x_min)*0.12,
u'zenith={:.0f}\u00B0'.format(licel.zenith_ang),
fontsize='x-large')
ax.text(y_min+(y_max-y_min)*0.65,x_min+(x_max_used-x_min)*0.05,
'{:.0f} shots'.format(licel.laser1['shots']),
fontsize='x-large')
ax.set_ylim(x_min,x_max_used)
ax.set_xlim(y_min,y_max)
else:
ax.text(x_min+(x_max_used-x_min)*0.45, y_min+(y_max-y_min)*0.92,
u'zenith={:.0f}\u00B0'.format(licel.zenith_ang),
fontsize='x-large')
ax.text(x_min+(x_max_used-x_min)*0.45, y_min+(y_max-y_min)*0.85,
'{:.0f} shots'.format(licel.laser1['shots']),
fontsize='x-large')
ax.set_xlim(x_min,x_max_used)
ax.set_ylim(y_min,y_max)
#plt.yscale('log')
height_magic = np.array(magic_tab.x/1000.)
sr_magic = np.array(magic_tab.y)
ids_magic = np.where(height_magic > min_magic)
if (vertical_axis):
ax.plot(sr_magic[ids_magic]+np.log(7.5/48.),height_magic[ids_magic],label=u"MAGIC, 21:09 UTC, 0.4\u00B0")
else:
ax.plot(height_magic[ids_magic],sr_magic[ids_magic]+np.log(7.5/48.),label=u"MAGIC, 21:09 UTC, 0.4\u00B0")
ax.legend(loc='best')
ax2.legend(loc='lower right')
ax.grid(axis='x')
ax.set_title('Barcelona Raman LIDAR - {:s}, {:s} UTC'.format(ddate,dtime))
if (vertical_axis):
ax.set_ylabel('Altitude above LIDAR (km)')
else:
ax.set_xlabel('Altitude above LIDAR (km)')
if (vertical_axis):
ax.set_xlabel(r'Logarithm range-corrected signal (log(p.e. $\cdot m) )$')
else:
ax.set_ylabel(r'Logarithm range-corrected signal (log(p.e. $\cdot m) )$')
plt.show()
if __name__ == "__main__":
for name in sorted(glob.glob(sys.argv[1])):
print(name)
show_file(name)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment