Mercurial > public > atmospheric_lidar / changeset
changeset
Initial files
Tue, 15 May 2012 16:43:26 +0200
- author
- ulalume3 <binietoglou@imaa.cnr.it>
- date
- Tue, 15 May 2012 16:43:26 +0200
- changeset 0
- 9d2b98ecf23d
- child 1
- 82b144ee09b2
Initial files
| .hgignore | file | annotate | diff | comparison | revisions | |
| __init__.py | file | annotate | diff | comparison | revisions | |
| ciao.py | file | annotate | diff | comparison | revisions | |
| generic.py | file | annotate | diff | comparison | revisions | |
| licel.py | file | annotate | diff | comparison | revisions | |
| musa.py | file | annotate | diff | comparison | revisions | |
| musa_2009_netcdf_parameters.py | file | annotate | diff | comparison | revisions | |
| musa_netcdf_parameters.py | file | annotate | diff | comparison | revisions | |
| pearl.py | file | annotate | diff | comparison | revisions | |
| pearl_netcdf_parameters.py | file | annotate | diff | comparison | revisions |
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.hgignore Tue May 15 16:43:26 2012 +0200 @@ -0,0 +1,3 @@ +syntax: glob +*.pyc +*.py~
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ciao.py Tue May 15 16:43:26 2012 +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/generic.py Tue May 15 16:43:26 2012 +0200 @@ -0,0 +1,543 @@ +# General imports +import datetime +from operator import itemgetter + +# Science imports +import numpy as np +import matplotlib as mpl +from matplotlib import pyplot as plt +import netCDF4 as netcdf + +# CNR-IMAA specific imports +import milos + + +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 (ex. 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 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 save_as_netcdf(self, filename): + """Saves the measurement in the netcdf format as required by the SCC. + Input: filename + """ + 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) + + + 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'] = '\'%s\'' % self.info['start_time'].strftime('%Y%m%d') + input_values['RawData_Start_Time_UT'] = '\'%s\'' % self.info['start_time'].strftime('%H%M%S') + input_values['RawData_Stop_Time_UT'] = '\'%s\'' % 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: + exec('f.%s = %s' % (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]]) + temp_v[:time_length, n] = 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','nb_of_time_scales')) + temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time','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 + + +class Lidar_channel: + + 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): + background = np.mean(self.matrix[:,4000:], 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__() + parameters_values = {'name': self.wavelength, + 'binwidth': self.binwidth, + 'data': subset_data[subset_time[0]],} + + c = Lidar_channel(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 av_profile(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'): + 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 plot(self, signal_type = 'rc', filename = None, zoom = [0,12000,0,-1], show_plot = True, cmap = plt.cm.jet): + #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) + 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]): + + 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]) + + ax1.set_ylabel('Altitude (km)') + ax1.set_xlabel('Time UTC') + #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(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 = 0, + vmin = data[:,10:400].max() * 0.1, + #vmax = 1.4*10**7, + vmax = data[:,10:400].max() * 0.9, + extent = [ts1,ts2,self.z[zoom[0]]/1000.0, self.z[hmax_idx]/1000.0], + ) + + cb1 = plt.colorbar(im1) + cb1.ax.set_ylabel('a.u.') + + 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/licel.py Tue May 15 16:43:26 2012 +0200 @@ -0,0 +1,207 @@ +import numpy as np +import datetime +from generic import BaseLidarMeasurement, Lidar_channel +import musa_netcdf_parameters +import musa_2009_netcdf_parameters + +licel_file_header_format = ['Filename', + 'Site StartDate StartTime EndDate EndTime Altitude Longtitude Latitude ZenithAngle', + '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 """ + + raw_info = {} + channels = {} + channel_info = [] + + f = open(filename, 'r') + + #Read the first 3 lines of the header + raw_info = {} + for c1 in range(3): + raw_info.update(match_lines(f.readline(), licel_file_header_format[c1])) + + 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)) + + # Check the complete header is read + a = f.readline() + + # Import the data + for current_channel_info in channel_info: + raw_data = np.fromfile(f, 'l', 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) + + 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 + f.close() + + self.raw_info = raw_info + self.channels = channels + + + +class LicelFileChannel: + + def __init__(self, raw_info = None, raw_data = None): + self.raw_info = raw_info + self.raw_data = raw_data + + @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 + + def import_file(self, filename): + if filename in self.files: + print "File has been imported already:" + filename + else: + current_file = LicelFile(filename) + for channel_name, channel in current_file.channels.items(): + if channel_name not in self.channels: + self.channels[channel_name] = Licel_channel(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]) + + +class Licel_channel(Lidar_channel): + + 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 microseconds + 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 = 60 + + 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.info['Wavelength'] + + 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/musa.py Tue May 15 16:43:26 2012 +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/musa_2009_netcdf_parameters.py Tue May 15 16:43:26 2012 +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} + +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/musa_netcdf_parameters.py Tue May 15 16:43:26 2012 +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} + +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/pearl.py Tue May 15 16:43:26 2012 +0200 @@ -0,0 +1,139 @@ +import datetime +import os +import glob + +import numpy as np + +from generic import BaseLidarMeasurement, Lidar_channel +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] = Lidar_channel(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/pearl_netcdf_parameters.py Tue May 15 16:43:26 2012 +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} + +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'} +''' + +
