Mercurial > public > atmospheric_lidar / changeset
changeset
Moving files (second attempt)
Mon, 22 Feb 2016 16:52:52 +0200
- author
- Ioannis <devnull@localhost>
- date
- Mon, 22 Feb 2016 16:52:52 +0200
- changeset 36
- a281a26f4626
- parent 35
- b1146c96f727
- child 37
- 7c76fdbdf1a3
Moving files (second attempt)
--- a/.hgignore Mon Feb 22 16:50:03 2016 +0200 +++ b/.hgignore Mon Feb 22 16:52:52 2016 +0200 @@ -3,3 +3,4 @@ *.py~ *.rst~ re:^readme\.md~$ +*.orig
--- a/aias_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -general_parameters = \ -{'System': '\'AIAS\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 36.992, - 'Longitude_degrees_east': 21.649, - 'Altitude_meter_asl': 3} - -channel_parameters = \ -{'00532.p_an': {'channel_ID': 327, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00532.p_ph': {'channel_ID': 328, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.06,}, - '00532.s_an': {'channel_ID': 329, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00532.s_ph': {'channel_ID': 330, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.06,}, - } \ No newline at end of file
--- a/at_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -general_parameters = \ -{'System': '\'LAMP Lidar\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 45.601039, - 'Longitude_degrees_east': 03.723771, - 'Altitude_meter_asl': 420} - -channel_parameters = \ -{'00387.o_ph': {'channel_ID': 316, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00355.p_ph': {'channel_ID': 315, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00355.s_an': {'channel_ID': 312, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0.17,}, - '00355.p_an': {'channel_ID': 314, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00355.s_ph': {'channel_ID': 313, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.17,}, - } - -#For testing. To be read from milos files. -''' -measurement_parameters = \ -{'Pressure_at_Lidar_Station': 930, - 'Temperature_at_Lidar_Station': 15, - 'Measurement_ID': '12345'} -''' - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/aias_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,38 @@ +general_parameters = \ +{'System': '\'AIAS\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 36.992, + 'Longitude_degrees_east': 21.649, + 'Altitude_meter_asl': 3} + +channel_parameters = \ +{'00532.p_an': {'channel_ID': 327, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00532.p_ph': {'channel_ID': 328, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.06,}, + '00532.s_an': {'channel_ID': 329, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00532.s_ph': {'channel_ID': 330, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.06,}, + } \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/at_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,55 @@ +general_parameters = \ +{'System': '\'LAMP Lidar\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 45.601039, + 'Longitude_degrees_east': 03.723771, + 'Altitude_meter_asl': 420} + +channel_parameters = \ +{'00387.o_ph': {'channel_ID': 316, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00355.p_ph': {'channel_ID': 315, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00355.s_an': {'channel_ID': 312, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0.17,}, + '00355.p_an': {'channel_ID': 314, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00355.s_ph': {'channel_ID': 313, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.17,}, + } + +#For testing. To be read from milos files. +''' +measurement_parameters = \ +{'Pressure_at_Lidar_Station': 930, + 'Temperature_at_Lidar_Station': 15, + 'Measurement_ID': '12345'} +''' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/cf_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,55 @@ +general_parameters = \ +{'System': '\'LAMP Lidar\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 45.601039, + 'Longitude_degrees_east': 03.723771, + 'Altitude_meter_asl': 420} + +channel_parameters = \ +{'00387.o_ph': {'channel_ID': 316, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00355.p_ph': {'channel_ID': 315, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00355.s_an': {'channel_ID': 312, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0.17,}, + '00355.p_an': {'channel_ID': 314, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00355.s_ph': {'channel_ID': 313, + 'Background_Low': 15000, + 'Background_High': 20000, + 'Laser_Shots': 600, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.17,}, + } + +#For testing. To be read from milos files. +''' +measurement_parameters = \ +{'Pressure_at_Lidar_Station': 930, + 'Temperature_at_Lidar_Station': 15, + 'Measurement_ID': '12345'} +''' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/cf_raymetrics.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,7 @@ +from licel import LicelLidarMeasurement + +import cf_netcdf_parameters + +class CfLidarMeasurement(LicelLidarMeasurement): + + extra_netcdf_parameters = cf_netcdf_parameters
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/ciao.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,20 @@ +import milos + +class CiaoMixin: + + def get_PT(self): + ''' Gets the pressure and temperature at station level from the Milos station. + The results are stored in the info dictionary. + ''' + + start_time = self.info['start_time'] + stop_time = self.info['stop_time'] + dt = stop_time - start_time + mean_time = start_time + dt/2 + + # this guarantees that more that half the measurement period is taken into account + atm = milos.Atmospheric_condition(mean_time) + temperature = atm.get_mean('Air_Temperature', start_time, stop_time) + pressure = atm.get_mean('Air_Pressure', start_time, stop_time) + self.info['Temperature'] = temperature + self.info['Pressure'] = pressure
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/eole.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,19 @@ +from licel import LicelLidarMeasurement +import eole_netcdf_parameters + + +class EoleLidarMeasurement(LicelLidarMeasurement): + extra_netcdf_parameters = eole_netcdf_parameters + + def get_PT(self): + ''' Sets the pressure and temperature at station level . + The results are stored in the info dictionary. + ''' + + self.info['Temperature'] = 25.0 + self.info['Pressure'] = 1020.0 + + + #def save_netcdf_extra(self, f): + # CHARMEX CLOUD MIN ALTITUDE + # temp_v=f.createVariable('max_altitude_m_asl', 'd', ('time', 'nb_of_time_scales'))
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/eole_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,67 @@ +general_parameters = \ +{'System': '\'EOLE\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 37.96, + 'Longitude_degrees_east': 23.78, + 'Altitude_meter_asl': 212.0} + +channel_parameters = \ +{'01064.o_an': {'channel_ID': 45, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00355.o_an': {'channel_ID': 41, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00355.o_ph': {'channel_ID': 42, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00387.o_ph': {'channel_ID': 46, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00532.o_an': {'channel_ID': 43, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00532.o_ph': {'channel_ID': 44, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00607.o_ph': {'channel_ID': 47, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00407.o_ph': {'channel_ID': 444, + 'Background_Low': 19000, + 'Background_High': 20000, + 'Laser_Shots': 1000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + } +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/generic.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,714 @@ +# General imports +import datetime +from operator import itemgetter + +# Science imports +import numpy as np +import matplotlib as mpl +from matplotlib.ticker import ScalarFormatter +from matplotlib import pyplot as plt +import netCDF4 as netcdf + +netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4' + + +class BaseLidarMeasurement(): + """ This is the general measurement object. + It is meant to become a general measurement object + independent of the input files. + + Each subclass should implement the following: + * the import_file method. + * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters. + + You can override the get_PT method to define a custom procedure to get ground temperature and pressure. + The one implemented by default is by using the MILOS meteorological station data. + + """ + + def __init__(self, filelist = None): + self.info = {} + self.dimensions = {} + self.variables = {} + self.channels = {} + self.attributes = {} + self.files = [] + self.dark_measurement = None + + if filelist: + self.import_files(filelist) + + def import_files(self, filelist): + for f in filelist: + self.import_file(f) + self.update() + + def import_file(self,filename): + raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.') + + def update(self): + ''' + Update the the info, variables and dimensions of the lidar measurement based + on the information found in the channels. + + Reading of the scan_angles parameter is not implemented. + ''' + + # Initialize + start_time =[] + stop_time = [] + points = [] + all_time_scales = [] + channel_name_list = [] + + # Get the information from all the channels + for channel_name, channel in self.channels.items(): + channel.update() + start_time.append(channel.start_time) + stop_time.append(channel.stop_time) + points.append(channel.points) + all_time_scales.append(channel.time) + channel_name_list.append(channel_name) + + # Find the unique time scales used in the channels + time_scales = set(all_time_scales) + + # Update the info dictionary + self.info['start_time'] = min(start_time) + self.info['stop_time'] = max(stop_time) + self.info['duration'] = self.info['stop_time'] - self.info['start_time'] + + # Update the dimensions dictionary + self.dimensions['points'] = max(points) + self.dimensions['channels'] = len(self.channels) + # self.dimensions['scan angles'] = 1 + self.dimensions['nb_of_time_scales'] = len(time_scales) + + # Update the variables dictionary + # Write time scales in seconds + raw_Data_Start_Time = [] + raw_Data_Stop_Time = [] + + for current_time_scale in list(time_scales): + raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time + raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds + raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list + # Check if this time scale has measurements every 30 or 60 seconds. + + duration = self._get_duration(raw_start_in_seconds) + + raw_stop_in_seconds = raw_start_in_seconds + duration + raw_Data_Stop_Time.append(raw_stop_in_seconds) + + self.variables['Raw_Data_Start_Time'] = raw_Data_Start_Time + self.variables['Raw_Data_Stop_Time'] = raw_Data_Stop_Time + + # Make a dictionary to match time scales and channels + channel_timescales = [] + for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales): + # The following lines are PEARL specific. The reason they are here is not clear. + # if channel_name =='1064BLR': + # channel_name = '1064' + for (ts,n) in zip(time_scales, range(len(time_scales))): + if current_time_scale == ts: + channel_timescales.append([channel_name,n]) + self.variables['id_timescale'] = dict(channel_timescales) + + def _get_duration(self, raw_start_in_seconds): + ''' Return the duration for a given time scale. In some files (e.g. Licel) this + can be specified from the files themselves. In others this must be guessed. + + ''' + # The old method, kept here for reference + #dt = np.mean(np.diff(raw_start_in_seconds)) + #for d in duration_list: + # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds + # duration = d + + duration = np.diff(raw_start_in_seconds)[0] + + return duration + + def subset_by_channels(self, channel_subset): + ''' Get a measurement object containing only the channels with names + contained in the channel_sublet list ''' + + m = self.__class__() # Create an object of the same type as this one. + m.channels = dict([(channel, self.channels[channel]) for channel + in channel_subset]) + m.update() + return m + + def subset_by_time(self, start_time, stop_time): + + if start_time > stop_time: + raise ValueError('Stop time should be after start time') + + if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']): + raise ValueError('The time interval specified is not part of the measurement') + + m = self.__class__() # Create an object of the same type as this one. + for (channel_name, channel) in self.channels.items(): + m.channels[channel_name] = channel.subset_by_time(start_time, stop_time) + m.update() + return m + + def subset_by_bins(self, b_min = 0, b_max = None): + """Remove some height bins from the file. This could be needed to + remove aquisition artifacts at the start or the end of the files. + """ + + m = self.__class__() # Create an object of the same type as this one. + + for (channel_name, channel) in self.channels.items(): + m.channels[channel_name] = channel.subset_by_bins(b_min, b_max) + + m.update() + + return m + + def r_plot(self): + #Make a basic plot of the data. + #Should include some dictionary with params to make plot stable. + pass + + def r_pdf(self): + # Create a pdf report using a basic plot and meta-data. + pass + + def save(self): + #Save the current state of the object to continue the analysis later. + pass + + def get_PT(self): + ''' Sets the pressure and temperature at station level . + The results are stored in the info dictionary. + ''' + + self.info['Temperature'] = 10.0 + self.info['Pressure'] = 930.0 + + def subtract_dark(self): + + if not self.dark_measurement: + raise IOError('No dark measurements have been imported yet.') + + for (channel_name, dark_channel) in self.dark_measurement.channels.iteritems(): + dark_profile = dark_channel.average_profile() + channel = self.channels[channel_name] + + for measurement_time, data in channel.data.iteritems(): + channel.data[measurement_time] = data - dark_profile + + channel.update() + + def save_as_netcdf(self, filename = None): + """Saves the measurement in the netcdf format as required by the SCC. + Input: filename. If no filename is provided <measurement_id>.nc will + be used. + """ + params = self.extra_netcdf_parameters + needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure'] + + for parameter in needed_parameters: + stored_value = self.info.get(parameter, None) + if stored_value is None: + raise ValueError('A value needs to be specified for %s' % parameter) + + if not filename: + filename = "%s.nc" % self.info['Measurement_ID'] + + dimensions = {'points': 1, + 'channels': 1, + 'time': None, + 'nb_of_time_scales': 1, + 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented + + global_att = {'Measurement_ID': None, + 'RawData_Start_Date': None, + 'RawData_Start_Time_UT': None, + 'RawData_Stop_Time_UT': None, + 'RawBck_Start_Date': None, + 'RawBck_Start_Time_UT': None, + 'RawBck_Stop_Time_UT': None, + 'Sounding_File_Name': None, + 'LR_File_Name': None, + 'Overlap_File_Name': None, + 'Location': None, + 'System': None, + 'Latitude_degrees_north': None, + 'Longitude_degrees_east': None, + 'Altitude_meter_asl': None} + + channel_variables = \ + {'channel_ID': (('channels', ), 'i'), + 'Background_Low': (('channels', ), 'd'), + 'Background_High': (('channels', ), 'd'), + 'LR_Input': (('channels', ), 'i'), + 'DAQ_Range': (('channels', ), 'd'), + 'Depolarization_Factor': (('channels', ), 'd'), } + + + channels = self.channels.keys() + + input_values = dict(self.dimensions, **self.variables) + + # Add some mandatory global attributes + input_values['Measurement_ID'] = self.info['Measurement_ID'] + input_values['RawData_Start_Date'] = self.info['start_time'].strftime('%Y%m%d') + input_values['RawData_Start_Time_UT'] = self.info['start_time'].strftime('%H%M%S') + input_values['RawData_Stop_Time_UT'] = self.info['stop_time'].strftime('%H%M%S') + + # Add some optional global attributes + input_values['System'] = params.general_parameters['System'] + input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north'] + input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east'] + input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl'] + + # Open a netCDF4 file + f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file + + # Create the dimensions in the file + for (d,v) in dimensions.iteritems(): + v = input_values.pop(d, v) + f.createDimension(d,v) + + # Create global attributes + for (attrib,value) in global_att.iteritems(): + val = input_values.pop(attrib,value) + if val: + setattr(f, attrib, val) + + """ Variables """ + # Write the values of fixes channel parameters + for (var,t) in channel_variables.iteritems(): + temp_v = f.createVariable(var,t[1],t[0]) + for (channel, n) in zip(channels, range(len(channels))): + temp_v[n] = params.channel_parameters[channel][var] + + # Write the id_timescale values + temp_id_timescale = f.createVariable('id_timescale','i',('channels',)) + for (channel, n) in zip(channels, range(len(channels))): + temp_id_timescale[n] = self.variables['id_timescale'][channel] + + # Laser pointing angle + temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',)) + temp_v[:] = params.general_parameters['Laser_Pointing_Angle'] + + # Molecular calculation + temp_v = f.createVariable('Molecular_Calc','i') + temp_v[:] = params.general_parameters['Molecular_Calc'] + + # Laser pointing angles of profiles + temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales')) + for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'], + range(len(self.variables['Raw_Data_Start_Time']))): + temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle + + # Raw data start/stop time + temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales')) + temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales')) + for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'], + self.variables['Raw_Data_Stop_Time'], + range(len(self.variables['Raw_Data_Start_Time']))): + temp_raw_start[:len(start_time),n] = start_time + temp_raw_stop[:len(stop_time),n] = stop_time + + #Laser shots + temp_v = f.createVariable('Laser_Shots','i',('time','channels')) + for (channel,n) in zip(channels, range(len(channels))): + time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]]) + # Array slicing stoped working as usual ex. temp_v[:10] = 100 does not work. ??? np.ones was added. + temp_v[:time_length, n] = np.ones(time_length) * params.channel_parameters[channel]['Laser_Shots'] + + #Raw lidar data + temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points')) + for (channel,n) in zip(channels, range(len(channels))): + c = self.channels[channel] + temp_v[:len(c.time),n, :c.points] = c.matrix + + self.add_dark_measurements_to_netcdf(f, channels) + + #Pressure at lidar station + temp_v = f.createVariable('Pressure_at_Lidar_Station','d') + temp_v[:] = self.info['Pressure'] + + #Temperature at lidar station + temp_v = f.createVariable('Temperature_at_Lidar_Station','d') + temp_v[:] = self.info['Temperature'] + + self.save_netcdf_extra(f) + f.close() + + def add_dark_measurements_to_netcdf(self, f, channels): + + # Get dark measurements. If it is not given in self.dark_measurement + # try to get it using the get_dark_measurements method. If none is found + # return without adding something. + if self.dark_measurement is None: + self.dark_measurement = self.get_dark_measurements() + + if self.dark_measurement is None: + return + + dark_measurement = self.dark_measurement + + # Calculate the length of the time_bck dimensions + number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()] + max_number_of_profiles = np.max(number_of_profiles) + + # Create the dimension + f.createDimension('time_bck', max_number_of_profiles) + + # Save the dark measurement data + temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points')) + for (channel,n) in zip(channels, range(len(channels))): + c = dark_measurement.channels[channel] + temp_v[:len(c.time),n, :c.points] = c.matrix + + # Dark profile start/stop time + temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time_bck','nb_of_time_scales')) + temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time_bck','nb_of_time_scales')) + for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'], + dark_measurement.variables['Raw_Data_Stop_Time'], + range(len(dark_measurement.variables['Raw_Data_Start_Time']))): + temp_raw_start[:len(start_time),n] = start_time + temp_raw_stop[:len(stop_time),n] = stop_time + + # Dark measurement start/stop time + f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d') + f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S') + f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S') + + + + def save_netcdf_extra(self, f): + pass + + def _gettime(self, date_str, time_str): + t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S') + return t + + def plot(self): + for channel in self.channels: + self.channels[channel].plot(show_plot = False) + plt.show() + + def get_dark_measurements(self): + return None + + @property + def mean_time(self): + start_time = self.info['start_time'] + stop_time = self.info['stop_time'] + dt = stop_time - start_time + t_mean = start_time + dt / 2 + return t_mean + + +class LidarChannel: + + def __init__(self, channel_parameters): + c = 299792458 # Speed of light + self.wavelength = channel_parameters['name'] + self.name = str(self.wavelength) + self.binwidth = float(channel_parameters['binwidth']) # in microseconds + self.data = {} + self.resolution = self.binwidth * c / 2 + self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution / 2.0 # Change: add half bin in the z + self.points = len(channel_parameters['data']) + self.rc = [] + self.duration = 60 + + def calculate_rc(self, idx_min = 4000, idx_max = 5000): + background = np.mean(self.matrix[:,idx_min:idx_max], axis = 1) # Calculate the background from 30000m and above + self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2) + + + def update(self): + self.start_time = min(self.data.keys()) + self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration) + self.time = tuple(sorted(self.data.keys())) + sorted_data = sorted(self.data.iteritems(), key=itemgetter(0)) + self.matrix = np.array(map(itemgetter(1),sorted_data)) + + def _nearest_dt(self,dtime): + margin = datetime.timedelta(seconds = 300) + if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time): + print "Requested date not covered in this file" + raise + dt = abs(self.time - np.array(dtime)) + dtmin = min(dt) + + if dtmin > datetime.timedelta(seconds = 60): + print "Nearest profile more than 60 seconds away. dt = %s." % dtmin + ind_t = np.where(dt == dtmin) + ind_a= ind_t[0] + if len(ind_a) > 1: + ind_a = ind_a[0] + chosen_time = self.time[ind_a] + return chosen_time, ind_a + + def subset_by_time(self, start_time, stop_time): + + time_array = np.array(self.time) + condition = (time_array >= start_time) & (time_array <= stop_time) + + subset_time = time_array[condition] + subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time]) + + #Create a list with the values needed by channel's __init__() + parameter_values = {'name': self.wavelength, + 'binwidth': self.binwidth, + 'data': subset_data[subset_time[0]],} + + # We should use __class__ instead of class name, so that this works properly + # with subclasses + # Ex: c = self.__class__(parameters_values) + # This does not currently work with Licel files though + c = LidarChannel(parameter_values) + c.data = subset_data + c.update() + return c + + def subset_by_bins(self, b_min = 0, b_max = None): + """Remove some height bins from the file. This could be needed to + remove aquisition artifacts at the start or the end of the files. + """ + + subset_data = {} + + for ctime, cdata in self.data.items(): + subset_data[ctime] = cdata[b_min:b_max] + + #Create a list with the values needed by channel's __init__() + parameters_values = {'name': self.wavelength, + 'binwidth': self.binwidth, + 'data': subset_data[subset_data.keys()[0]],} # We just need any array with the correct length + + c = LidarChannel(parameters_values) + c.data = subset_data + c.update() + return c + + def profile(self,dtime, signal_type = 'rc'): + t, idx = self._nearest_dt(dtime) + if signal_type == 'rc': + data = self.rc + else: + data = self.matrix + + prof = data[idx,:][0] + return prof, t + + def get_slice(self, starttime, endtime, signal_type = 'rc'): + if signal_type == 'rc': + data = self.rc + else: + data = self.matrix + tim = np.array(self.time) + starttime = self._nearest_dt(starttime)[0] + endtime = self._nearest_dt(endtime)[0] + condition = (tim >= starttime) & (tim <= endtime) + sl = data[condition, :] + t = tim[condition] + return sl,t + + def profile_for_duration(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'): + """ Calculates the profile around a specific time (tim). """ + starttime = tim - duration/2 + endtime = tim + duration/2 + d,t = self.get_slice(starttime, endtime, signal_type = signal_type) + prof = np.mean(d, axis = 0) + tmin = min(t) + tmax = max(t) + tav = tmin + (tmax-tmin)/2 + return prof,(tav, tmin,tmax) + + def average_profile(self): + """ Return the average profile (NOT range corrected) for all the duration of the measurement. """ + prof = np.mean(self.matrix, axis = 0) + return prof + + def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet, z0 = None, title = None, vmin = 0, vmax = 1.3 * 10 ** 7): + #if filename is not None: + # matplotlib.use('Agg') + + fig = plt.figure() + ax1 = fig.add_subplot(111) + self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom, z0 = z0, vmin = vmin, vmax = vmax) + + if title: + ax1.set_title(title) + else: + ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name)) + + if filename is not None: + pass + #plt.savefig(filename) + else: + if show_plot: + plt.show() + #plt.close() ??? + + def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', + zoom = [0,12000,0,-1], z0 = None, + add_colorbar = True, cmap_label = 'a.u.', cb_format = None, + vmin = 0, vmax = 1.3 * 10 ** 7): + + if signal_type == 'rc': + if len(self.rc) == 0: + self.calculate_rc() + data = self.rc + else: + data = self.matrix + + hmax_idx = self.index_at_height(zoom[1]) + + # If z0 is given, then the plot is a.s.l. + if z0: + ax1.set_ylabel('Altitude a.s.l. [km]') + else: + ax1.set_ylabel('Altitude a.g.l. [km]') + z0 = 0 + + ax1.set_xlabel('Time UTC [hh:mm]') + #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 + #dateFormatter = mpl.dates.DateFormatter('%H.%M') + #hourlocator = mpl.dates.HourLocator() + + #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m') + #daylocator = mpl.dates.DayLocator() + hourFormatter = mpl.dates.DateFormatter('%H:%M') + hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) + + + #ax1.axes.xaxis.set_major_formatter(dayFormatter) + #ax1.axes.xaxis.set_major_locator(daylocator) + ax1.axes.xaxis.set_major_formatter(hourFormatter) + ax1.axes.xaxis.set_major_locator(hourlocator) + + + ts1 = mpl.dates.date2num(self.start_time) + ts2 = mpl.dates.date2num(self.stop_time) + + + im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]], + aspect = 'auto', + origin = 'lower', + cmap = cmap, + vmin = vmin, + #vmin = data[:,10:400].max() * 0.1, + vmax = vmax, + #vmax = data[:,10:400].max() * 0.9, + extent = [ts1,ts2,self.z[zoom[0]]/1000.0 + z0/1000., self.z[hmax_idx]/1000.0 + z0/1000.], + ) + + if add_colorbar: + if cb_format: + cb1 = plt.colorbar(im1, format = cb_format) + else: + cb1 = plt.colorbar(im1) + cb1.ax.set_ylabel(cmap_label) + + # Make the ticks of the colorbar smaller, two points smaller than the default font size + cb_font_size = mpl.rcParams['font.size'] - 2 + for ticklabels in cb1.ax.get_yticklabels(): + ticklabels.set_fontsize(cb_font_size) + cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size) + + + def draw_plot_new(self, ax1, cmap = plt.cm.jet, signal_type = 'rc', + zoom = [0, 12000, 0, None], z0 = None, + add_colorbar = True, cmap_label = 'a.u.', + cb_format = None, power_limits = (-2, 2), + date_labels = False, + vmin = 0, vmax = 1.3 * 10 ** 7): + + if signal_type == 'rc': + if len(self.rc) == 0: + self.calculate_rc() + data = self.rc + else: + data = self.matrix + + hmax_idx = self.index_at_height(zoom[1]) + hmin_idx = self.index_at_height(zoom[0]) + + # If z0 is given, then the plot is a.s.l. + if z0: + ax1.set_ylabel('Altitude a.s.l. [km]') + else: + ax1.set_ylabel('Altitude a.g.l. [km]') + z0 = 0 + + ax1.set_xlabel('Time UTC [hh:mm]') + #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 + #dateFormatter = mpl.dates.DateFormatter('%H.%M') + #hourlocator = mpl.dates.HourLocator() + + + if date_labels: + dayFormatter = mpl.dates.DateFormatter('%H:%M\n%d/%m/%Y') + daylocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) + ax1.axes.xaxis.set_major_formatter(dayFormatter) + ax1.axes.xaxis.set_major_locator(daylocator) + else: + hourFormatter = mpl.dates.DateFormatter('%H:%M') + hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) + ax1.axes.xaxis.set_major_formatter(hourFormatter) + ax1.axes.xaxis.set_major_locator(hourlocator) + + + # Get the values of the time axis + dt = datetime.timedelta(seconds = self.duration) + + time_cut = self.time[zoom[2]:zoom[3]] + time_last = time_cut[-1] + dt # The last element needed for pcolormesh + + time_all = time_cut + (time_last,) + + t_axis = mpl.dates.date2num(time_all) + + # Get the values of the z axis + z_cut = self.z[hmin_idx:hmax_idx] - self.resolution / 2. + z_last = z_cut[-1] + self.resolution + + z_axis = np.append(z_cut, z_last) / 1000. + z0 / 1000. # Convert to km + + # Plot + im1 = ax1.pcolormesh(t_axis, z_axis, data.T[hmin_idx:hmax_idx, zoom[2]:zoom[3]], + cmap = cmap, + vmin = vmin, + vmax = vmax, + ) + + if add_colorbar: + if cb_format: + cb1 = plt.colorbar(im1, format = cb_format) + else: + cb_formatter = ScalarFormatter() + cb_formatter.set_powerlimits(power_limits) + cb1 = plt.colorbar(im1, format = cb_formatter) + cb1.ax.set_ylabel(cmap_label) + + # Make the ticks of the colorbar smaller, two points smaller than the default font size + cb_font_size = mpl.rcParams['font.size'] - 2 + for ticklabels in cb1.ax.get_yticklabels(): + ticklabels.set_fontsize(cb_font_size) + cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size) + + def index_at_height(self, height): + idx = np.array(np.abs(self.z - height).argmin()) + if idx.size > 1: + idx =idx[0] + return idx + +def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID): + #Read the lidar files and select channels + temp_m = LidarClass(files) + m = temp_m.subset_by_channels(channels) + m.get_PT() + m.info['Measurement_ID'] = measurement_ID + m.save_as_netcdf(filename) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/licel.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,266 @@ +import numpy as np +import datetime +from generic import BaseLidarMeasurement, LidarChannel +import musa_netcdf_parameters +import musa_2009_netcdf_parameters + +licel_file_header_format = ['Filename', + 'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', # Appart from Site that is read manually + 'LS1 Rate1 LS2 Rate2 DataSets', ] +licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID' + + +class LicelFile: + + def __init__(self, filename): + self.start_time = None + self.stop_time = None + self.import_file(filename) + self.calculate_physical() + self.filename = filename + + def calculate_physical(self): + for channel in self.channels.itervalues(): + channel.calculate_physical() + + def import_file(self, filename): + """Imports a licel file. + Input: filename + Output: object """ + + channels = {} + + with open(filename, 'rb') as f: + + self.read_header(f) + + # Check the complete header is read + a = f.readline() + + # Import the data + for current_channel_info in self.channel_info: + raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints'])) + a = np.fromfile(f, 'b', 1) + b = np.fromfile(f, 'b', 1) + + if (a[0] != 13) | (b[0] != 10): + print "Warning: No end of line found after record. File could be corrupt" + channel = LicelFileChannel(current_channel_info, raw_data, self.duration()) + + channel_name = channel.channel_name + + if channel_name in channels.keys(): + # If the analog/photon naming scheme is not enough, find a new one! + raise IOError('Trying to import two channels with the same name') + + channels[channel_name] = channel + + self.channels = channels + + + def read_header(self, f): + """ Read the header of a open file f. + + Returns raw_info and channel_info. Updates some object properties. """ + + #Read the first 3 lines of the header + raw_info = {} + channel_info = [] + + # Read first line + raw_info['Filename'] = f.readline().strip() + + # Read second line + second_line = f.readline() + + # Many licel files don't follow the licel standard. Specifically, the + # measurement site is not always 8 characters, and can include white + # spaces. For this, the site name is detect everything before the first + # date. For efficiency, the first date is found by the first '/'. + # e.g. assuming a string like 'Site name 01/01/2010 ...' + + site_name = second_line.split('/')[0][:-2] + clean_site_name = site_name.strip() + raw_info['Site'] = clean_site_name + raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1])) + + # Read third line + third_line = f.readline() + raw_info.update(match_lines(third_line, licel_file_header_format[2])) + + # Update the object properties based on the raw info + start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) + stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) + date_format = '%d/%m/%Y %H:%M:%S' + + self.start_time = datetime.datetime.strptime(start_string, date_format) + self.stop_time = datetime.datetime.strptime(stop_string, date_format) + self.latitude = float(raw_info['Latitude']) + self.longitude = float(raw_info['Longtitude']) + + # Read the rest of the header. + for c1 in range(int(raw_info['DataSets'])): + channel_info.append(match_lines(f.readline(), licel_file_channel_format)) + + self.raw_info = raw_info + self.channel_info = channel_info + + def duration(self): + """ Return the duration of the file. """ + dt = self.stop_time - self.start_time + return dt.seconds + + +class LicelFileChannel: + + def __init__(self, raw_info = None, raw_data = None, duration = None): + self.raw_info = raw_info + self.raw_data = raw_data + self.duration = duration + + @property + def wavelength(self): + if self.raw_info is not None: + wave_str = self.raw_info['Wavelength'] + wavelength = wave_str.split('.')[0] + return int(wavelength) + else: + return None + + @property + def channel_name(self): + ''' + Construct the channel name adding analog photon info to avoid duplicates + ''' + acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton']) + channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type) + return channel_name + + def analog_photon_string(self, analog_photon_number): + if analog_photon_number == '0': + string = 'an' + else: + string = 'ph' + return string + + def calculate_physical(self): + data = self.raw_data + + number_of_shots = float(self.raw_info['NShots']) + norm = data / number_of_shots + dz = float(self.raw_info['BinW']) + + if self.raw_info['AnalogPhoton']=='0': + # If the channel is in analog mode + ADCb = int(self.raw_info['ADCbits']) + ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV + channel_data = norm*ADCrange/((2**ADCb)-1) + + # print ADCb, ADCRange,cdata,norm + else: + # If the channel is in photoncounting mode + # Frequency deduced from range resolution! (is this ok?) + # c = 300 # The result will be in MHZ + # SR = c/(2*dz) # To account for pulse folding + # channel_data = norm*SR + # CHANGE: + # For the SCC the data are needed in photons + channel_data = norm *number_of_shots + #print res,c,cdata,norm + + #Calculate Z + number_of_bins = int(self.raw_info['DataPoints']) + self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) + self.dz = dz + self.number_of_bins = number_of_bins + self.data = channel_data + + +class LicelLidarMeasurement(BaseLidarMeasurement): + ''' + + ''' + extra_netcdf_parameters = musa_netcdf_parameters + raw_info = {} # Keep the raw info from the files + durations = {} # Keep the duration of the files + + def import_file(self, filename): + if filename in self.files: + print "File has been imported already:" + filename + else: + current_file = LicelFile(filename) + self.raw_info[current_file.filename] = current_file.raw_info + self.durations[current_file.filename] = current_file.duration() + + for channel_name, channel in current_file.channels.items(): + if channel_name not in self.channels: + self.channels[channel_name] = LicelChannel(channel) + self.channels[channel_name].data[current_file.start_time] = channel.data + self.files.append(current_file.filename) + + def append(self, other): + + self.start_times.extend(other.start_times) + self.stop_times.extend(other.stop_times) + + for channel_name, channel in self.channels.items(): + channel.append(other.channels[channel_name]) + + def _get_duration(self, raw_start_in_seconds): + ''' Return the duration for a given time scale. If only a single + file is imported, then this cannot be guessed from the time difference + and the raw_info of the file are checked. + ''' + + if len(raw_start_in_seconds) == 1: # If only one file imported + duration = self.durations.itervalues().next() # Get the first (and only) raw_info + duration_sec = duration + else: + duration_sec = np.diff(raw_start_in_seconds)[0] + + return duration_sec + + +class LicelChannel(LidarChannel): + + def __init__(self, channel_file): + c = 299792458.0 # Speed of light + self.wavelength = channel_file.wavelength + self.name = channel_file.channel_name + self.binwidth = channel_file.dz * 2 / c # in seconds + self.data = {} + self.resolution = channel_file.dz + self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z + self.points = channel_file.number_of_bins + self.rc = [] + self.duration = channel_file.duration + + def append(self, other): + if self.info != other.info: + raise ValueError('Channel info are different. Data can not be combined.') + + self.data = np.vstack([self.data, other.data]) + + def __unicode__(self): + return "<Licel channel: %s>" % self.name + + def __str__(self): + return unicode(self).encode('utf-8') + + +class Licel2009LidarMeasurement(LicelLidarMeasurement): + extra_netcdf_parameters = musa_2009_netcdf_parameters + + +def match_lines(f1,f2): + list1 = f1.split() + list2 = f2.split() + + if len(list1) != len(list2): + print "Warning: Combining lists of different lengths." + print "List 1: %s" % list1 + print "List 2: %s" % list2 + combined = zip(list2,list1) + combined = dict(combined) + return combined +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/musa.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,5 @@ +from licel import LicelLidarMeasurement +from ciao import CiaoMixin + +class MusaLidarMeasurement(CiaoMixin, LicelLidarMeasurement): + pass
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/musa_2009_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,97 @@ +general_parameters = \ +{'System': '\'MUSA\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 40.601039, + 'Longitude_degrees_east': 15.723771, + 'Altitude_meter_asl': 760.0} # This should be float + +channel_parameters = \ +{'01064.o_an': {'channel_ID': 203, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '00355.o_an': {'channel_ID': 193, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0,}, + '00355.o_ph': {'channel_ID': 194, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00387.o_an': {'channel_ID': 195, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':20.0, + 'Depolarization_Factor': 0,}, + '00387.o_ph': {'channel_ID': 196, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00532.p_an': {'channel_ID': 197, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0,}, + '00532.p_ph': {'channel_ID': 198, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '00532.s_an': {'channel_ID': 199, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0.0441,}, + '00532.s_ph': {'channel_ID': 200, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.0441,}, + '00607.o_an': {'channel_ID': 201, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':20.0, + 'Depolarization_Factor': 0,}, + '00607.o_ph': {'channel_ID': 202, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + } + +#For testing. To be read from milos files. +''' +measurement_parameters = \ +{'Pressure_at_Lidar_Station': 930, + 'Temperature_at_Lidar_Station': 15, + 'Measurement_ID': '12345'} +''' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/musa_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,97 @@ +general_parameters = \ +{'System': '\'MUSA\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 40.601039, + 'Longitude_degrees_east': 15.723771, + 'Altitude_meter_asl': 760.0} # This should be float + +channel_parameters = \ +{'1064.o_an': {'channel_ID': 203, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':500.0, + 'Depolarization_Factor': 0,}, + '355.o_an': {'channel_ID': 193, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0,}, + '355.o_ph': {'channel_ID': 194, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '387.o_an': {'channel_ID': 195, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':20.0, + 'Depolarization_Factor': 0,}, + '387.o_ph': {'channel_ID': 196, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532.p_an': {'channel_ID': 197, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0,}, + '532.p_ph': {'channel_ID': 198, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532.s_an': {'channel_ID': 199, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':100.0, + 'Depolarization_Factor': 0.0441,}, + '532.s_ph': {'channel_ID': 200, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0.0441,}, + '607.o_an': {'channel_ID': 201, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':20.0, + 'Depolarization_Factor': 0,}, + '607.o_ph': {'channel_ID': 202, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1200, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + } + +#For testing. To be read from milos files. +''' +measurement_parameters = \ +{'Pressure_at_Lidar_Station': 930, + 'Temperature_at_Lidar_Station': 15, + 'Measurement_ID': '12345'} +''' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/pearl.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,139 @@ +import datetime +import os +import glob + +import numpy as np + +from generic import BaseLidarMeasurement, LidarChannel +from ciao import CiaoMixin + +import pearl_netcdf_parameters +from report_file import Report_file + + +repository = '/mnt/storage/lidar_data/pearl/' + + +class PearlLidarMeasurement(CiaoMixin, BaseLidarMeasurement): + + extra_netcdf_parameters = pearl_netcdf_parameters + + def import_file(self,filename): + ''' Import a pearl file. ''' + + if filename in self.files: + print "File has been imported already:" + filename + else: + parameters, channels_dict = self.read_pearl_data(filename) + start_time = self._gettime(parameters['Acq_date'],parameters['Acq_start_time']) + + for channel_info in channels_dict.itervalues(): + + if channel_info['name'] == '1064ALR': + name = '1064' + tm = start_time + elif channel_info['name'] == '1064BLR': + name = '1064' + tm = start_time + datetime.timedelta(seconds = 30) + else: + name = channel_info['name'] + tm = start_time + if name not in self.channels: + self.channels[name] = LidarChannel(channel_info) + self.channels[name].data[tm] = channel_info['data'] + self.files.append(filename) + + def read_pearl_data(self, filename): + ''' + Reads a pearl file. + + Returns: + parameters - a dictionary of general parameters + channels - a dictionary with keys the channel number and values lists + [channel name, channel bin width, channel data]. + ''' + f = open(filename,'r') # Open the file + s = f.read(26) # Read the first 26 bytes + + #Get the values in a dictionary + parameters = {} + parameters['Acq_date'] = s[0:10] # First 10 bytes are the acquisition date. + parameters['Acq_start_time'] = s[10:20].strip() # Next 10 bytes are start time. Strip from trailing spaces. + parameters['Channel_no'] = np.fromstring(s[20:22], dtype = np.int16) # Next 2 bytes are the number of channels. Short integer. + parameters['Point_no'] = np.fromstring(s[22:26], dtype = np.int32) # Next 4 bytes are the number of points. Integer. + p = parameters # Just for less typing + + # Read the channel parameters + len = 20*p['Channel_no'] + s = f.read(len) + channels = {} + for (c1,n) in zip(range(0,len, 20),range(p['Channel_no'])): + channels[str(n)] = {'name' : s[c1+10:c1+20].strip(), + 'binwidth' : s[c1:c1+10].strip()} + + #Read the data + data = np.fromfile(f,dtype = np.float32) + #print filename + ': ' + str(data.size) +',' + str(p['Point_no']) +str(p['Channel_no']) + data = data.reshape(p['Point_no'],p['Channel_no']) + for ch in channels.iterkeys(): + channels[ch]['data'] = data[:,int(ch)] + #Close the file + f.close() + return parameters,channels + + +def get_measurement_for_interval(start_time, stop_time): + ''' Searches for a pearl measurement based on a time interval + ''' + + correct_series = None + day = datetime.timedelta(hours = 24) + + if start_time > stop_time: + raise ValueError('Stop time should be after start time') + + + + #The list of directories based on the given time. Same, previous, Next day + possible_paths = [get_path(t) for t in [start_time - day, start_time, start_time + day] + if get_path(t) is not None] + for path in possible_paths: + try: + rf = Report_file(path) + except: + rf = None + + if rf is not None: + for serie in rf.series: + if (start_time > serie.starttime) and (stop_time < serie.endtime): + correct_series = serie + + if correct_series: + files = correct_series.files.get('apd', []) + correct_series.files.get('mcb', []) + m_series = PearlLidarMeasurement(files) + m_subset = m_series.subset_by_time(start_time, stop_time) + return m_subset + else: + return None + + +def get_channel(tim, channel = '1064'): + if channel =='1064': + extension = '*.apd' + else: + extension = '*.mcb' + + dirstr = get_path(tim) + + if not os.path.isdir(dirstr): + raise IOError('No measurement for that date (directory does not exist.).') + #No measurement for that date (directory does not exist.). + files = glob.glob(dirstr + extension) + m = PearlLidarMeasurement(files) + c = m.channels[channel] + return c + + +def get_path(tim): + dirstr = repository + tim.strftime('%Y')+ '/' +tim.strftime('%d%m%Y') + '/' + return dirstr
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/pearl_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,98 @@ +general_parameters = \ +{'System': '\'PEARL\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 40.601039, + 'Longitude_degrees_east': 15.723771, + 'Altitude_meter_asl': 760.0} # This should be float + +channel_parameters = \ +{'355HR': {'channel_ID': 8, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '355LR': {'channel_ID': 14, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '386HR': {'channel_ID': 9, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '386LR': {'channel_ID': 15, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '407HR': {'channel_ID': 17, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532HR': {'channel_ID': 7, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532LR': {'channel_ID': 13, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532SHR': {'channel_ID': 12, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '532PHR': {'channel_ID': 11, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '607HR': {'channel_ID': 10, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':0, + 'Depolarization_Factor': 0,}, + '1064': {'channel_ID': 6, + 'Background_Low': 30000, + 'Background_High': 50000, + 'Laser_Shots': 1500, + 'LR_Input':1, + 'DAQ_Range':100, + 'Depolarization_Factor': 0,}, + + } + +#For testing. To be read from milos files. +''' +measurement_parameters = \ +{'Pressure_at_Lidar_Station': 930, + 'Temperature_at_Lidar_Station': 15, + 'Measurement_ID': '12345'} +''' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/rali.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,33 @@ +import radiometer + +from licel import LicelLidarMeasurement + +import rali_netcdf_parameters + +class RaliLidarMeasurement(LicelLidarMeasurement): + extra_netcdf_parameters = rali_netcdf_parameters + + def get_PT(self): + ''' Gets the pressure and temperature from Radiometer data. + If no data file is found, mean values from past measurements are + used. + ''' + + start_time = self.info['start_time'] + stop_time = self.info['stop_time'] + dt = stop_time - start_time + mean_time = start_time + dt/2 + + meteo_triplet = radiometer.get_mean_PT(mean_time) + + if meteo_triplet: + pressure, temperature, humidity = meteo_triplet + else: + print "Radiometer meteo data not available. Using past values." + pressure = radiometer.P_mean[mean_time.month - 1, mean_time.hour] + temperature = radiometer.T_mean[mean_time.month - 1, mean_time.hour] + + self.info['Temperature'] = temperature - 273.15 + self.info['Pressure'] = pressure + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/rali_netcdf_parameters.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,95 @@ +general_parameters = \ +{'System': '\'RALI\'', + 'Laser_Pointing_Angle': 0, + 'Molecular_Calc': 0, + 'Latitude_degrees_north': 44.348, + 'Longitude_degrees_east': 26.029, + 'Altitude_meter_asl': 93.0} # This should be float + +channel_parameters = \ +{'01064.o_an': {'channel_ID': 89, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 100.0, + 'Depolarization_Factor': 0,}, + '00355.o_an': {'channel_ID': 98, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 100.0, + 'Depolarization_Factor': 0,}, + '00355.o_ph': {'channel_ID': 99, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 0, + 'Depolarization_Factor': 0,}, + '00387.o_an': {'channel_ID': 90, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 100.0, + 'Depolarization_Factor': 0,}, + '00387.o_ph': {'channel_ID': 91, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 0, + 'Depolarization_Factor': 0,}, + '00532.p_an': {'channel_ID': 94, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 100.0, + 'Depolarization_Factor': 0,}, + '00532.p_ph': {'channel_ID': 95, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 0, + 'Depolarization_Factor': 0,}, + '00532.s_an': {'channel_ID': 96, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 20.0, + 'Depolarization_Factor': 0.115,}, + '00532.s_ph': {'channel_ID': 97, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 0, + 'Depolarization_Factor': 0.115,}, + '00607.o_an': {'channel_ID': 92, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 20.0, + 'Depolarization_Factor': 0,}, + '00607.o_ph': {'channel_ID': 93, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input':1, + 'DAQ_Range':20.0, + 'Depolarization_Factor': 0,}, + '00408.o_ph': {'channel_ID': 170, + 'Background_Low': 50000, + 'Background_High': 60000, + 'Laser_Shots': 3000, + 'LR_Input': 1, + 'DAQ_Range': 0, + 'Depolarization_Factor': 0,}, + } +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/requirements.txt Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,1 @@ +.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/atmospheric_lidar/setup.py Mon Feb 22 16:52:52 2016 +0200 @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +from setuptools import setup + +# Read the long description from the readmefile +with open("README.rst", "rb") as f: + long_descr = f.read().decode("utf-8") + +# Run setup +setup(name='cloudmask', + packages=['cloudmask', ], + version='0.1.1', + description='Scripts for applying a cloud mask to uncalibrated lidar signals', + long_description = long_descr, + author='Ioannis Binietoglou', + author_email='ioannis@inoe.ro', + install_requires=[ + "numpy", + "matplotlib", + "sphinx", + "scikit-image", + "pytest" + ], + ) +
--- a/cf_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -general_parameters = \ -{'System': '\'LAMP Lidar\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 45.601039, - 'Longitude_degrees_east': 03.723771, - 'Altitude_meter_asl': 420} - -channel_parameters = \ -{'00387.o_ph': {'channel_ID': 316, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00355.p_ph': {'channel_ID': 315, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00355.s_an': {'channel_ID': 312, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0.17,}, - '00355.p_an': {'channel_ID': 314, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00355.s_ph': {'channel_ID': 313, - 'Background_Low': 15000, - 'Background_High': 20000, - 'Laser_Shots': 600, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.17,}, - } - -#For testing. To be read from milos files. -''' -measurement_parameters = \ -{'Pressure_at_Lidar_Station': 930, - 'Temperature_at_Lidar_Station': 15, - 'Measurement_ID': '12345'} -''' - -
--- a/cf_raymetrics.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -from licel import LicelLidarMeasurement - -import cf_netcdf_parameters - -class CfLidarMeasurement(LicelLidarMeasurement): - - extra_netcdf_parameters = cf_netcdf_parameters
--- a/ciao.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,20 +0,0 @@ -import milos - -class CiaoMixin: - - def get_PT(self): - ''' Gets the pressure and temperature at station level from the Milos station. - The results are stored in the info dictionary. - ''' - - start_time = self.info['start_time'] - stop_time = self.info['stop_time'] - dt = stop_time - start_time - mean_time = start_time + dt/2 - - # this guarantees that more that half the measurement period is taken into account - atm = milos.Atmospheric_condition(mean_time) - temperature = atm.get_mean('Air_Temperature', start_time, stop_time) - pressure = atm.get_mean('Air_Pressure', start_time, stop_time) - self.info['Temperature'] = temperature - self.info['Pressure'] = pressure
--- a/eole.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -from licel import LicelLidarMeasurement -import eole_netcdf_parameters - - -class EoleLidarMeasurement(LicelLidarMeasurement): - extra_netcdf_parameters = eole_netcdf_parameters - - def get_PT(self): - ''' Sets the pressure and temperature at station level . - The results are stored in the info dictionary. - ''' - - self.info['Temperature'] = 25.0 - self.info['Pressure'] = 1020.0 - - - #def save_netcdf_extra(self, f): - # CHARMEX CLOUD MIN ALTITUDE - # temp_v=f.createVariable('max_altitude_m_asl', 'd', ('time', 'nb_of_time_scales'))
--- a/eole_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ -general_parameters = \ -{'System': '\'EOLE\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 37.96, - 'Longitude_degrees_east': 23.78, - 'Altitude_meter_asl': 212.0} - -channel_parameters = \ -{'01064.o_an': {'channel_ID': 45, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00355.o_an': {'channel_ID': 41, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00355.o_ph': {'channel_ID': 42, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00387.o_ph': {'channel_ID': 46, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00532.o_an': {'channel_ID': 43, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00532.o_ph': {'channel_ID': 44, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00607.o_ph': {'channel_ID': 47, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00407.o_ph': {'channel_ID': 444, - 'Background_Low': 19000, - 'Background_High': 20000, - 'Laser_Shots': 1000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - } -
--- a/generic.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,714 +0,0 @@ -# General imports -import datetime -from operator import itemgetter - -# Science imports -import numpy as np -import matplotlib as mpl -from matplotlib.ticker import ScalarFormatter -from matplotlib import pyplot as plt -import netCDF4 as netcdf - -netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4' - - -class BaseLidarMeasurement(): - """ This is the general measurement object. - It is meant to become a general measurement object - independent of the input files. - - Each subclass should implement the following: - * the import_file method. - * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters. - - You can override the get_PT method to define a custom procedure to get ground temperature and pressure. - The one implemented by default is by using the MILOS meteorological station data. - - """ - - def __init__(self, filelist = None): - self.info = {} - self.dimensions = {} - self.variables = {} - self.channels = {} - self.attributes = {} - self.files = [] - self.dark_measurement = None - - if filelist: - self.import_files(filelist) - - def import_files(self, filelist): - for f in filelist: - self.import_file(f) - self.update() - - def import_file(self,filename): - raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.') - - def update(self): - ''' - Update the the info, variables and dimensions of the lidar measurement based - on the information found in the channels. - - Reading of the scan_angles parameter is not implemented. - ''' - - # Initialize - start_time =[] - stop_time = [] - points = [] - all_time_scales = [] - channel_name_list = [] - - # Get the information from all the channels - for channel_name, channel in self.channels.items(): - channel.update() - start_time.append(channel.start_time) - stop_time.append(channel.stop_time) - points.append(channel.points) - all_time_scales.append(channel.time) - channel_name_list.append(channel_name) - - # Find the unique time scales used in the channels - time_scales = set(all_time_scales) - - # Update the info dictionary - self.info['start_time'] = min(start_time) - self.info['stop_time'] = max(stop_time) - self.info['duration'] = self.info['stop_time'] - self.info['start_time'] - - # Update the dimensions dictionary - self.dimensions['points'] = max(points) - self.dimensions['channels'] = len(self.channels) - # self.dimensions['scan angles'] = 1 - self.dimensions['nb_of_time_scales'] = len(time_scales) - - # Update the variables dictionary - # Write time scales in seconds - raw_Data_Start_Time = [] - raw_Data_Stop_Time = [] - - for current_time_scale in list(time_scales): - raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time - raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds - raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list - # Check if this time scale has measurements every 30 or 60 seconds. - - duration = self._get_duration(raw_start_in_seconds) - - raw_stop_in_seconds = raw_start_in_seconds + duration - raw_Data_Stop_Time.append(raw_stop_in_seconds) - - self.variables['Raw_Data_Start_Time'] = raw_Data_Start_Time - self.variables['Raw_Data_Stop_Time'] = raw_Data_Stop_Time - - # Make a dictionary to match time scales and channels - channel_timescales = [] - for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales): - # The following lines are PEARL specific. The reason they are here is not clear. - # if channel_name =='1064BLR': - # channel_name = '1064' - for (ts,n) in zip(time_scales, range(len(time_scales))): - if current_time_scale == ts: - channel_timescales.append([channel_name,n]) - self.variables['id_timescale'] = dict(channel_timescales) - - def _get_duration(self, raw_start_in_seconds): - ''' Return the duration for a given time scale. In some files (e.g. Licel) this - can be specified from the files themselves. In others this must be guessed. - - ''' - # The old method, kept here for reference - #dt = np.mean(np.diff(raw_start_in_seconds)) - #for d in duration_list: - # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds - # duration = d - - duration = np.diff(raw_start_in_seconds)[0] - - return duration - - def subset_by_channels(self, channel_subset): - ''' Get a measurement object containing only the channels with names - contained in the channel_sublet list ''' - - m = self.__class__() # Create an object of the same type as this one. - m.channels = dict([(channel, self.channels[channel]) for channel - in channel_subset]) - m.update() - return m - - def subset_by_time(self, start_time, stop_time): - - if start_time > stop_time: - raise ValueError('Stop time should be after start time') - - if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']): - raise ValueError('The time interval specified is not part of the measurement') - - m = self.__class__() # Create an object of the same type as this one. - for (channel_name, channel) in self.channels.items(): - m.channels[channel_name] = channel.subset_by_time(start_time, stop_time) - m.update() - return m - - def subset_by_bins(self, b_min = 0, b_max = None): - """Remove some height bins from the file. This could be needed to - remove aquisition artifacts at the start or the end of the files. - """ - - m = self.__class__() # Create an object of the same type as this one. - - for (channel_name, channel) in self.channels.items(): - m.channels[channel_name] = channel.subset_by_bins(b_min, b_max) - - m.update() - - return m - - def r_plot(self): - #Make a basic plot of the data. - #Should include some dictionary with params to make plot stable. - pass - - def r_pdf(self): - # Create a pdf report using a basic plot and meta-data. - pass - - def save(self): - #Save the current state of the object to continue the analysis later. - pass - - def get_PT(self): - ''' Sets the pressure and temperature at station level . - The results are stored in the info dictionary. - ''' - - self.info['Temperature'] = 10.0 - self.info['Pressure'] = 930.0 - - def subtract_dark(self): - - if not self.dark_measurement: - raise IOError('No dark measurements have been imported yet.') - - for (channel_name, dark_channel) in self.dark_measurement.channels.iteritems(): - dark_profile = dark_channel.average_profile() - channel = self.channels[channel_name] - - for measurement_time, data in channel.data.iteritems(): - channel.data[measurement_time] = data - dark_profile - - channel.update() - - def save_as_netcdf(self, filename = None): - """Saves the measurement in the netcdf format as required by the SCC. - Input: filename. If no filename is provided <measurement_id>.nc will - be used. - """ - params = self.extra_netcdf_parameters - needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure'] - - for parameter in needed_parameters: - stored_value = self.info.get(parameter, None) - if stored_value is None: - raise ValueError('A value needs to be specified for %s' % parameter) - - if not filename: - filename = "%s.nc" % self.info['Measurement_ID'] - - dimensions = {'points': 1, - 'channels': 1, - 'time': None, - 'nb_of_time_scales': 1, - 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented - - global_att = {'Measurement_ID': None, - 'RawData_Start_Date': None, - 'RawData_Start_Time_UT': None, - 'RawData_Stop_Time_UT': None, - 'RawBck_Start_Date': None, - 'RawBck_Start_Time_UT': None, - 'RawBck_Stop_Time_UT': None, - 'Sounding_File_Name': None, - 'LR_File_Name': None, - 'Overlap_File_Name': None, - 'Location': None, - 'System': None, - 'Latitude_degrees_north': None, - 'Longitude_degrees_east': None, - 'Altitude_meter_asl': None} - - channel_variables = \ - {'channel_ID': (('channels', ), 'i'), - 'Background_Low': (('channels', ), 'd'), - 'Background_High': (('channels', ), 'd'), - 'LR_Input': (('channels', ), 'i'), - 'DAQ_Range': (('channels', ), 'd'), - 'Depolarization_Factor': (('channels', ), 'd'), } - - - channels = self.channels.keys() - - input_values = dict(self.dimensions, **self.variables) - - # Add some mandatory global attributes - input_values['Measurement_ID'] = self.info['Measurement_ID'] - input_values['RawData_Start_Date'] = self.info['start_time'].strftime('%Y%m%d') - input_values['RawData_Start_Time_UT'] = self.info['start_time'].strftime('%H%M%S') - input_values['RawData_Stop_Time_UT'] = self.info['stop_time'].strftime('%H%M%S') - - # Add some optional global attributes - input_values['System'] = params.general_parameters['System'] - input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north'] - input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east'] - input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl'] - - # Open a netCDF4 file - f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file - - # Create the dimensions in the file - for (d,v) in dimensions.iteritems(): - v = input_values.pop(d, v) - f.createDimension(d,v) - - # Create global attributes - for (attrib,value) in global_att.iteritems(): - val = input_values.pop(attrib,value) - if val: - setattr(f, attrib, val) - - """ Variables """ - # Write the values of fixes channel parameters - for (var,t) in channel_variables.iteritems(): - temp_v = f.createVariable(var,t[1],t[0]) - for (channel, n) in zip(channels, range(len(channels))): - temp_v[n] = params.channel_parameters[channel][var] - - # Write the id_timescale values - temp_id_timescale = f.createVariable('id_timescale','i',('channels',)) - for (channel, n) in zip(channels, range(len(channels))): - temp_id_timescale[n] = self.variables['id_timescale'][channel] - - # Laser pointing angle - temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',)) - temp_v[:] = params.general_parameters['Laser_Pointing_Angle'] - - # Molecular calculation - temp_v = f.createVariable('Molecular_Calc','i') - temp_v[:] = params.general_parameters['Molecular_Calc'] - - # Laser pointing angles of profiles - temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales')) - for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'], - range(len(self.variables['Raw_Data_Start_Time']))): - temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle - - # Raw data start/stop time - temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales')) - temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales')) - for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'], - self.variables['Raw_Data_Stop_Time'], - range(len(self.variables['Raw_Data_Start_Time']))): - temp_raw_start[:len(start_time),n] = start_time - temp_raw_stop[:len(stop_time),n] = stop_time - - #Laser shots - temp_v = f.createVariable('Laser_Shots','i',('time','channels')) - for (channel,n) in zip(channels, range(len(channels))): - time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]]) - # Array slicing stoped working as usual ex. temp_v[:10] = 100 does not work. ??? np.ones was added. - temp_v[:time_length, n] = np.ones(time_length) * params.channel_parameters[channel]['Laser_Shots'] - - #Raw lidar data - temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points')) - for (channel,n) in zip(channels, range(len(channels))): - c = self.channels[channel] - temp_v[:len(c.time),n, :c.points] = c.matrix - - self.add_dark_measurements_to_netcdf(f, channels) - - #Pressure at lidar station - temp_v = f.createVariable('Pressure_at_Lidar_Station','d') - temp_v[:] = self.info['Pressure'] - - #Temperature at lidar station - temp_v = f.createVariable('Temperature_at_Lidar_Station','d') - temp_v[:] = self.info['Temperature'] - - self.save_netcdf_extra(f) - f.close() - - def add_dark_measurements_to_netcdf(self, f, channels): - - # Get dark measurements. If it is not given in self.dark_measurement - # try to get it using the get_dark_measurements method. If none is found - # return without adding something. - if self.dark_measurement is None: - self.dark_measurement = self.get_dark_measurements() - - if self.dark_measurement is None: - return - - dark_measurement = self.dark_measurement - - # Calculate the length of the time_bck dimensions - number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()] - max_number_of_profiles = np.max(number_of_profiles) - - # Create the dimension - f.createDimension('time_bck', max_number_of_profiles) - - # Save the dark measurement data - temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points')) - for (channel,n) in zip(channels, range(len(channels))): - c = dark_measurement.channels[channel] - temp_v[:len(c.time),n, :c.points] = c.matrix - - # Dark profile start/stop time - temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time_bck','nb_of_time_scales')) - temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time_bck','nb_of_time_scales')) - for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'], - dark_measurement.variables['Raw_Data_Stop_Time'], - range(len(dark_measurement.variables['Raw_Data_Start_Time']))): - temp_raw_start[:len(start_time),n] = start_time - temp_raw_stop[:len(stop_time),n] = stop_time - - # Dark measurement start/stop time - f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d') - f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S') - f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S') - - - - def save_netcdf_extra(self, f): - pass - - def _gettime(self, date_str, time_str): - t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S') - return t - - def plot(self): - for channel in self.channels: - self.channels[channel].plot(show_plot = False) - plt.show() - - def get_dark_measurements(self): - return None - - @property - def mean_time(self): - start_time = self.info['start_time'] - stop_time = self.info['stop_time'] - dt = stop_time - start_time - t_mean = start_time + dt / 2 - return t_mean - - -class LidarChannel: - - def __init__(self, channel_parameters): - c = 299792458 # Speed of light - self.wavelength = channel_parameters['name'] - self.name = str(self.wavelength) - self.binwidth = float(channel_parameters['binwidth']) # in microseconds - self.data = {} - self.resolution = self.binwidth * c / 2 - self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution / 2.0 # Change: add half bin in the z - self.points = len(channel_parameters['data']) - self.rc = [] - self.duration = 60 - - def calculate_rc(self, idx_min = 4000, idx_max = 5000): - background = np.mean(self.matrix[:,idx_min:idx_max], axis = 1) # Calculate the background from 30000m and above - self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2) - - - def update(self): - self.start_time = min(self.data.keys()) - self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration) - self.time = tuple(sorted(self.data.keys())) - sorted_data = sorted(self.data.iteritems(), key=itemgetter(0)) - self.matrix = np.array(map(itemgetter(1),sorted_data)) - - def _nearest_dt(self,dtime): - margin = datetime.timedelta(seconds = 300) - if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time): - print "Requested date not covered in this file" - raise - dt = abs(self.time - np.array(dtime)) - dtmin = min(dt) - - if dtmin > datetime.timedelta(seconds = 60): - print "Nearest profile more than 60 seconds away. dt = %s." % dtmin - ind_t = np.where(dt == dtmin) - ind_a= ind_t[0] - if len(ind_a) > 1: - ind_a = ind_a[0] - chosen_time = self.time[ind_a] - return chosen_time, ind_a - - def subset_by_time(self, start_time, stop_time): - - time_array = np.array(self.time) - condition = (time_array >= start_time) & (time_array <= stop_time) - - subset_time = time_array[condition] - subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time]) - - #Create a list with the values needed by channel's __init__() - parameter_values = {'name': self.wavelength, - 'binwidth': self.binwidth, - 'data': subset_data[subset_time[0]],} - - # We should use __class__ instead of class name, so that this works properly - # with subclasses - # Ex: c = self.__class__(parameters_values) - # This does not currently work with Licel files though - c = LidarChannel(parameter_values) - c.data = subset_data - c.update() - return c - - def subset_by_bins(self, b_min = 0, b_max = None): - """Remove some height bins from the file. This could be needed to - remove aquisition artifacts at the start or the end of the files. - """ - - subset_data = {} - - for ctime, cdata in self.data.items(): - subset_data[ctime] = cdata[b_min:b_max] - - #Create a list with the values needed by channel's __init__() - parameters_values = {'name': self.wavelength, - 'binwidth': self.binwidth, - 'data': subset_data[subset_data.keys()[0]],} # We just need any array with the correct length - - c = LidarChannel(parameters_values) - c.data = subset_data - c.update() - return c - - def profile(self,dtime, signal_type = 'rc'): - t, idx = self._nearest_dt(dtime) - if signal_type == 'rc': - data = self.rc - else: - data = self.matrix - - prof = data[idx,:][0] - return prof, t - - def get_slice(self, starttime, endtime, signal_type = 'rc'): - if signal_type == 'rc': - data = self.rc - else: - data = self.matrix - tim = np.array(self.time) - starttime = self._nearest_dt(starttime)[0] - endtime = self._nearest_dt(endtime)[0] - condition = (tim >= starttime) & (tim <= endtime) - sl = data[condition, :] - t = tim[condition] - return sl,t - - def profile_for_duration(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'): - """ Calculates the profile around a specific time (tim). """ - starttime = tim - duration/2 - endtime = tim + duration/2 - d,t = self.get_slice(starttime, endtime, signal_type = signal_type) - prof = np.mean(d, axis = 0) - tmin = min(t) - tmax = max(t) - tav = tmin + (tmax-tmin)/2 - return prof,(tav, tmin,tmax) - - def average_profile(self): - """ Return the average profile (NOT range corrected) for all the duration of the measurement. """ - prof = np.mean(self.matrix, axis = 0) - return prof - - def plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet, z0 = None, title = None, vmin = 0, vmax = 1.3 * 10 ** 7): - #if filename is not None: - # matplotlib.use('Agg') - - fig = plt.figure() - ax1 = fig.add_subplot(111) - self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom, z0 = z0, vmin = vmin, vmax = vmax) - - if title: - ax1.set_title(title) - else: - ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name)) - - if filename is not None: - pass - #plt.savefig(filename) - else: - if show_plot: - plt.show() - #plt.close() ??? - - def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', - zoom = [0,12000,0,-1], z0 = None, - add_colorbar = True, cmap_label = 'a.u.', cb_format = None, - vmin = 0, vmax = 1.3 * 10 ** 7): - - if signal_type == 'rc': - if len(self.rc) == 0: - self.calculate_rc() - data = self.rc - else: - data = self.matrix - - hmax_idx = self.index_at_height(zoom[1]) - - # If z0 is given, then the plot is a.s.l. - if z0: - ax1.set_ylabel('Altitude a.s.l. [km]') - else: - ax1.set_ylabel('Altitude a.g.l. [km]') - z0 = 0 - - ax1.set_xlabel('Time UTC [hh:mm]') - #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 - #dateFormatter = mpl.dates.DateFormatter('%H.%M') - #hourlocator = mpl.dates.HourLocator() - - #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m') - #daylocator = mpl.dates.DayLocator() - hourFormatter = mpl.dates.DateFormatter('%H:%M') - hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) - - - #ax1.axes.xaxis.set_major_formatter(dayFormatter) - #ax1.axes.xaxis.set_major_locator(daylocator) - ax1.axes.xaxis.set_major_formatter(hourFormatter) - ax1.axes.xaxis.set_major_locator(hourlocator) - - - ts1 = mpl.dates.date2num(self.start_time) - ts2 = mpl.dates.date2num(self.stop_time) - - - im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]], - aspect = 'auto', - origin = 'lower', - cmap = cmap, - vmin = vmin, - #vmin = data[:,10:400].max() * 0.1, - vmax = vmax, - #vmax = data[:,10:400].max() * 0.9, - extent = [ts1,ts2,self.z[zoom[0]]/1000.0 + z0/1000., self.z[hmax_idx]/1000.0 + z0/1000.], - ) - - if add_colorbar: - if cb_format: - cb1 = plt.colorbar(im1, format = cb_format) - else: - cb1 = plt.colorbar(im1) - cb1.ax.set_ylabel(cmap_label) - - # Make the ticks of the colorbar smaller, two points smaller than the default font size - cb_font_size = mpl.rcParams['font.size'] - 2 - for ticklabels in cb1.ax.get_yticklabels(): - ticklabels.set_fontsize(cb_font_size) - cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size) - - - def draw_plot_new(self, ax1, cmap = plt.cm.jet, signal_type = 'rc', - zoom = [0, 12000, 0, None], z0 = None, - add_colorbar = True, cmap_label = 'a.u.', - cb_format = None, power_limits = (-2, 2), - date_labels = False, - vmin = 0, vmax = 1.3 * 10 ** 7): - - if signal_type == 'rc': - if len(self.rc) == 0: - self.calculate_rc() - data = self.rc - else: - data = self.matrix - - hmax_idx = self.index_at_height(zoom[1]) - hmin_idx = self.index_at_height(zoom[0]) - - # If z0 is given, then the plot is a.s.l. - if z0: - ax1.set_ylabel('Altitude a.s.l. [km]') - else: - ax1.set_ylabel('Altitude a.g.l. [km]') - z0 = 0 - - ax1.set_xlabel('Time UTC [hh:mm]') - #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 - #dateFormatter = mpl.dates.DateFormatter('%H.%M') - #hourlocator = mpl.dates.HourLocator() - - - if date_labels: - dayFormatter = mpl.dates.DateFormatter('%H:%M\n%d/%m/%Y') - daylocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) - ax1.axes.xaxis.set_major_formatter(dayFormatter) - ax1.axes.xaxis.set_major_locator(daylocator) - else: - hourFormatter = mpl.dates.DateFormatter('%H:%M') - hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) - ax1.axes.xaxis.set_major_formatter(hourFormatter) - ax1.axes.xaxis.set_major_locator(hourlocator) - - - # Get the values of the time axis - dt = datetime.timedelta(seconds = self.duration) - - time_cut = self.time[zoom[2]:zoom[3]] - time_last = time_cut[-1] + dt # The last element needed for pcolormesh - - time_all = time_cut + (time_last,) - - t_axis = mpl.dates.date2num(time_all) - - # Get the values of the z axis - z_cut = self.z[hmin_idx:hmax_idx] - self.resolution / 2. - z_last = z_cut[-1] + self.resolution - - z_axis = np.append(z_cut, z_last) / 1000. + z0 / 1000. # Convert to km - - # Plot - im1 = ax1.pcolormesh(t_axis, z_axis, data.T[hmin_idx:hmax_idx, zoom[2]:zoom[3]], - cmap = cmap, - vmin = vmin, - vmax = vmax, - ) - - if add_colorbar: - if cb_format: - cb1 = plt.colorbar(im1, format = cb_format) - else: - cb_formatter = ScalarFormatter() - cb_formatter.set_powerlimits(power_limits) - cb1 = plt.colorbar(im1, format = cb_formatter) - cb1.ax.set_ylabel(cmap_label) - - # Make the ticks of the colorbar smaller, two points smaller than the default font size - cb_font_size = mpl.rcParams['font.size'] - 2 - for ticklabels in cb1.ax.get_yticklabels(): - ticklabels.set_fontsize(cb_font_size) - cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size) - - def index_at_height(self, height): - idx = np.array(np.abs(self.z - height).argmin()) - if idx.size > 1: - idx =idx[0] - return idx - -def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID): - #Read the lidar files and select channels - temp_m = LidarClass(files) - m = temp_m.subset_by_channels(channels) - m.get_PT() - m.info['Measurement_ID'] = measurement_ID - m.save_as_netcdf(filename) -
--- a/licel.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,266 +0,0 @@ -import numpy as np -import datetime -from generic import BaseLidarMeasurement, LidarChannel -import musa_netcdf_parameters -import musa_2009_netcdf_parameters - -licel_file_header_format = ['Filename', - 'StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', # Appart from Site that is read manually - 'LS1 Rate1 LS2 Rate2 DataSets', ] -licel_file_channel_format = 'Active AnalogPhoton LaserUsed DataPoints 1 HV BinW Wavelength d1 d2 d3 d4 ADCbits NShots Discriminator ID' - - -class LicelFile: - - def __init__(self, filename): - self.start_time = None - self.stop_time = None - self.import_file(filename) - self.calculate_physical() - self.filename = filename - - def calculate_physical(self): - for channel in self.channels.itervalues(): - channel.calculate_physical() - - def import_file(self, filename): - """Imports a licel file. - Input: filename - Output: object """ - - channels = {} - - with open(filename, 'rb') as f: - - self.read_header(f) - - # Check the complete header is read - a = f.readline() - - # Import the data - for current_channel_info in self.channel_info: - raw_data = np.fromfile(f, 'i4', int(current_channel_info['DataPoints'])) - a = np.fromfile(f, 'b', 1) - b = np.fromfile(f, 'b', 1) - - if (a[0] != 13) | (b[0] != 10): - print "Warning: No end of line found after record. File could be corrupt" - channel = LicelFileChannel(current_channel_info, raw_data, self.duration()) - - channel_name = channel.channel_name - - if channel_name in channels.keys(): - # If the analog/photon naming scheme is not enough, find a new one! - raise IOError('Trying to import two channels with the same name') - - channels[channel_name] = channel - - self.channels = channels - - - def read_header(self, f): - """ Read the header of a open file f. - - Returns raw_info and channel_info. Updates some object properties. """ - - #Read the first 3 lines of the header - raw_info = {} - channel_info = [] - - # Read first line - raw_info['Filename'] = f.readline().strip() - - # Read second line - second_line = f.readline() - - # Many licel files don't follow the licel standard. Specifically, the - # measurement site is not always 8 characters, and can include white - # spaces. For this, the site name is detect everything before the first - # date. For efficiency, the first date is found by the first '/'. - # e.g. assuming a string like 'Site name 01/01/2010 ...' - - site_name = second_line.split('/')[0][:-2] - clean_site_name = site_name.strip() - raw_info['Site'] = clean_site_name - raw_info.update(match_lines(second_line[len(clean_site_name) + 1:], licel_file_header_format[1])) - - # Read third line - third_line = f.readline() - raw_info.update(match_lines(third_line, licel_file_header_format[2])) - - # Update the object properties based on the raw info - start_string = '%s %s' % (raw_info['StartDate'], raw_info['StartTime']) - stop_string = '%s %s' % (raw_info['EndDate'], raw_info['EndTime']) - date_format = '%d/%m/%Y %H:%M:%S' - - self.start_time = datetime.datetime.strptime(start_string, date_format) - self.stop_time = datetime.datetime.strptime(stop_string, date_format) - self.latitude = float(raw_info['Latitude']) - self.longitude = float(raw_info['Longtitude']) - - # Read the rest of the header. - for c1 in range(int(raw_info['DataSets'])): - channel_info.append(match_lines(f.readline(), licel_file_channel_format)) - - self.raw_info = raw_info - self.channel_info = channel_info - - def duration(self): - """ Return the duration of the file. """ - dt = self.stop_time - self.start_time - return dt.seconds - - -class LicelFileChannel: - - def __init__(self, raw_info = None, raw_data = None, duration = None): - self.raw_info = raw_info - self.raw_data = raw_data - self.duration = duration - - @property - def wavelength(self): - if self.raw_info is not None: - wave_str = self.raw_info['Wavelength'] - wavelength = wave_str.split('.')[0] - return int(wavelength) - else: - return None - - @property - def channel_name(self): - ''' - Construct the channel name adding analog photon info to avoid duplicates - ''' - acquisition_type = self.analog_photon_string(self.raw_info['AnalogPhoton']) - channel_name = "%s_%s" % (self.raw_info['Wavelength'], acquisition_type) - return channel_name - - def analog_photon_string(self, analog_photon_number): - if analog_photon_number == '0': - string = 'an' - else: - string = 'ph' - return string - - def calculate_physical(self): - data = self.raw_data - - number_of_shots = float(self.raw_info['NShots']) - norm = data / number_of_shots - dz = float(self.raw_info['BinW']) - - if self.raw_info['AnalogPhoton']=='0': - # If the channel is in analog mode - ADCb = int(self.raw_info['ADCbits']) - ADCrange = float(self.raw_info['Discriminator'])*1000 # Value in mV - channel_data = norm*ADCrange/((2**ADCb)-1) - - # print ADCb, ADCRange,cdata,norm - else: - # If the channel is in photoncounting mode - # Frequency deduced from range resolution! (is this ok?) - # c = 300 # The result will be in MHZ - # SR = c/(2*dz) # To account for pulse folding - # channel_data = norm*SR - # CHANGE: - # For the SCC the data are needed in photons - channel_data = norm *number_of_shots - #print res,c,cdata,norm - - #Calculate Z - number_of_bins = int(self.raw_info['DataPoints']) - self.z = np.array([dz*bin_number + dz/2.0 for bin_number in range(number_of_bins)]) - self.dz = dz - self.number_of_bins = number_of_bins - self.data = channel_data - - -class LicelLidarMeasurement(BaseLidarMeasurement): - ''' - - ''' - extra_netcdf_parameters = musa_netcdf_parameters - raw_info = {} # Keep the raw info from the files - durations = {} # Keep the duration of the files - - def import_file(self, filename): - if filename in self.files: - print "File has been imported already:" + filename - else: - current_file = LicelFile(filename) - self.raw_info[current_file.filename] = current_file.raw_info - self.durations[current_file.filename] = current_file.duration() - - for channel_name, channel in current_file.channels.items(): - if channel_name not in self.channels: - self.channels[channel_name] = LicelChannel(channel) - self.channels[channel_name].data[current_file.start_time] = channel.data - self.files.append(current_file.filename) - - def append(self, other): - - self.start_times.extend(other.start_times) - self.stop_times.extend(other.stop_times) - - for channel_name, channel in self.channels.items(): - channel.append(other.channels[channel_name]) - - def _get_duration(self, raw_start_in_seconds): - ''' Return the duration for a given time scale. If only a single - file is imported, then this cannot be guessed from the time difference - and the raw_info of the file are checked. - ''' - - if len(raw_start_in_seconds) == 1: # If only one file imported - duration = self.durations.itervalues().next() # Get the first (and only) raw_info - duration_sec = duration - else: - duration_sec = np.diff(raw_start_in_seconds)[0] - - return duration_sec - - -class LicelChannel(LidarChannel): - - def __init__(self, channel_file): - c = 299792458.0 # Speed of light - self.wavelength = channel_file.wavelength - self.name = channel_file.channel_name - self.binwidth = channel_file.dz * 2 / c # in seconds - self.data = {} - self.resolution = channel_file.dz - self.z = np.arange(channel_file.number_of_bins) * self.resolution + self.resolution/2.0 #Change: add half bin in the z - self.points = channel_file.number_of_bins - self.rc = [] - self.duration = channel_file.duration - - def append(self, other): - if self.info != other.info: - raise ValueError('Channel info are different. Data can not be combined.') - - self.data = np.vstack([self.data, other.data]) - - def __unicode__(self): - return "<Licel channel: %s>" % self.name - - def __str__(self): - return unicode(self).encode('utf-8') - - -class Licel2009LidarMeasurement(LicelLidarMeasurement): - extra_netcdf_parameters = musa_2009_netcdf_parameters - - -def match_lines(f1,f2): - list1 = f1.split() - list2 = f2.split() - - if len(list1) != len(list2): - print "Warning: Combining lists of different lengths." - print "List 1: %s" % list1 - print "List 2: %s" % list2 - combined = zip(list2,list1) - combined = dict(combined) - return combined -
--- a/musa.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -from licel import LicelLidarMeasurement -from ciao import CiaoMixin - -class MusaLidarMeasurement(CiaoMixin, LicelLidarMeasurement): - pass
--- a/musa_2009_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -general_parameters = \ -{'System': '\'MUSA\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 40.601039, - 'Longitude_degrees_east': 15.723771, - 'Altitude_meter_asl': 760.0} # This should be float - -channel_parameters = \ -{'01064.o_an': {'channel_ID': 203, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '00355.o_an': {'channel_ID': 193, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0,}, - '00355.o_ph': {'channel_ID': 194, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00387.o_an': {'channel_ID': 195, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':20.0, - 'Depolarization_Factor': 0,}, - '00387.o_ph': {'channel_ID': 196, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00532.p_an': {'channel_ID': 197, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0,}, - '00532.p_ph': {'channel_ID': 198, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '00532.s_an': {'channel_ID': 199, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0.0441,}, - '00532.s_ph': {'channel_ID': 200, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.0441,}, - '00607.o_an': {'channel_ID': 201, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':20.0, - 'Depolarization_Factor': 0,}, - '00607.o_ph': {'channel_ID': 202, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - } - -#For testing. To be read from milos files. -''' -measurement_parameters = \ -{'Pressure_at_Lidar_Station': 930, - 'Temperature_at_Lidar_Station': 15, - 'Measurement_ID': '12345'} -''' - -
--- a/musa_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -general_parameters = \ -{'System': '\'MUSA\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 40.601039, - 'Longitude_degrees_east': 15.723771, - 'Altitude_meter_asl': 760.0} # This should be float - -channel_parameters = \ -{'1064.o_an': {'channel_ID': 203, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':500.0, - 'Depolarization_Factor': 0,}, - '355.o_an': {'channel_ID': 193, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0,}, - '355.o_ph': {'channel_ID': 194, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '387.o_an': {'channel_ID': 195, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':20.0, - 'Depolarization_Factor': 0,}, - '387.o_ph': {'channel_ID': 196, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532.p_an': {'channel_ID': 197, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0,}, - '532.p_ph': {'channel_ID': 198, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532.s_an': {'channel_ID': 199, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':100.0, - 'Depolarization_Factor': 0.0441,}, - '532.s_ph': {'channel_ID': 200, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0.0441,}, - '607.o_an': {'channel_ID': 201, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':20.0, - 'Depolarization_Factor': 0,}, - '607.o_ph': {'channel_ID': 202, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1200, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - } - -#For testing. To be read from milos files. -''' -measurement_parameters = \ -{'Pressure_at_Lidar_Station': 930, - 'Temperature_at_Lidar_Station': 15, - 'Measurement_ID': '12345'} -''' - -
--- a/pearl.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,139 +0,0 @@ -import datetime -import os -import glob - -import numpy as np - -from generic import BaseLidarMeasurement, LidarChannel -from ciao import CiaoMixin - -import pearl_netcdf_parameters -from report_file import Report_file - - -repository = '/mnt/storage/lidar_data/pearl/' - - -class PearlLidarMeasurement(CiaoMixin, BaseLidarMeasurement): - - extra_netcdf_parameters = pearl_netcdf_parameters - - def import_file(self,filename): - ''' Import a pearl file. ''' - - if filename in self.files: - print "File has been imported already:" + filename - else: - parameters, channels_dict = self.read_pearl_data(filename) - start_time = self._gettime(parameters['Acq_date'],parameters['Acq_start_time']) - - for channel_info in channels_dict.itervalues(): - - if channel_info['name'] == '1064ALR': - name = '1064' - tm = start_time - elif channel_info['name'] == '1064BLR': - name = '1064' - tm = start_time + datetime.timedelta(seconds = 30) - else: - name = channel_info['name'] - tm = start_time - if name not in self.channels: - self.channels[name] = LidarChannel(channel_info) - self.channels[name].data[tm] = channel_info['data'] - self.files.append(filename) - - def read_pearl_data(self, filename): - ''' - Reads a pearl file. - - Returns: - parameters - a dictionary of general parameters - channels - a dictionary with keys the channel number and values lists - [channel name, channel bin width, channel data]. - ''' - f = open(filename,'r') # Open the file - s = f.read(26) # Read the first 26 bytes - - #Get the values in a dictionary - parameters = {} - parameters['Acq_date'] = s[0:10] # First 10 bytes are the acquisition date. - parameters['Acq_start_time'] = s[10:20].strip() # Next 10 bytes are start time. Strip from trailing spaces. - parameters['Channel_no'] = np.fromstring(s[20:22], dtype = np.int16) # Next 2 bytes are the number of channels. Short integer. - parameters['Point_no'] = np.fromstring(s[22:26], dtype = np.int32) # Next 4 bytes are the number of points. Integer. - p = parameters # Just for less typing - - # Read the channel parameters - len = 20*p['Channel_no'] - s = f.read(len) - channels = {} - for (c1,n) in zip(range(0,len, 20),range(p['Channel_no'])): - channels[str(n)] = {'name' : s[c1+10:c1+20].strip(), - 'binwidth' : s[c1:c1+10].strip()} - - #Read the data - data = np.fromfile(f,dtype = np.float32) - #print filename + ': ' + str(data.size) +',' + str(p['Point_no']) +str(p['Channel_no']) - data = data.reshape(p['Point_no'],p['Channel_no']) - for ch in channels.iterkeys(): - channels[ch]['data'] = data[:,int(ch)] - #Close the file - f.close() - return parameters,channels - - -def get_measurement_for_interval(start_time, stop_time): - ''' Searches for a pearl measurement based on a time interval - ''' - - correct_series = None - day = datetime.timedelta(hours = 24) - - if start_time > stop_time: - raise ValueError('Stop time should be after start time') - - - - #The list of directories based on the given time. Same, previous, Next day - possible_paths = [get_path(t) for t in [start_time - day, start_time, start_time + day] - if get_path(t) is not None] - for path in possible_paths: - try: - rf = Report_file(path) - except: - rf = None - - if rf is not None: - for serie in rf.series: - if (start_time > serie.starttime) and (stop_time < serie.endtime): - correct_series = serie - - if correct_series: - files = correct_series.files.get('apd', []) + correct_series.files.get('mcb', []) - m_series = PearlLidarMeasurement(files) - m_subset = m_series.subset_by_time(start_time, stop_time) - return m_subset - else: - return None - - -def get_channel(tim, channel = '1064'): - if channel =='1064': - extension = '*.apd' - else: - extension = '*.mcb' - - dirstr = get_path(tim) - - if not os.path.isdir(dirstr): - raise IOError('No measurement for that date (directory does not exist.).') - #No measurement for that date (directory does not exist.). - files = glob.glob(dirstr + extension) - m = PearlLidarMeasurement(files) - c = m.channels[channel] - return c - - -def get_path(tim): - dirstr = repository + tim.strftime('%Y')+ '/' +tim.strftime('%d%m%Y') + '/' - return dirstr
--- a/pearl_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,98 +0,0 @@ -general_parameters = \ -{'System': '\'PEARL\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 40.601039, - 'Longitude_degrees_east': 15.723771, - 'Altitude_meter_asl': 760.0} # This should be float - -channel_parameters = \ -{'355HR': {'channel_ID': 8, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '355LR': {'channel_ID': 14, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '386HR': {'channel_ID': 9, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '386LR': {'channel_ID': 15, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '407HR': {'channel_ID': 17, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532HR': {'channel_ID': 7, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532LR': {'channel_ID': 13, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532SHR': {'channel_ID': 12, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '532PHR': {'channel_ID': 11, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '607HR': {'channel_ID': 10, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':0, - 'Depolarization_Factor': 0,}, - '1064': {'channel_ID': 6, - 'Background_Low': 30000, - 'Background_High': 50000, - 'Laser_Shots': 1500, - 'LR_Input':1, - 'DAQ_Range':100, - 'Depolarization_Factor': 0,}, - - } - -#For testing. To be read from milos files. -''' -measurement_parameters = \ -{'Pressure_at_Lidar_Station': 930, - 'Temperature_at_Lidar_Station': 15, - 'Measurement_ID': '12345'} -''' - -
--- a/rali.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -import radiometer - -from licel import LicelLidarMeasurement - -import rali_netcdf_parameters - -class RaliLidarMeasurement(LicelLidarMeasurement): - extra_netcdf_parameters = rali_netcdf_parameters - - def get_PT(self): - ''' Gets the pressure and temperature from Radiometer data. - If no data file is found, mean values from past measurements are - used. - ''' - - start_time = self.info['start_time'] - stop_time = self.info['stop_time'] - dt = stop_time - start_time - mean_time = start_time + dt/2 - - meteo_triplet = radiometer.get_mean_PT(mean_time) - - if meteo_triplet: - pressure, temperature, humidity = meteo_triplet - else: - print "Radiometer meteo data not available. Using past values." - pressure = radiometer.P_mean[mean_time.month - 1, mean_time.hour] - temperature = radiometer.T_mean[mean_time.month - 1, mean_time.hour] - - self.info['Temperature'] = temperature - 273.15 - self.info['Pressure'] = pressure - -
--- a/rali_netcdf_parameters.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,95 +0,0 @@ -general_parameters = \ -{'System': '\'RALI\'', - 'Laser_Pointing_Angle': 0, - 'Molecular_Calc': 0, - 'Latitude_degrees_north': 44.348, - 'Longitude_degrees_east': 26.029, - 'Altitude_meter_asl': 93.0} # This should be float - -channel_parameters = \ -{'01064.o_an': {'channel_ID': 89, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 100.0, - 'Depolarization_Factor': 0,}, - '00355.o_an': {'channel_ID': 98, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 100.0, - 'Depolarization_Factor': 0,}, - '00355.o_ph': {'channel_ID': 99, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 0, - 'Depolarization_Factor': 0,}, - '00387.o_an': {'channel_ID': 90, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 100.0, - 'Depolarization_Factor': 0,}, - '00387.o_ph': {'channel_ID': 91, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 0, - 'Depolarization_Factor': 0,}, - '00532.p_an': {'channel_ID': 94, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 100.0, - 'Depolarization_Factor': 0,}, - '00532.p_ph': {'channel_ID': 95, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 0, - 'Depolarization_Factor': 0,}, - '00532.s_an': {'channel_ID': 96, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 20.0, - 'Depolarization_Factor': 0.115,}, - '00532.s_ph': {'channel_ID': 97, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 0, - 'Depolarization_Factor': 0.115,}, - '00607.o_an': {'channel_ID': 92, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 20.0, - 'Depolarization_Factor': 0,}, - '00607.o_ph': {'channel_ID': 93, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input':1, - 'DAQ_Range':20.0, - 'Depolarization_Factor': 0,}, - '00408.o_ph': {'channel_ID': 170, - 'Background_Low': 50000, - 'Background_High': 60000, - 'Laser_Shots': 3000, - 'LR_Input': 1, - 'DAQ_Range': 0, - 'Depolarization_Factor': 0,}, - } -
--- a/requirements.txt Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -.
--- a/setup.py Mon Feb 22 16:50:03 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -#!/usr/bin/env python - -from setuptools import setup - -# Read the long description from the readmefile -with open("README.rst", "rb") as f: - long_descr = f.read().decode("utf-8") - -# Run setup -setup(name='cloudmask', - packages=['cloudmask', ], - version='0.1.1', - description='Scripts for applying a cloud mask to uncalibrated lidar signals', - long_description = long_descr, - author='Ioannis Binietoglou', - author_email='ioannis@inoe.ro', - install_requires=[ - "numpy", - "matplotlib", - "sphinx", - "scikit-image", - "pytest" - ], - ) -
