Sun, 23 Nov 2014 23:24:32 +0200
New plots for alex.
binietoglou@0 | 1 | import datetime |
binietoglou@0 | 2 | import os |
binietoglou@0 | 3 | import glob |
binietoglou@0 | 4 | |
binietoglou@0 | 5 | import numpy as np |
binietoglou@0 | 6 | |
ulalume3@27 | 7 | from generic import BaseLidarMeasurement, LidarChannel |
binietoglou@0 | 8 | from ciao import CiaoMixin |
binietoglou@0 | 9 | |
binietoglou@0 | 10 | import pearl_netcdf_parameters |
binietoglou@0 | 11 | from report_file import Report_file |
binietoglou@0 | 12 | |
binietoglou@0 | 13 | |
binietoglou@0 | 14 | repository = '/mnt/storage/lidar_data/pearl/' |
binietoglou@0 | 15 | |
binietoglou@0 | 16 | |
binietoglou@0 | 17 | class PearlLidarMeasurement(CiaoMixin, BaseLidarMeasurement): |
binietoglou@0 | 18 | |
binietoglou@0 | 19 | extra_netcdf_parameters = pearl_netcdf_parameters |
binietoglou@0 | 20 | |
binietoglou@0 | 21 | def import_file(self,filename): |
binietoglou@0 | 22 | ''' Import a pearl file. ''' |
binietoglou@0 | 23 | |
binietoglou@0 | 24 | if filename in self.files: |
binietoglou@0 | 25 | print "File has been imported already:" + filename |
binietoglou@0 | 26 | else: |
binietoglou@0 | 27 | parameters, channels_dict = self.read_pearl_data(filename) |
binietoglou@0 | 28 | start_time = self._gettime(parameters['Acq_date'],parameters['Acq_start_time']) |
binietoglou@0 | 29 | |
binietoglou@0 | 30 | for channel_info in channels_dict.itervalues(): |
binietoglou@0 | 31 | |
binietoglou@0 | 32 | if channel_info['name'] == '1064ALR': |
binietoglou@0 | 33 | name = '1064' |
binietoglou@0 | 34 | tm = start_time |
binietoglou@0 | 35 | elif channel_info['name'] == '1064BLR': |
binietoglou@0 | 36 | name = '1064' |
binietoglou@0 | 37 | tm = start_time + datetime.timedelta(seconds = 30) |
binietoglou@0 | 38 | else: |
binietoglou@0 | 39 | name = channel_info['name'] |
binietoglou@0 | 40 | tm = start_time |
binietoglou@0 | 41 | if name not in self.channels: |
ulalume3@27 | 42 | self.channels[name] = LidarChannel(channel_info) |
binietoglou@0 | 43 | self.channels[name].data[tm] = channel_info['data'] |
binietoglou@0 | 44 | self.files.append(filename) |
binietoglou@0 | 45 | |
binietoglou@0 | 46 | def read_pearl_data(self, filename): |
binietoglou@0 | 47 | ''' |
binietoglou@0 | 48 | Reads a pearl file. |
binietoglou@0 | 49 | |
binietoglou@0 | 50 | Returns: |
binietoglou@0 | 51 | parameters - a dictionary of general parameters |
binietoglou@0 | 52 | channels - a dictionary with keys the channel number and values lists |
binietoglou@0 | 53 | [channel name, channel bin width, channel data]. |
binietoglou@0 | 54 | ''' |
binietoglou@0 | 55 | f = open(filename,'r') # Open the file |
binietoglou@0 | 56 | s = f.read(26) # Read the first 26 bytes |
binietoglou@0 | 57 | |
binietoglou@0 | 58 | #Get the values in a dictionary |
binietoglou@0 | 59 | parameters = {} |
binietoglou@0 | 60 | parameters['Acq_date'] = s[0:10] # First 10 bytes are the acquisition date. |
binietoglou@0 | 61 | parameters['Acq_start_time'] = s[10:20].strip() # Next 10 bytes are start time. Strip from trailing spaces. |
binietoglou@0 | 62 | parameters['Channel_no'] = np.fromstring(s[20:22], dtype = np.int16) # Next 2 bytes are the number of channels. Short integer. |
binietoglou@0 | 63 | parameters['Point_no'] = np.fromstring(s[22:26], dtype = np.int32) # Next 4 bytes are the number of points. Integer. |
binietoglou@0 | 64 | p = parameters # Just for less typing |
binietoglou@0 | 65 | |
binietoglou@0 | 66 | # Read the channel parameters |
binietoglou@0 | 67 | len = 20*p['Channel_no'] |
binietoglou@0 | 68 | s = f.read(len) |
binietoglou@0 | 69 | channels = {} |
binietoglou@0 | 70 | for (c1,n) in zip(range(0,len, 20),range(p['Channel_no'])): |
binietoglou@0 | 71 | channels[str(n)] = {'name' : s[c1+10:c1+20].strip(), |
binietoglou@0 | 72 | 'binwidth' : s[c1:c1+10].strip()} |
binietoglou@0 | 73 | |
binietoglou@0 | 74 | #Read the data |
binietoglou@0 | 75 | data = np.fromfile(f,dtype = np.float32) |
binietoglou@0 | 76 | #print filename + ': ' + str(data.size) +',' + str(p['Point_no']) +str(p['Channel_no']) |
binietoglou@0 | 77 | data = data.reshape(p['Point_no'],p['Channel_no']) |
binietoglou@0 | 78 | for ch in channels.iterkeys(): |
binietoglou@0 | 79 | channels[ch]['data'] = data[:,int(ch)] |
binietoglou@0 | 80 | #Close the file |
binietoglou@0 | 81 | f.close() |
binietoglou@0 | 82 | return parameters,channels |
binietoglou@0 | 83 | |
binietoglou@0 | 84 | |
binietoglou@0 | 85 | def get_measurement_for_interval(start_time, stop_time): |
binietoglou@0 | 86 | ''' Searches for a pearl measurement based on a time interval |
binietoglou@0 | 87 | ''' |
binietoglou@0 | 88 | |
binietoglou@0 | 89 | correct_series = None |
binietoglou@0 | 90 | day = datetime.timedelta(hours = 24) |
binietoglou@0 | 91 | |
binietoglou@0 | 92 | if start_time > stop_time: |
binietoglou@0 | 93 | raise ValueError('Stop time should be after start time') |
binietoglou@0 | 94 | |
binietoglou@0 | 95 | |
binietoglou@0 | 96 | |
binietoglou@0 | 97 | #The list of directories based on the given time. Same, previous, Next day |
binietoglou@0 | 98 | possible_paths = [get_path(t) for t in [start_time - day, start_time, start_time + day] |
binietoglou@0 | 99 | if get_path(t) is not None] |
binietoglou@0 | 100 | for path in possible_paths: |
binietoglou@0 | 101 | try: |
binietoglou@0 | 102 | rf = Report_file(path) |
binietoglou@0 | 103 | except: |
binietoglou@0 | 104 | rf = None |
binietoglou@0 | 105 | |
binietoglou@0 | 106 | if rf is not None: |
binietoglou@0 | 107 | for serie in rf.series: |
binietoglou@0 | 108 | if (start_time > serie.starttime) and (stop_time < serie.endtime): |
binietoglou@0 | 109 | correct_series = serie |
binietoglou@0 | 110 | |
binietoglou@0 | 111 | if correct_series: |
binietoglou@0 | 112 | files = correct_series.files.get('apd', []) + correct_series.files.get('mcb', []) |
binietoglou@0 | 113 | m_series = PearlLidarMeasurement(files) |
binietoglou@0 | 114 | m_subset = m_series.subset_by_time(start_time, stop_time) |
binietoglou@0 | 115 | return m_subset |
binietoglou@0 | 116 | else: |
binietoglou@0 | 117 | return None |
binietoglou@0 | 118 | |
binietoglou@0 | 119 | |
binietoglou@0 | 120 | def get_channel(tim, channel = '1064'): |
binietoglou@0 | 121 | if channel =='1064': |
binietoglou@0 | 122 | extension = '*.apd' |
binietoglou@0 | 123 | else: |
binietoglou@0 | 124 | extension = '*.mcb' |
binietoglou@0 | 125 | |
binietoglou@0 | 126 | dirstr = get_path(tim) |
binietoglou@0 | 127 | |
binietoglou@0 | 128 | if not os.path.isdir(dirstr): |
binietoglou@0 | 129 | raise IOError('No measurement for that date (directory does not exist.).') |
binietoglou@0 | 130 | #No measurement for that date (directory does not exist.). |
binietoglou@0 | 131 | files = glob.glob(dirstr + extension) |
binietoglou@0 | 132 | m = PearlLidarMeasurement(files) |
binietoglou@0 | 133 | c = m.channels[channel] |
binietoglou@0 | 134 | return c |
binietoglou@0 | 135 | |
binietoglou@0 | 136 | |
binietoglou@0 | 137 | def get_path(tim): |
binietoglou@0 | 138 | dirstr = repository + tim.strftime('%Y')+ '/' +tim.strftime('%d%m%Y') + '/' |
binietoglou@0 | 139 | return dirstr |