pearl.py

Tue, 12 Feb 2013 16:56:59 +0100

author
ulalume3 <binietoglou@imaa.cnr.it>
date
Tue, 12 Feb 2013 16:56:59 +0100
changeset 13
da3b70fcebf6
parent 0
9d2b98ecf23d
child 27
74f7617f5356
permissions
-rw-r--r--

Changed altitude parameter from int to float.

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
binietoglou@0 7 from generic import BaseLidarMeasurement, Lidar_channel
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:
binietoglou@0 42 self.channels[name] = Lidar_channel(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

mercurial