Mon, 12 May 2014 12:41:22 +0300
Auto filename when saving netcdf
binietoglou@0 | 1 | # General imports |
binietoglou@0 | 2 | import datetime |
binietoglou@0 | 3 | from operator import itemgetter |
binietoglou@0 | 4 | |
binietoglou@0 | 5 | # Science imports |
binietoglou@0 | 6 | import numpy as np |
binietoglou@0 | 7 | import matplotlib as mpl |
binietoglou@0 | 8 | from matplotlib import pyplot as plt |
binietoglou@0 | 9 | import netCDF4 as netcdf |
binietoglou@0 | 10 | |
binietoglou@1 | 11 | netcdf_format = 'NETCDF3_CLASSIC' # choose one of 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC' and 'NETCDF4' |
binietoglou@0 | 12 | |
binietoglou@0 | 13 | |
binietoglou@0 | 14 | class BaseLidarMeasurement(): |
binietoglou@0 | 15 | """ This is the general measurement object. |
binietoglou@0 | 16 | It is meant to become a general measurement object |
binietoglou@0 | 17 | independent of the input files. |
binietoglou@0 | 18 | |
binietoglou@0 | 19 | Each subclass should implement the following: |
binietoglou@0 | 20 | * the import_file method. |
binietoglou@0 | 21 | * set the "extra_netcdf_parameters" variable to a dictionary that includes the appropriate parameters. |
binietoglou@0 | 22 | |
binietoglou@0 | 23 | You can override the get_PT method to define a custom procedure to get ground temperature and pressure. |
binietoglou@0 | 24 | The one implemented by default is by using the MILOS meteorological station data. |
binietoglou@0 | 25 | |
binietoglou@0 | 26 | """ |
binietoglou@0 | 27 | |
binietoglou@0 | 28 | def __init__(self, filelist= None): |
binietoglou@0 | 29 | self.info = {} |
binietoglou@0 | 30 | self.dimensions = {} |
binietoglou@0 | 31 | self.variables = {} |
binietoglou@0 | 32 | self.channels = {} |
binietoglou@0 | 33 | self.attributes = {} |
binietoglou@0 | 34 | self.files = [] |
binietoglou@0 | 35 | self.dark_measurement = None |
binietoglou@0 | 36 | if filelist: |
binietoglou@0 | 37 | self.import_files(filelist) |
binietoglou@0 | 38 | |
binietoglou@0 | 39 | def import_files(self,filelist): |
binietoglou@0 | 40 | for f in filelist: |
binietoglou@0 | 41 | self.import_file(f) |
binietoglou@0 | 42 | self.update() |
binietoglou@0 | 43 | |
binietoglou@0 | 44 | def import_file(self,filename): |
binietoglou@0 | 45 | raise NotImplementedError('Importing files should be defined in the instrument-specific subclass.') |
binietoglou@0 | 46 | |
binietoglou@0 | 47 | def update(self): |
binietoglou@0 | 48 | ''' |
binietoglou@0 | 49 | Update the the info, variables and dimensions of the lidar measurement based |
binietoglou@0 | 50 | on the information found in the channels. |
binietoglou@0 | 51 | |
binietoglou@0 | 52 | Reading of the scan_angles parameter is not implemented. |
binietoglou@0 | 53 | ''' |
binietoglou@0 | 54 | |
binietoglou@0 | 55 | # Initialize |
binietoglou@0 | 56 | start_time =[] |
binietoglou@0 | 57 | stop_time = [] |
binietoglou@0 | 58 | points = [] |
binietoglou@0 | 59 | all_time_scales = [] |
binietoglou@0 | 60 | channel_name_list = [] |
binietoglou@0 | 61 | |
binietoglou@0 | 62 | # Get the information from all the channels |
binietoglou@0 | 63 | for channel_name, channel in self.channels.items(): |
binietoglou@0 | 64 | channel.update() |
binietoglou@0 | 65 | start_time.append(channel.start_time) |
binietoglou@0 | 66 | stop_time.append(channel.stop_time) |
binietoglou@0 | 67 | points.append(channel.points) |
binietoglou@0 | 68 | all_time_scales.append(channel.time) |
binietoglou@0 | 69 | channel_name_list.append(channel_name) |
binietoglou@0 | 70 | |
binietoglou@0 | 71 | # Find the unique time scales used in the channels |
binietoglou@0 | 72 | time_scales = set(all_time_scales) |
binietoglou@0 | 73 | |
binietoglou@0 | 74 | # Update the info dictionary |
binietoglou@0 | 75 | self.info['start_time'] = min(start_time) |
binietoglou@0 | 76 | self.info['stop_time'] = max(stop_time) |
binietoglou@0 | 77 | self.info['duration'] = self.info['stop_time'] - self.info['start_time'] |
binietoglou@0 | 78 | |
binietoglou@0 | 79 | # Update the dimensions dictionary |
binietoglou@0 | 80 | self.dimensions['points'] = max(points) |
binietoglou@0 | 81 | self.dimensions['channels'] = len(self.channels) |
binietoglou@0 | 82 | # self.dimensions['scan angles'] = 1 |
binietoglou@0 | 83 | self.dimensions['nb_of_time_scales'] = len(time_scales) |
binietoglou@0 | 84 | |
binietoglou@0 | 85 | # Update the variables dictionary |
binietoglou@0 | 86 | # Write time scales in seconds |
binietoglou@0 | 87 | raw_Data_Start_Time = [] |
binietoglou@0 | 88 | raw_Data_Stop_Time = [] |
binietoglou@0 | 89 | |
binietoglou@0 | 90 | for current_time_scale in list(time_scales): |
binietoglou@0 | 91 | raw_start_time = np.array(current_time_scale) - min(start_time) # Time since start_time |
binietoglou@0 | 92 | raw_start_in_seconds = np.array([t.seconds for t in raw_start_time]) # Convert in seconds |
binietoglou@0 | 93 | raw_Data_Start_Time.append(raw_start_in_seconds) # And add to the list |
binietoglou@0 | 94 | # Check if this time scale has measurements every 30 or 60 seconds. |
binietoglou@0 | 95 | |
binietoglou@0 | 96 | duration = self._get_duration(raw_start_in_seconds) |
binietoglou@0 | 97 | |
binietoglou@0 | 98 | raw_stop_in_seconds = raw_start_in_seconds + duration |
binietoglou@0 | 99 | raw_Data_Stop_Time.append(raw_stop_in_seconds) |
binietoglou@0 | 100 | |
ioannis@22 | 101 | self.variables['Raw_Data_Start_Time'] = raw_Data_Start_Time |
ioannis@22 | 102 | self.variables['Raw_Data_Stop_Time'] = raw_Data_Stop_Time |
binietoglou@0 | 103 | |
binietoglou@0 | 104 | # Make a dictionary to match time scales and channels |
binietoglou@0 | 105 | channel_timescales = [] |
binietoglou@0 | 106 | for (channel_name, current_time_scale) in zip(channel_name_list, all_time_scales): |
binietoglou@0 | 107 | # The following lines are PEARL specific. The reason they are here is not clear. |
binietoglou@0 | 108 | # if channel_name =='1064BLR': |
binietoglou@0 | 109 | # channel_name = '1064' |
binietoglou@0 | 110 | for (ts,n) in zip(time_scales, range(len(time_scales))): |
binietoglou@0 | 111 | if current_time_scale == ts: |
binietoglou@0 | 112 | channel_timescales.append([channel_name,n]) |
binietoglou@0 | 113 | self.variables['id_timescale'] = dict(channel_timescales) |
binietoglou@0 | 114 | |
binietoglou@0 | 115 | def _get_duration(self, raw_start_in_seconds): |
binietoglou@0 | 116 | ''' Return the duration for a given time scale. In some files (ex. Licel) this |
binietoglou@0 | 117 | can be specified from the files themselves. In others this must be guessed. |
binietoglou@0 | 118 | |
binietoglou@0 | 119 | ''' |
binietoglou@0 | 120 | # The old method, kept here for reference |
binietoglou@0 | 121 | #dt = np.mean(np.diff(raw_start_in_seconds)) |
binietoglou@0 | 122 | #for d in duration_list: |
binietoglou@0 | 123 | # if abs(dt - d) <15: #If the difference of measuremetns is 10s near the(30 or 60) seconds |
binietoglou@0 | 124 | # duration = d |
binietoglou@0 | 125 | |
binietoglou@0 | 126 | duration = np.diff(raw_start_in_seconds)[0] |
binietoglou@0 | 127 | |
binietoglou@0 | 128 | return duration |
binietoglou@0 | 129 | |
binietoglou@0 | 130 | def subset_by_channels(self, channel_subset): |
binietoglou@0 | 131 | ''' Get a measurement object containing only the channels with names |
binietoglou@0 | 132 | contained in the channel_sublet list ''' |
binietoglou@0 | 133 | |
binietoglou@0 | 134 | m = self.__class__() # Create an object of the same type as this one. |
binietoglou@0 | 135 | m.channels = dict([(channel, self.channels[channel]) for channel |
binietoglou@0 | 136 | in channel_subset]) |
binietoglou@0 | 137 | m.update() |
binietoglou@0 | 138 | return m |
binietoglou@0 | 139 | |
binietoglou@0 | 140 | def subset_by_time(self, start_time, stop_time): |
binietoglou@0 | 141 | |
binietoglou@0 | 142 | if start_time > stop_time: |
binietoglou@0 | 143 | raise ValueError('Stop time should be after start time') |
binietoglou@0 | 144 | |
binietoglou@0 | 145 | if (start_time < self.info['start_time']) or (stop_time > self.info['stop_time']): |
binietoglou@0 | 146 | raise ValueError('The time interval specified is not part of the measurement') |
binietoglou@0 | 147 | |
binietoglou@0 | 148 | m = self.__class__() # Create an object of the same type as this one. |
binietoglou@0 | 149 | for (channel_name, channel) in self.channels.items(): |
binietoglou@0 | 150 | m.channels[channel_name] = channel.subset_by_time(start_time, stop_time) |
binietoglou@0 | 151 | m.update() |
binietoglou@0 | 152 | return m |
binietoglou@0 | 153 | |
binietoglou@0 | 154 | def r_plot(self): |
binietoglou@0 | 155 | #Make a basic plot of the data. |
binietoglou@0 | 156 | #Should include some dictionary with params to make plot stable. |
binietoglou@0 | 157 | pass |
binietoglou@0 | 158 | |
binietoglou@0 | 159 | def r_pdf(self): |
binietoglou@0 | 160 | # Create a pdf report using a basic plot and meta-data. |
binietoglou@0 | 161 | pass |
binietoglou@0 | 162 | |
binietoglou@0 | 163 | def save(self): |
binietoglou@0 | 164 | #Save the current state of the object to continue the analysis later. |
binietoglou@0 | 165 | pass |
binietoglou@0 | 166 | |
binietoglou@0 | 167 | def get_PT(self): |
binietoglou@0 | 168 | ''' Sets the pressure and temperature at station level . |
binietoglou@0 | 169 | The results are stored in the info dictionary. |
binietoglou@0 | 170 | ''' |
binietoglou@0 | 171 | |
binietoglou@0 | 172 | self.info['Temperature'] = 10.0 |
binietoglou@0 | 173 | self.info['Pressure'] = 930.0 |
binietoglou@0 | 174 | |
binietoglou@15 | 175 | def subtract_dark(self): |
binietoglou@15 | 176 | |
binietoglou@15 | 177 | if not self.dark_measurement: |
binietoglou@15 | 178 | raise IOError('No dark measurements have been imported yet.') |
binietoglou@15 | 179 | |
binietoglou@15 | 180 | for (channel_name, dark_channel) in self.dark_measurement.channels.iteritems(): |
binietoglou@15 | 181 | dark_profile = dark_channel.average_profile() |
binietoglou@15 | 182 | channel = self.channels[channel_name] |
binietoglou@15 | 183 | |
binietoglou@15 | 184 | for measurement_time, data in channel.data.iteritems(): |
binietoglou@15 | 185 | channel.data[measurement_time] = data - dark_profile |
binietoglou@15 | 186 | |
binietoglou@15 | 187 | channel.update() |
binietoglou@15 | 188 | |
ulalume3@23 | 189 | def save_as_netcdf(self, filename = None): |
binietoglou@0 | 190 | """Saves the measurement in the netcdf format as required by the SCC. |
ulalume3@23 | 191 | Input: filename. If no filename is provided <measurement_id>.nc will |
ulalume3@23 | 192 | be used. |
binietoglou@0 | 193 | """ |
binietoglou@0 | 194 | params = self.extra_netcdf_parameters |
binietoglou@0 | 195 | needed_parameters = ['Measurement_ID', 'Temperature', 'Pressure'] |
binietoglou@0 | 196 | |
binietoglou@0 | 197 | for parameter in needed_parameters: |
binietoglou@0 | 198 | stored_value = self.info.get(parameter, None) |
binietoglou@0 | 199 | if stored_value is None: |
binietoglou@0 | 200 | raise ValueError('A value needs to be specified for %s' % parameter) |
binietoglou@0 | 201 | |
ulalume3@23 | 202 | if not filename: |
ulalume3@23 | 203 | filename = "%s.nc" % self.info['Measurement_ID'] |
binietoglou@0 | 204 | |
binietoglou@0 | 205 | dimensions = {'points': 1, |
binietoglou@0 | 206 | 'channels': 1, |
binietoglou@0 | 207 | 'time': None, |
binietoglou@0 | 208 | 'nb_of_time_scales': 1, |
binietoglou@0 | 209 | 'scan_angles': 1,} # Mandatory dimensions. Time bck not implemented |
binietoglou@0 | 210 | |
binietoglou@0 | 211 | global_att = {'Measurement_ID': None, |
binietoglou@0 | 212 | 'RawData_Start_Date': None, |
binietoglou@0 | 213 | 'RawData_Start_Time_UT': None, |
binietoglou@0 | 214 | 'RawData_Stop_Time_UT': None, |
binietoglou@0 | 215 | 'RawBck_Start_Date': None, |
binietoglou@0 | 216 | 'RawBck_Start_Time_UT': None, |
binietoglou@0 | 217 | 'RawBck_Stop_Time_UT': None, |
binietoglou@0 | 218 | 'Sounding_File_Name': None, |
binietoglou@0 | 219 | 'LR_File_Name': None, |
binietoglou@0 | 220 | 'Overlap_File_Name': None, |
binietoglou@0 | 221 | 'Location': None, |
binietoglou@0 | 222 | 'System': None, |
binietoglou@0 | 223 | 'Latitude_degrees_north': None, |
binietoglou@0 | 224 | 'Longitude_degrees_east': None, |
binietoglou@0 | 225 | 'Altitude_meter_asl': None} |
binietoglou@0 | 226 | |
binietoglou@0 | 227 | channel_variables = \ |
binietoglou@0 | 228 | {'channel_ID': (('channels', ), 'i'), |
binietoglou@0 | 229 | 'Background_Low': (('channels', ), 'd'), |
binietoglou@0 | 230 | 'Background_High': (('channels', ), 'd'), |
binietoglou@0 | 231 | 'LR_Input': (('channels', ), 'i'), |
binietoglou@0 | 232 | 'DAQ_Range': (('channels', ), 'd'), |
binietoglou@0 | 233 | 'Depolarization_Factor': (('channels', ), 'd'), } |
binietoglou@0 | 234 | |
binietoglou@0 | 235 | |
binietoglou@0 | 236 | channels = self.channels.keys() |
binietoglou@0 | 237 | |
binietoglou@0 | 238 | input_values = dict(self.dimensions, **self.variables) |
binietoglou@0 | 239 | |
binietoglou@0 | 240 | # Add some mandatory global attributes |
binietoglou@0 | 241 | input_values['Measurement_ID'] = self.info['Measurement_ID'] |
ioannis@21 | 242 | input_values['RawData_Start_Date'] = self.info['start_time'].strftime('%Y%m%d') |
ioannis@21 | 243 | input_values['RawData_Start_Time_UT'] = self.info['start_time'].strftime('%H%M%S') |
ioannis@21 | 244 | input_values['RawData_Stop_Time_UT'] = self.info['stop_time'].strftime('%H%M%S') |
binietoglou@0 | 245 | |
binietoglou@0 | 246 | # Add some optional global attributes |
binietoglou@0 | 247 | input_values['System'] = params.general_parameters['System'] |
binietoglou@0 | 248 | input_values['Latitude_degrees_north'] = params.general_parameters['Latitude_degrees_north'] |
binietoglou@0 | 249 | input_values['Longitude_degrees_east'] = params.general_parameters['Longitude_degrees_east'] |
binietoglou@0 | 250 | input_values['Altitude_meter_asl'] = params.general_parameters['Altitude_meter_asl'] |
binietoglou@0 | 251 | |
binietoglou@0 | 252 | # Open a netCDF4 file |
binietoglou@0 | 253 | f = netcdf.Dataset(filename,'w', format = netcdf_format) # the format is specified in the begining of the file |
binietoglou@0 | 254 | |
binietoglou@0 | 255 | # Create the dimensions in the file |
binietoglou@0 | 256 | for (d,v) in dimensions.iteritems(): |
binietoglou@0 | 257 | v = input_values.pop(d, v) |
binietoglou@0 | 258 | f.createDimension(d,v) |
binietoglou@0 | 259 | |
binietoglou@0 | 260 | # Create global attributes |
binietoglou@0 | 261 | for (attrib,value) in global_att.iteritems(): |
binietoglou@0 | 262 | val = input_values.pop(attrib,value) |
binietoglou@0 | 263 | if val: |
ioannis@21 | 264 | setattr(f, attrib, val) |
binietoglou@0 | 265 | |
binietoglou@0 | 266 | """ Variables """ |
binietoglou@0 | 267 | # Write the values of fixes channel parameters |
binietoglou@0 | 268 | for (var,t) in channel_variables.iteritems(): |
binietoglou@0 | 269 | temp_v = f.createVariable(var,t[1],t[0]) |
binietoglou@0 | 270 | for (channel, n) in zip(channels, range(len(channels))): |
binietoglou@0 | 271 | temp_v[n] = params.channel_parameters[channel][var] |
binietoglou@0 | 272 | |
binietoglou@0 | 273 | # Write the id_timescale values |
binietoglou@0 | 274 | temp_id_timescale = f.createVariable('id_timescale','i',('channels',)) |
binietoglou@0 | 275 | for (channel, n) in zip(channels, range(len(channels))): |
binietoglou@0 | 276 | temp_id_timescale[n] = self.variables['id_timescale'][channel] |
binietoglou@0 | 277 | |
binietoglou@0 | 278 | # Laser pointing angle |
binietoglou@0 | 279 | temp_v = f.createVariable('Laser_Pointing_Angle','d',('scan_angles',)) |
binietoglou@0 | 280 | temp_v[:] = params.general_parameters['Laser_Pointing_Angle'] |
binietoglou@0 | 281 | |
binietoglou@0 | 282 | # Molecular calculation |
binietoglou@0 | 283 | temp_v = f.createVariable('Molecular_Calc','i') |
binietoglou@0 | 284 | temp_v[:] = params.general_parameters['Molecular_Calc'] |
binietoglou@0 | 285 | |
binietoglou@0 | 286 | # Laser pointing angles of profiles |
binietoglou@0 | 287 | temp_v = f.createVariable('Laser_Pointing_Angle_of_Profiles','i',('time','nb_of_time_scales')) |
binietoglou@0 | 288 | for (time_scale,n) in zip(self.variables['Raw_Data_Start_Time'], |
binietoglou@0 | 289 | range(len(self.variables['Raw_Data_Start_Time']))): |
binietoglou@0 | 290 | temp_v[:len(time_scale), n] = 0 # The lidar has only one laser pointing angle |
binietoglou@0 | 291 | |
binietoglou@0 | 292 | # Raw data start/stop time |
binietoglou@0 | 293 | temp_raw_start = f.createVariable('Raw_Data_Start_Time','i',('time','nb_of_time_scales')) |
binietoglou@0 | 294 | temp_raw_stop = f.createVariable('Raw_Data_Stop_Time','i',('time','nb_of_time_scales')) |
binietoglou@0 | 295 | for (start_time, stop_time,n) in zip(self.variables['Raw_Data_Start_Time'], |
binietoglou@0 | 296 | self.variables['Raw_Data_Stop_Time'], |
binietoglou@0 | 297 | range(len(self.variables['Raw_Data_Start_Time']))): |
binietoglou@0 | 298 | temp_raw_start[:len(start_time),n] = start_time |
binietoglou@0 | 299 | temp_raw_stop[:len(stop_time),n] = stop_time |
binietoglou@0 | 300 | |
binietoglou@0 | 301 | #Laser shots |
binietoglou@0 | 302 | temp_v = f.createVariable('Laser_Shots','i',('time','channels')) |
binietoglou@0 | 303 | for (channel,n) in zip(channels, range(len(channels))): |
binietoglou@0 | 304 | time_length = len(self.variables['Raw_Data_Start_Time'][self.variables['id_timescale'][channel]]) |
binietoglou@0 | 305 | temp_v[:time_length, n] = params.channel_parameters[channel]['Laser_Shots'] |
binietoglou@0 | 306 | |
binietoglou@0 | 307 | #Raw lidar data |
binietoglou@0 | 308 | temp_v = f.createVariable('Raw_Lidar_Data','d',('time', 'channels','points')) |
binietoglou@0 | 309 | for (channel,n) in zip(channels, range(len(channels))): |
binietoglou@0 | 310 | c = self.channels[channel] |
binietoglou@0 | 311 | temp_v[:len(c.time),n, :c.points] = c.matrix |
binietoglou@0 | 312 | |
binietoglou@0 | 313 | self.add_dark_measurements_to_netcdf(f, channels) |
binietoglou@0 | 314 | |
binietoglou@0 | 315 | #Pressure at lidar station |
binietoglou@0 | 316 | temp_v = f.createVariable('Pressure_at_Lidar_Station','d') |
binietoglou@0 | 317 | temp_v[:] = self.info['Pressure'] |
binietoglou@0 | 318 | |
binietoglou@0 | 319 | #Temperature at lidar station |
binietoglou@0 | 320 | temp_v = f.createVariable('Temperature_at_Lidar_Station','d') |
binietoglou@0 | 321 | temp_v[:] = self.info['Temperature'] |
binietoglou@0 | 322 | |
binietoglou@0 | 323 | self.save_netcdf_extra(f) |
binietoglou@0 | 324 | f.close() |
binietoglou@0 | 325 | |
binietoglou@0 | 326 | def add_dark_measurements_to_netcdf(self, f, channels): |
binietoglou@0 | 327 | |
binietoglou@0 | 328 | # Get dark measurements. If it is not given in self.dark_measurement |
binietoglou@0 | 329 | # try to get it using the get_dark_measurements method. If none is found |
binietoglou@0 | 330 | # return without adding something. |
binietoglou@0 | 331 | if self.dark_measurement is None: |
binietoglou@0 | 332 | self.dark_measurement = self.get_dark_measurements() |
binietoglou@0 | 333 | |
binietoglou@0 | 334 | if self.dark_measurement is None: |
binietoglou@0 | 335 | return |
binietoglou@0 | 336 | |
binietoglou@0 | 337 | dark_measurement = self.dark_measurement |
binietoglou@0 | 338 | |
binietoglou@0 | 339 | # Calculate the length of the time_bck dimensions |
binietoglou@0 | 340 | number_of_profiles = [len(c.time) for c in dark_measurement.channels.values()] |
binietoglou@0 | 341 | max_number_of_profiles = np.max(number_of_profiles) |
binietoglou@0 | 342 | |
binietoglou@0 | 343 | # Create the dimension |
binietoglou@0 | 344 | f.createDimension('time_bck', max_number_of_profiles) |
binietoglou@0 | 345 | |
binietoglou@0 | 346 | # Save the dark measurement data |
binietoglou@0 | 347 | temp_v = f.createVariable('Background_Profile','d',('time_bck', 'channels', 'points')) |
binietoglou@0 | 348 | for (channel,n) in zip(channels, range(len(channels))): |
binietoglou@0 | 349 | c = dark_measurement.channels[channel] |
binietoglou@0 | 350 | temp_v[:len(c.time),n, :c.points] = c.matrix |
binietoglou@0 | 351 | |
binietoglou@0 | 352 | # Dark profile start/stop time |
binietoglou@14 | 353 | temp_raw_start = f.createVariable('Raw_Bck_Start_Time','i',('time_bck','nb_of_time_scales')) |
binietoglou@14 | 354 | temp_raw_stop = f.createVariable('Raw_Bck_Stop_Time','i',('time_bck','nb_of_time_scales')) |
binietoglou@0 | 355 | for (start_time, stop_time,n) in zip(dark_measurement.variables['Raw_Data_Start_Time'], |
binietoglou@0 | 356 | dark_measurement.variables['Raw_Data_Stop_Time'], |
binietoglou@0 | 357 | range(len(dark_measurement.variables['Raw_Data_Start_Time']))): |
binietoglou@0 | 358 | temp_raw_start[:len(start_time),n] = start_time |
binietoglou@0 | 359 | temp_raw_stop[:len(stop_time),n] = stop_time |
binietoglou@0 | 360 | |
binietoglou@0 | 361 | # Dark measurement start/stop time |
binietoglou@0 | 362 | f.RawBck_Start_Date = dark_measurement.info['start_time'].strftime('%Y%m%d') |
binietoglou@0 | 363 | f.RawBck_Start_Time_UT = dark_measurement.info['start_time'].strftime('%H%M%S') |
binietoglou@0 | 364 | f.RawBck_Stop_Time_UT = dark_measurement.info['stop_time'].strftime('%H%M%S') |
binietoglou@0 | 365 | |
binietoglou@0 | 366 | |
binietoglou@0 | 367 | |
binietoglou@0 | 368 | def save_netcdf_extra(self, f): |
binietoglou@0 | 369 | pass |
binietoglou@0 | 370 | |
binietoglou@0 | 371 | def _gettime(self, date_str, time_str): |
binietoglou@0 | 372 | t = datetime.datetime.strptime(date_str+time_str,'%d/%m/%Y%H.%M.%S') |
binietoglou@0 | 373 | return t |
binietoglou@0 | 374 | |
binietoglou@0 | 375 | def plot(self): |
binietoglou@0 | 376 | for channel in self.channels: |
binietoglou@0 | 377 | self.channels[channel].plot(show_plot = False) |
binietoglou@0 | 378 | plt.show() |
binietoglou@0 | 379 | |
binietoglou@0 | 380 | def get_dark_measurements(self): |
binietoglou@0 | 381 | return None |
ioannis@22 | 382 | |
ioannis@22 | 383 | @property |
ioannis@22 | 384 | def mean_time(self): |
ioannis@22 | 385 | start_time = self.info['start_time'] |
ioannis@22 | 386 | stop_time = self.info['stop_time'] |
ioannis@22 | 387 | dt = stop_time - start_time |
ioannis@22 | 388 | t_mean = start_time + dt / 2 |
ioannis@22 | 389 | return t_mean |
binietoglou@0 | 390 | |
binietoglou@0 | 391 | |
binietoglou@0 | 392 | class Lidar_channel: |
binietoglou@0 | 393 | |
binietoglou@0 | 394 | def __init__(self,channel_parameters): |
binietoglou@0 | 395 | c = 299792458 #Speed of light |
binietoglou@0 | 396 | self.wavelength = channel_parameters['name'] |
binietoglou@0 | 397 | self.name = str(self.wavelength) |
binietoglou@0 | 398 | self.binwidth = float(channel_parameters['binwidth']) # in microseconds |
binietoglou@0 | 399 | self.data = {} |
binietoglou@0 | 400 | self.resolution = self.binwidth * c / 2 |
binietoglou@0 | 401 | self.z = np.arange(len(channel_parameters['data'])) * self.resolution + self.resolution/2.0 # Change: add half bin in the z |
binietoglou@0 | 402 | self.points = len(channel_parameters['data']) |
binietoglou@0 | 403 | self.rc = [] |
binietoglou@0 | 404 | self.duration = 60 |
binietoglou@0 | 405 | |
binietoglou@0 | 406 | def calculate_rc(self): |
binietoglou@0 | 407 | background = np.mean(self.matrix[:,4000:], axis = 1) #Calculate the background from 30000m and above |
binietoglou@0 | 408 | self.rc = (self.matrix.transpose()- background).transpose() * (self.z **2) |
binietoglou@0 | 409 | |
binietoglou@0 | 410 | |
binietoglou@0 | 411 | def update(self): |
binietoglou@0 | 412 | self.start_time = min(self.data.keys()) |
binietoglou@0 | 413 | self.stop_time = max(self.data.keys()) + datetime.timedelta(seconds = self.duration) |
binietoglou@0 | 414 | self.time = tuple(sorted(self.data.keys())) |
binietoglou@0 | 415 | sorted_data = sorted(self.data.iteritems(), key=itemgetter(0)) |
binietoglou@0 | 416 | self.matrix = np.array(map(itemgetter(1),sorted_data)) |
binietoglou@0 | 417 | |
binietoglou@0 | 418 | def _nearest_dt(self,dtime): |
binietoglou@0 | 419 | margin = datetime.timedelta(seconds = 300) |
binietoglou@0 | 420 | if ((dtime + margin) < self.start_time)| ((dtime - margin) > self.stop_time): |
binietoglou@0 | 421 | print "Requested date not covered in this file" |
binietoglou@0 | 422 | raise |
binietoglou@0 | 423 | dt = abs(self.time - np.array(dtime)) |
binietoglou@0 | 424 | dtmin = min(dt) |
binietoglou@0 | 425 | |
binietoglou@0 | 426 | if dtmin > datetime.timedelta(seconds = 60): |
binietoglou@0 | 427 | print "Nearest profile more than 60 seconds away. dt = %s." % dtmin |
binietoglou@0 | 428 | ind_t = np.where(dt == dtmin) |
binietoglou@0 | 429 | ind_a= ind_t[0] |
binietoglou@0 | 430 | if len(ind_a) > 1: |
binietoglou@0 | 431 | ind_a = ind_a[0] |
binietoglou@0 | 432 | chosen_time = self.time[ind_a] |
binietoglou@0 | 433 | return chosen_time, ind_a |
binietoglou@0 | 434 | |
binietoglou@0 | 435 | def subset_by_time(self, start_time, stop_time): |
binietoglou@0 | 436 | |
binietoglou@0 | 437 | time_array = np.array(self.time) |
binietoglou@0 | 438 | condition = (time_array >= start_time) & (time_array <= stop_time) |
binietoglou@0 | 439 | |
binietoglou@0 | 440 | subset_time = time_array[condition] |
binietoglou@0 | 441 | subset_data = dict([(c_time, self.data[c_time]) for c_time in subset_time]) |
binietoglou@0 | 442 | |
binietoglou@0 | 443 | #Create a list with the values needed by channel's __init__() |
binietoglou@0 | 444 | parameters_values = {'name': self.wavelength, |
binietoglou@0 | 445 | 'binwidth': self.binwidth, |
binietoglou@0 | 446 | 'data': subset_data[subset_time[0]],} |
binietoglou@0 | 447 | |
binietoglou@0 | 448 | c = Lidar_channel(parameters_values) |
binietoglou@0 | 449 | c.data = subset_data |
binietoglou@0 | 450 | c.update() |
binietoglou@0 | 451 | return c |
binietoglou@0 | 452 | |
binietoglou@0 | 453 | |
binietoglou@0 | 454 | def profile(self,dtime, signal_type = 'rc'): |
binietoglou@0 | 455 | t, idx = self._nearest_dt(dtime) |
binietoglou@0 | 456 | if signal_type == 'rc': |
binietoglou@0 | 457 | data = self.rc |
binietoglou@0 | 458 | else: |
binietoglou@0 | 459 | data = self.matrix |
binietoglou@0 | 460 | |
binietoglou@0 | 461 | prof = data[idx,:][0] |
binietoglou@0 | 462 | return prof, t |
binietoglou@0 | 463 | |
binietoglou@0 | 464 | def get_slice(self, starttime, endtime, signal_type = 'rc'): |
binietoglou@0 | 465 | if signal_type == 'rc': |
binietoglou@0 | 466 | data = self.rc |
binietoglou@0 | 467 | else: |
binietoglou@0 | 468 | data = self.matrix |
binietoglou@0 | 469 | tim = np.array(self.time) |
binietoglou@0 | 470 | starttime = self._nearest_dt(starttime)[0] |
binietoglou@0 | 471 | endtime = self._nearest_dt(endtime)[0] |
binietoglou@0 | 472 | condition = (tim >= starttime) & (tim <= endtime) |
binietoglou@0 | 473 | sl = data[condition, :] |
binietoglou@0 | 474 | t = tim[condition] |
binietoglou@0 | 475 | return sl,t |
binietoglou@0 | 476 | |
binietoglou@15 | 477 | def profile_for_duration(self, tim, duration = datetime.timedelta(seconds = 0), signal_type = 'rc'): |
binietoglou@15 | 478 | """ Calculates the profile around a specific time (tim). """ |
binietoglou@0 | 479 | starttime = tim - duration/2 |
binietoglou@0 | 480 | endtime = tim + duration/2 |
binietoglou@0 | 481 | d,t = self.get_slice(starttime, endtime, signal_type = signal_type) |
binietoglou@0 | 482 | prof = np.mean(d, axis = 0) |
binietoglou@0 | 483 | tmin = min(t) |
binietoglou@0 | 484 | tmax = max(t) |
binietoglou@0 | 485 | tav = tmin + (tmax-tmin)/2 |
binietoglou@0 | 486 | return prof,(tav, tmin,tmax) |
binietoglou@15 | 487 | |
binietoglou@15 | 488 | def average_profile(self): |
binietoglou@15 | 489 | """ Return the average profile (NOT range corrected) for all the duration of the measurement. """ |
binietoglou@15 | 490 | prof = np.mean(self.matrix, axis = 0) |
binietoglou@15 | 491 | return prof |
binietoglou@15 | 492 | |
binietoglou@17 | 493 | 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): |
binietoglou@0 | 494 | #if filename is not None: |
binietoglou@0 | 495 | # matplotlib.use('Agg') |
binietoglou@0 | 496 | |
binietoglou@0 | 497 | fig = plt.figure() |
binietoglou@0 | 498 | ax1 = fig.add_subplot(111) |
binietoglou@17 | 499 | self.draw_plot(ax1, cmap = cmap, signal_type = signal_type, zoom = zoom, z0 = z0, vmin = vmin, vmax = vmax) |
binietoglou@17 | 500 | |
binietoglou@17 | 501 | if title: |
binietoglou@17 | 502 | ax1.set_title(title) |
binietoglou@17 | 503 | else: |
binietoglou@17 | 504 | ax1.set_title("%s signal - %s" % (signal_type.upper(), self.name)) |
binietoglou@0 | 505 | |
binietoglou@0 | 506 | if filename is not None: |
binietoglou@0 | 507 | pass |
binietoglou@0 | 508 | #plt.savefig(filename) |
binietoglou@0 | 509 | else: |
binietoglou@0 | 510 | if show_plot: |
binietoglou@0 | 511 | plt.show() |
binietoglou@0 | 512 | #plt.close() ??? |
binietoglou@0 | 513 | |
binietoglou@17 | 514 | def draw_plot(self,ax1, cmap = plt.cm.jet, signal_type = 'rc', zoom = [0,12000,0,-1], z0 = None, cmap_label = 'a.u.', vmin = 0, vmax = 1.3 * 10 ** 7): |
binietoglou@0 | 515 | |
binietoglou@0 | 516 | if signal_type == 'rc': |
binietoglou@0 | 517 | if len(self.rc) == 0: |
binietoglou@0 | 518 | self.calculate_rc() |
binietoglou@0 | 519 | data = self.rc |
binietoglou@0 | 520 | else: |
binietoglou@0 | 521 | data = self.matrix |
binietoglou@0 | 522 | |
binietoglou@0 | 523 | hmax_idx = self.index_at_height(zoom[1]) |
binietoglou@0 | 524 | |
binietoglou@17 | 525 | # If z0 is given, then the plot is a.s.l. |
binietoglou@17 | 526 | if z0: |
binietoglou@17 | 527 | ax1.set_ylabel('Altitude a.s.l. [km]') |
binietoglou@17 | 528 | else: |
binietoglou@17 | 529 | ax1.set_ylabel('Altitude a.g.l. [km]') |
binietoglou@17 | 530 | z0 = 0 |
binietoglou@17 | 531 | |
binietoglou@17 | 532 | ax1.set_xlabel('Time UTC [hh:mm]') |
binietoglou@0 | 533 | #y axis in km, xaxis /2 to make 30s measurements in minutes. Only for 1064 |
binietoglou@0 | 534 | #dateFormatter = mpl.dates.DateFormatter('%H.%M') |
binietoglou@0 | 535 | #hourlocator = mpl.dates.HourLocator() |
binietoglou@0 | 536 | |
binietoglou@0 | 537 | #dayFormatter = mpl.dates.DateFormatter('\n\n%d/%m') |
binietoglou@0 | 538 | #daylocator = mpl.dates.DayLocator() |
binietoglou@17 | 539 | hourFormatter = mpl.dates.DateFormatter('%H:%M') |
binietoglou@19 | 540 | hourlocator = mpl.dates.AutoDateLocator(minticks = 3, maxticks = 8, interval_multiples=True) |
binietoglou@0 | 541 | |
binietoglou@0 | 542 | |
binietoglou@0 | 543 | #ax1.axes.xaxis.set_major_formatter(dayFormatter) |
binietoglou@0 | 544 | #ax1.axes.xaxis.set_major_locator(daylocator) |
binietoglou@0 | 545 | ax1.axes.xaxis.set_major_formatter(hourFormatter) |
binietoglou@0 | 546 | ax1.axes.xaxis.set_major_locator(hourlocator) |
binietoglou@0 | 547 | |
binietoglou@0 | 548 | |
binietoglou@0 | 549 | ts1 = mpl.dates.date2num(self.start_time) |
binietoglou@0 | 550 | ts2 = mpl.dates.date2num(self.stop_time) |
binietoglou@0 | 551 | |
binietoglou@0 | 552 | |
binietoglou@0 | 553 | im1 = ax1.imshow(data.transpose()[zoom[0]:hmax_idx,zoom[2]:zoom[3]], |
binietoglou@0 | 554 | aspect = 'auto', |
binietoglou@0 | 555 | origin = 'lower', |
binietoglou@0 | 556 | cmap = cmap, |
binietoglou@17 | 557 | vmin = vmin, |
binietoglou@17 | 558 | #vmin = data[:,10:400].max() * 0.1, |
binietoglou@17 | 559 | vmax = vmax, |
binietoglou@17 | 560 | #vmax = data[:,10:400].max() * 0.9, |
binietoglou@17 | 561 | extent = [ts1,ts2,self.z[zoom[0]]/1000.0 + z0/1000., self.z[hmax_idx]/1000.0 + z0/1000.], |
binietoglou@0 | 562 | ) |
binietoglou@0 | 563 | |
binietoglou@0 | 564 | cb1 = plt.colorbar(im1) |
binietoglou@17 | 565 | cb1.ax.set_ylabel(cmap_label) |
binietoglou@19 | 566 | |
binietoglou@19 | 567 | # Make the ticks of the colorbar smaller, two points smaller than the default font size |
binietoglou@19 | 568 | cb_font_size = mpl.rcParams['font.size'] - 2 |
binietoglou@19 | 569 | for ticklabels in cb1.ax.get_yticklabels(): |
binietoglou@19 | 570 | ticklabels.set_fontsize(cb_font_size) |
binietoglou@19 | 571 | cb1.ax.yaxis.get_offset_text().set_fontsize(cb_font_size) |
binietoglou@0 | 572 | |
binietoglou@0 | 573 | def index_at_height(self, height): |
binietoglou@0 | 574 | idx = np.array(np.abs(self.z - height).argmin()) |
binietoglou@0 | 575 | if idx.size >1: |
binietoglou@0 | 576 | idx =idx[0] |
binietoglou@0 | 577 | return idx |
binietoglou@0 | 578 | |
binietoglou@0 | 579 | def netcdf_from_files(LidarClass, filename, files, channels, measurement_ID): |
binietoglou@0 | 580 | #Read the lidar files and select channels |
binietoglou@0 | 581 | temp_m = LidarClass(files) |
binietoglou@0 | 582 | m = temp_m.subset_by_channels(channels) |
binietoglou@0 | 583 | m.get_PT() |
binietoglou@0 | 584 | m.info['Measurement_ID'] = measurement_ID |
binietoglou@0 | 585 | m.save_as_netcdf(filename) |
binietoglou@0 | 586 |