Commit 4c021168 authored by Markus Gaug's avatar Markus Gaug
Browse files

new file data_analyzer.py

parent 34605980
Pipeline #6877 failed with stages
in 0 seconds
import numpy as np
from licel_reader import LICEL
from lidar_analysis import *
from matplotlib import pyplot as plt
import sys
ffile = sys.argv[1]
#ffile = sys.argv[1]
is_raw = False
is_groundline = False
bg_range = 5000 # Background subtraction (both atmospheric and electronic)
if (is_groundline):
bg_range = 0
sig_range = 10000
offsets = [ 6 , -2 , 14, 2, 0, 0, -10, -10 ]
rebin = 0
np.set_printoptions(precision=3,threshold=50000)
rebin_thres = 0
labels = [ '355 nm (analog)', '355 nm (photon counting)', '387 nm (analog)', '387 nm (photon counting)' ,
'532 nm (analog)', '532 nm (photon counting)' , 'near range (analog)', 'near range (photon counting)']
is_glue = False
molfile = 'my_ecmwf_interim_data_201807.txt'
lambda_elastic = 355 # wavelength of elastic channel
lambda_raman = 387 # wavelength of raman channel
angstrom_k = 1.14 # assumed angstrom exponent for aerosol extinction
altitude_asl = 135 # altitude of the LIDAR a.s.l., in meters
elastic_id = 1 # ID of the elastic file in previous declaration
raman_id = 0 # ID of the Raman file in previous declaration
offset_bin = 21 # general offset bin for any analysis
pc_thres = 0.1 # threshold from where on the PC signal is used for the Raman signal
pc_sum = 20 # photon counting sum in each (new) bin of the sample, after re-binning
z0 = 6100 # z0 in beta inversion (referene point with assumed mol. content only)
deadtime = 0.0062 # dead time of the photon counting channel, needed for gluing
deltae = 1.5 # relative increase of parabola to determine final window, used for gluing
#day = licel.start_datetime.date().day
#hour = licel.start_datetime.time().hour//6 * 6 # ECMWF hours are given in time intervals of 6 hours only
testplot = False # plot a test plot to control the gluing
plt.figure(figsize=[15,7])
if (is_groundline):
plt.subplot(421)
#for (file_id,ffile) in enumerate(files):
nr = 0
licel = LICEL(ffile)
licel.print_header()
id_amp = 6
id_pc = 7
rebin_glue = 0 # do not rebin for gluing
deltaR = 1000 # initial gluing window
fmin = 4.
fmax = 180.
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 (not print_only):
if (is_glue):
(a0, b0) = glue_lange(licel,id_amp,id_pc,rebin_glue,bg_range,deltaR,deltae,fmin,fmax,testplot) #2000.,7500.,deadtime,testplot)
if (testplot):
plt.show()
id_amp = 2
id_pc = 3
deltaR = 500 # initial gluing window
fmin = 1.
fmax = 180.
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):
(a1, b1) = glue_lange(licel,id_amp,id_pc,rebin_glue,bg_range,deltaR,deltae,fmin,fmax,testplot) #2000.,7500.,deadtime,testplot)
if (testplot):
plt.show()
rebin_thres = 0.
for (channel_id,channel) in enumerate(licel.channels):
channel.offset_range = offsets[channel_id] # correct for the different zero-offset in channel_id 1
channel.print_header() # print channel header
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=rebin_thres,rebin_seq=50,rebin_fac=0,bg_range=bg_range) # plot data
if (is_glue and channel.analog ): # apply the gluing here!!
print ('apply glueing here: ', a0, b0)
y = y*a0 + b0
#print ('\n\n\n',x,y)
# This is the threshold from where on the PC signal is used for the Raman signal
ids = np.where(y > pc_thres)
#pc_idx = (np.where(y <= pc_thres))[0][0]
#print ('Found pc_idx: ' ,pc_idx, ' at: ', x[pc_idx])
if (is_glue and channel_id == 2 ): # apply the gluing here!!
print ('apply glueing here: ', a1, b1)
y = y*a1 + b1
#print ('\n\n\n',x,y)
# This is the threshold from where on the PC signal is used for the Raman signal
ids = np.where(y <= pc_thres)
#pc_idx = (np.where(y <= pc_thres))[0][0]
#print ('Found pc_idx: ' ,pc_idx, ' at: ', x[pc_idx])
#print (y[0:5000])
#if (channel_id == 3):# or channel_id == 3):
#if (channel_id > 1):
# ids = np.where(y > 0.2)
#ids = np.where(y > 0.002)
if (is_groundline):
print ('nr = ',nr)
ax = plt.subplot(4,2,nr+1)
ax.plot(x,y,label=labels[nr])
ax.set_xlim(0.,sig_range)
ax.set_xlabel('distance (m)')
if (channel.analog):
ax.set_ylabel('raw signal (mV)')
if (bg_range == 0):
ax.set_ylim(-0.3,10.)
else:
if (channel_id == 6):
ax.set_ylim(-0.1,0.1)
else:
ax.set_ylim(-0.1,0.1)
else:
ax.set_ylabel('raw signal (MHz)')
if (bg_range == 0):
ax.set_ylim(-1.,2.5)
else:
ax.set_ylim(-1.,2.5)
ax.axhline(y=0., color='k')
ax.legend(loc='best')
else:
if (is_raw):
plt.plot(x,y,label=labels[nr])
else:
plt.plot(x,np.log(y*x**2),label=labels[nr])
#if (channel_id == 1):# or channel_id == 3):
# plt.plot(x,y,label=labels[channel_id])
print ('channel_id = ',channel_id)
#if (channel_id == 3):
# print (y)
nr = nr + 1
if (not is_groundline):
if (is_raw):
plt.xlim(0.,sig_range)
plt.yscale('log')
else:
plt.xlim(0.,sig_range)
plt.ylim(0.,19.)
plt.xlabel('distance (m)')
plt.ylabel(r'$\log(S*r^{2})$ ')
plt.grid()
plt.legend(loc='best')
#plt.yscale('log')
ttit = '_logS'.format(sig_range)
if (is_raw):
ttit = '_raw{:d}'.format(sig_range)
plt.xlabel('distance (m)')
plt.ylabel('raw signal (mV or MHz)')
if (is_groundline):
ttit = '_groundline_BG{:d}'.format(bg_range)
plt.savefig('Signal_{:s}_{:s}.pdf'.format(ffile,ttit), bbox_inches='tight')
plt.show()
markusgaug@cie-49-198.uab.es.364
\ No newline at end of file
......@@ -52,11 +52,11 @@ class AtmosphericDataReader(ABC):
"rh", # Relative humidity
]
name : str
filenames: str
pressure_levels: int
atmospheric_data: dict
#name : str
#filenames: str
#pressure_levels: int
#atmospheric_data: dict
@abstractmethod
def _read_data(self) -> None:
pass
\ No newline at end of file
pass
......@@ -582,12 +582,12 @@ class AtmosphericProfiles():
unit = x_data.unit
except:
unit = u.dimensionless_unscaled
ax.set_xlabel(plt_xlabel + rf"[{unit:FITS}]")
ax.set_xlabel(plt_xlabel + '[{unit:FITS}]')
try:
unit = y_data.unit
except:
unit = u.dimensionless_unscaled
ax.set_ylabel(plt_ylabel + rf"[{unit:FITS}]")
ax.set_ylabel(plt_ylabel + '[{unit:FITS}]')
if isinstance(x_data, Quantity):
x_data = x_data.to_value()
......@@ -894,4 +894,4 @@ class AtmosphericProfiles():
print(data.to_markdown())
print("Model info")
print(self.entry_data[date]["models"])
np.set_printoptions(**oldoptions)
\ No newline at end of file
np.set_printoptions(**oldoptions)
import glob
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import json
import sys
from os.path import basename
from licel_reader import LICEL
from observatory import Observatory
from gdas_reader import GDASReader
from atmospheric_profiles import AtmosphericProfiles
gdas_config = 'config/entry_data_gdas.json'
obs_config = 'config/ObservatoryData.json'
x_min = 0.
x_max = 5000.
y_min = 1e-4
y_max = 1e3
is_groundline = False
sig_range = 140000
labels = [ '355 nm (analog)', '355 nm (photon counting)',
'389 nm (analog)', '389 nm (photon counting)',
'532 nm (analog)', '532 nm (photon counting)',
'near range (analog)', 'near range (photon counting)' ]
rebin = -1
rebin_thres = 0
offsets = [ 2 , -7 , 10, -5 , -7, -9, -3, -5 ]
bg_range = 5000 # 5000 # Background subtraction (both atmospheric and electronic)
if (is_groundline):
bg_range = 5000
def show_file(filename):
licel = LICEL(filename)
licel.print_header()
nr = 0
plt.figure(figsize=[15,7])
if (is_groundline):
plt.subplot(421)
for (channel_id,channel) in enumerate(licel.channels):
#if (channel_id > 0):
# continue
channel.offset_range = offsets[nr] # + global_offset # correct for the different zero-offset
if (channel_id > 3):
channel.bin_width = 3.75
channel.print_header() # print channel header
if (channel_id == 2):
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=0,rebin_seq=50,rebin_fac=0,bg_range=bg_range+(int)(channel.num_data/2),bg_range2=(int)(channel.num_data/2),verbose=False)
else:
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=0,rebin_seq=50,rebin_fac=0,bg_range=bg_range,verbose=False)
#plt.plot(x,np.log(y*x**2),label=labels[channel_id])
if (channel_id == 0):
x = x - 2 # move x by 2 meters to the left
if (channel.analog == False):
y = y * 0.008
if (is_groundline):
ax = plt.subplot(4,2,channel_id+1)
ax.plot(x,y,label=labels[channel_id])
ax.set_xlim(0.,sig_range)
if (channel_id > 3):
ax.set_xlim(0.,sig_range/2.)
ax.set_xlabel('distance (m)')
if (channel.analog):
ax.set_ylabel('raw signal (mV)')
if (bg_range == 0):
ax.set_ylim(4.,6.)
else:
if (channel_id == 6):
ax.set_ylim(-0.01,0.01)
else:
ax.set_ylim(-0.03,0.03)
else:
ax.set_ylabel('raw signal (MHz)')
if (bg_range == 0):
ax.set_ylim(-0.1,0.5)
else:
ax.set_ylim(-1.,2.5)
ax.axhline(y=0., color='k')
ax.legend(loc='best')
else:
plt.plot(x,y,label=labels[channel_id])
plt.xlim(x_min,x_max)
plt.ylim(y_min,y_max)
plt.yscale('log')
nr = nr + 1
plt.legend(loc='best')
plt.grid(axis='x')
plt.title(basename(filename))
plt.show()
gdas_file = [ '/Users/markusgaug/CTA/LIDAR/Measurements/LaPalma/GDAS/gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name)
gdas = GDASReader(gdas_file, obs)
ap = AtmosphericProfiles(gdas_config, gdas)
ap.calculate_results_for_all_profiles(plot=True, save_to_file=True, metrics=True)
#for name in sorted(glob.glob(sys.argv[1])):
# print(name)
# show_file(name)
......@@ -77,9 +77,12 @@ class GDASReader(AtmosphericDataReader):
for i in range(0, len(df), self.pressure_levels + 1):
#Parse 4 first columns numbers to datetime format
#date = df.iloc[i, :4].astype('uint32').to_string(index=False).replace('\n', '/').replace(' ', ' 0').replace(' ', '') #datetime.datetime.strptime(
#print(date)
date = datetime.datetime.strptime(
df.iloc[i, :4].astype('uint32').to_string(index=False)
.replace('\n', '/').replace(' ', '0'),
df.iloc[i, :4].astype('uint32').to_string(index=False)
.replace('\n', '/').replace(' ', ' 0').replace(' ', ''),
'%y/%m/%d/%H')
#Read header
......@@ -99,4 +102,4 @@ class GDASReader(AtmosphericDataReader):
block_data['geopotential_height'].to_numpy()
) / 1000
self.atmospheric_data[date] = (header_data, block_data)
\ No newline at end of file
self.atmospheric_data[date] = (header_data, block_data)
import sys
import json
import bson
filename = sys.argv[1]
if not filename.endswith(".json"):
print ("Argument name must end with .json, cannot convert")
sys.exit()
outfile = filename.replace(".json",".bson")
def read_json(filename):
with open(filename, 'r') as f:
j = json.load(f)
return j
def dump_json(j, filename):
with open(filename, 'w') as f:
bson.dump(j, f, default=str)
def read_bson(filename="out.bson"):
with open(filename, 'rb') as f:
j = bson.loads(f.read())
return j
def writetofile(b, filename="out.bson"):
with open(filename, "wb") as bf:
bf.write(b)
starting_json = read_json(filename)
bson_obj = bson.dumps(starting_json)
writetofile(bson_obj, filename=outfile)
#final_json = read_bson(outfile)
#print(final_json)
......@@ -221,7 +221,7 @@ ret1 = testInefficientPcChannels(licel)
ret2 = testOscillations(licel)
offsets = [ 6 , -3 , 14, 2, 0, 0, 0, 0 ]
offsets = [ 6 , -3 , 14, 2, 0, 0, 0, 0]
licel.channels[0].offset_range = offsets[0] # correct for the different zero-offset in channel_id 1
licel.channels[1].offset_range = offsets[1] # correct for the different zero-offset in channel_id 1
......
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