Mercurial > public > atmospheric_lidar / changeset

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/atmospheric_lidar/licel_depol.py	Tue Jun 07 17:06:12 2016 +0300
@@ -0,0 +1,136 @@
+import datetime
+
+import numpy as np
+
+from licel import LicelLidarMeasurement
+
+
+class LicelCalibrationMeasurement(LicelLidarMeasurement):
+
+    def __init__(self, plus45_files=None, minus45_files=None):
+        # Setup the empty class
+        super(LicelCalibrationMeasurement, self).__init__()
+
+        self.plus45_files = plus45_files
+        self.minus45_files = minus45_files
+
+        if (plus45_files and minus45_files):
+            self.check_equal_length()
+
+            self.read_channel_data()
+            self.update()
+
+    def update(self):
+        """ Correct timescales after each update.
+        """
+        super(LicelCalibrationMeasurement, self).update()
+        self.correct_timescales()
+
+    def check_equal_length(self):
+        """
+        Check if input time series have equal lengths.
+        """
+        len_plus = len(self.plus45_files)
+        len_minus = len(self.minus45_files)
+        if len_plus != len_minus:
+            raise self.UnequalMeasurementLengthError(
+                            "Input timeseries have different length: %s vs %s." % (len_plus, len_minus))
+
+    def read_channel_data(self):
+
+        # Read plus and minus 45 measurements
+        self.plus45_measurement = LicelLidarMeasurement(self.plus45_files)
+        self.plus45_measurement.rename_channel(suffix='_p45')
+
+        self.minus45_measurement = LicelLidarMeasurement(self.minus45_files)
+        self.minus45_measurement.rename_channel(suffix='_m45')
+
+        # Combine them in this object
+        self.channels = {}
+        self.channels.update(self.plus45_measurement.channels)
+        self.channels.update(self.minus45_measurement.channels)
+
+    def correct_timescales(self):
+        self.check_timescales_are_two()
+        self.combine_scales()
+
+    def check_timescales_are_two(self):
+        no_timescales = len(self.variables['Raw_Data_Start_Time'])
+        if no_timescales != 2:
+            raise self.WrongNumberOfTimescalesError("Wrong number of timescales: %s instead of 2." % no_timescales)
+
+    def combine_scales(self):
+        start_times, end_times = self.get_ordered_timescales()
+        new_start_time = start_times[0]
+        new_stop_time = end_times[1]
+        self.variables['Raw_Data_Start_Time'] = [new_start_time, ]
+        self.variables['Raw_Data_Stop_Time'] = [new_stop_time, ]
+        self.reset_timescale_id()
+
+    def reset_timescale_id(self):
+        """
+        Set all timescales to 0
+        :return:
+        """
+        timescale_dict = self.variables['id_timescale']
+        self.variables['id_timescale'] = dict.fromkeys(timescale_dict, 0)
+
+    def get_ordered_timescales(self):
+        scale_start_1, scale_start_2 = self.variables['Raw_Data_Start_Time']
+        scale_end_1, scale_end_2 = self.variables['Raw_Data_Stop_Time']
+
+        if scale_start_1[0] > scale_start_2[0]:
+            scale_start_1, scale_start_2 = scale_start_2, scale_start_1
+
+        if scale_end_1[0] > scale_end_2[0]:
+            scale_end_1, scale_end_2 = scale_end_2, scale_end_1
+
+        return (scale_start_1, scale_start_2), (scale_end_1, scale_end_2)
+
+    def add_fake_measurements(self, no_profiles, variation=0.1):
+        """
+        Add a number of fake measurements. This is done to allow testing with single analog profiles.
+
+        Adds a predefined variation in each new profile.
+        """
+        duration = self.info['duration']
+        for channel_name, channel in self.channels.items():
+            base_time = channel.data.keys()[0]
+            base_data = channel.data[base_time]
+
+            for n in range(no_profiles):
+                random_variation = base_data * (np.random.rand(len(base_data)) * 2 - 1) * variation
+
+                new_time = base_time + n * duration
+                new_data = channel.data[base_time].copy() + random_variation
+                if 'ph' in channel_name:
+                    new_data = new_data.astype('int')
+                channel.data[new_time] = new_data
+
+        self.update()
+
+    def subset_by_netcdf_channels(self):
+        """
+        Subset the measurement based on the available netcdf channels in the parameters file.
+        """
+        channels = self.extra_netcdf_parameters.channel_parameters.keys()
+        new_measurement = self.subset_by_channels(channels)
+        return new_measurement
+
+    def subset_photoncounting(self):
+        """
+        Subset photoncounting channels.
+        """
+        ph_channels = [channel for channel in self.channels.keys() if 'ph' in channel]
+        new_measurement = self.subset_by_channels(ph_channels)
+        return new_measurement
+
+    class UnequalMeasurementLengthError(RuntimeError):
+        """ Raised when the plus and minus files have different length.
+        """
+        pass
+
+    class WrongNumberOfTimescalesError(RuntimeError):
+        """ Raised when timescales are not two.
+        """
+        pass

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/atmospheric_lidar/rali_depol.py	Tue Jun 07 17:06:12 2016 +0300
@@ -0,0 +1,5 @@
+from licel_depol import LicelCalibrationMeasurement
+import rali_depolarization_parameters
+
+class RALICalibrationMeasurement(LicelCalibrationMeasurement):
+    extra_netcdf_parameters = rali_depolarization_parameters

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/atmospheric_lidar/rali_depolarization_parameters.py	Tue Jun 07 17:06:12 2016 +0300
@@ -0,0 +1,83 @@
+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 = \
+{ '00532.p_an_p45': {'channel_ID': 383,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 100.0,
+         'Depolarization_Factor': 0,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.p_ph_p45': {'channel_ID': 378,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 0,
+         'Depolarization_Factor': 0,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.s_an_p45': {'channel_ID': 385 ,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 20.0,
+         'Depolarization_Factor': 0.115,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.s_ph_p45': {'channel_ID': 380,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 0,
+         'Depolarization_Factor': 0.115,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.p_an_m45': {'channel_ID': 384,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 100.0,
+         'Depolarization_Factor': 0,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.p_ph_m45': {'channel_ID': 379,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 0,
+         'Depolarization_Factor': 0,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.s_an_m45': {'channel_ID': 386,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 20.0,
+         'Depolarization_Factor': 0.115,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+ '00532.s_ph_m45': {'channel_ID': 382,
+         'Background_Low': 50000,
+         'Background_High': 60000,
+         'Laser_Shots': 3000,
+         'LR_Input': 1,
+         'DAQ_Range': 0,
+         'Depolarization_Factor': 0.115,
+         'Pol_Calib_Range_Min': 1000,
+         'Pol_Calib_Range_Max': 3000},
+         }
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/readme.rst	Tue Jun 07 17:06:12 2016 +0300
@@ -0,0 +1,78 @@
+Basic instructions
+==================
+
+These classes can be used to handle lidar input data. They are still in a very initial stage. There are many features probably not working, but they work for some specific tasks. They work with Licel input files (and also with the Raymetrics modified format).
+
+
+Set up
+------
+
+### Parameter file
+
+Before using the classes you need to setup some channel parameters, that are used when converting the lidar data to Single Calculus Chain format. 
+
+All the parameters are read from an external file stored in the same folder as the code. You can start by changing the file “cf_netcdf_parameters.py” that describe such parameters for the Clermont Ferrand  lidar. 
+
+### System class
+
+The next thing you need to create a class that describes you system. This is very simple if your lidar data are in the Licel format, as you only need to specify the external file with the extra SCC parameters. You can use as an example the file “cf_raymetrics.py”:
+   
+    :::python
+    from licel import LicelLidarMeasurement
+    import cf_netcdf_parameters
+
+    class CfLidarMeasurement(LicelLidarMeasurement):
+       extra_netcdf_parameters = cf_netcdf_parameters
+
+
+Using the class
+---------------
+
+Once you have made the above setup you can start using it. The best way to understand how it works is through an interactive shell (I suggest [ipython](http://ipython.org/)).  In the following example I use the cf_raymetrics setup:
+   
+    :::python
+    import glob # This is needed to read a list of filenames
+    from lidar import cf_raymetrics #If you have saved the files in a directrory called “lidar”
+
+    # Go to the folder where you files are stored
+    cd /path/to/lidar/files
+
+    # Read the filenames
+    files  = glob.glob("*") # The * reads all the files in the folder.
+
+    #Read the files
+    my_measurement = cf_raymetrics.CfLidarMeasurement(files)
+
+    #Now the data have been read, and you have a measurement object to work with:
+    # See what channels are present
+    print my_measurement.channels
+
+    #Quicklooks of all the channels
+    my_measurements.plot()
+
+
+Converting to SCC format
+------------------------
+
+There are some extra info you need to put in before converting to SCC format, "Measurement_ID", "Temperature", "Pressure":
+   
+    :::python
+    my_measurement.info["Measurement_ID"] = "20101229op00"
+    my_measurement.info["Temperature"] = "14"
+    my_measurement.info["Pressure"] = "1010"
+
+You can use standard values of temperature and pressure by just calling:
+   
+    :::python
+    my_measurement.get_PT() 
+
+The standard values can be changed in generic.py. Search the  get_PT method and change of what is appropriate for your station. If you have an external source of temperature and pressure information (a meteorological station) you can automate this by overriding the get_PT method in your system"s class (in our example in the cf_raymetrics.py file).
+
+
+After you have used this extra input, you save the file using this command:
+
+    :::python
+    my_measurement.save_as_netcdf("filename")
+
+where you change the filename to the filename you want to use.
+