# Copyright (c) 2011-2014, B.I.Stepanov Institute of Physics, National Academy # of Sciences of Belarus. # # This program is free software; you can redistribute it and/or modify it under # the terms of the GNU General Public License as published by the Free Software # Foundation; either version 2 of the License, or (at your option) any later # version. # # This program is distributed in the hope that it will be useful, but WITHOUT # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. import datetime import math import numpy import numpy.linalg import os import os.path from common.utils import txt from common.DataRecord import * __all__ = ['PhotometerInput'] # ***************************************************************************** class PhotometerInput(DataRecord): """Sun photometer input data for the aerosol profile retrieval algorithm. Attributes to be read from an Excel file: - 'date', 'time': measurement date and time. - 'aot675': total aerosol optical thickness at 675 nanometers. - 'aot675Corr': manually corrected value for 'aot675'. The same correction (with correction factor defined as 'aot675Corr / aot675') has to be applied to 'vFine', 'vCoarse' and 'ext<...>'. - 'vFine', 'vCoarse': volumetric concentrations of fine and coarse aerosol modes, as retrieved with Sun photometer, in micrometers (with both sperical and spheroidal particles taken into account). - 'sphericity': fraction of spherical aerosol particles, in percents, in contrast to spheroidal particles, as retrieved with Sun photometer. - 'bFine<...>', bSpherical<...>', 'bSpheroid<...>', bCoarse<...>: aerosol backscatter coefficients for different aerosol modes at differnt wavelengths and polarizations. Cross-polarized and parallel-polarized backscatters are denoted with additional 'C' and 'P' suffixes respectively; total backscatters are denoted with just the wavelength in nanometers. 'Fine' aerosol mode includes both spherical and spheroidal particles; 'Coarse' mode is a sum of 'Spherical' and 'Spheroid' modes. - 'aFine<...>', 'aSpherical<...>', 'aSpheroid<...>', aCoarse<...>: aerosol attenuation coefficients for different aerosol modes at different wavelengths. - 'wl<...>Um': wavelengths, at which recalculated aerosol parameters are given, in micrometers. 'wl355Um' has to be equal to '0.355', and so on. - 're<...>': real part of the complex refractive index at the given wavelength. - 'im<...>': imaginary part of the complex refractive index at the given wavelength. - 'ext<...>': aerosol optical thickness at the given wavelength. - 'ssa<...>': single scattering albedo at the given wavelength. - 'f11_<...>': main component of the aerosol scattering matrix at the given wavelength. - 'f22by11_<...>': second diagonal element of the aerosol scattering matrix divided by its main component, at the given wavelength. Attributes that are calculated manually in 'onDataLoaded': - 'vFineCorr', 'vCoarseCorr': corrected values for 'vFine' and 'vCoarse' (see 'aot675Corr'). - 'vSphericalCorr', 'vSpheroidCorr': the value of 'vCoarseCorr', divided into two parts according to the value of 'sphericity' attribute. - 'wl<...>': same as 'wl<...>Um', but converted into nanometers. - 'ext<...>Corr': corrected values for 'ext<...>' (see 'aot675Corr'). - 'evRatio2Modes', 'evRatio3Modes': Auxiliary attribute (used in the photometer input widget only): - 'aot675CorrOriginal': the stored value for 'aot675Corr', in contrast to the value modified by the user interactively.""" # ---- Class attributes --------------------------------------------------- # Name of the Excel file worksheet that holds the data. tableName = 'Photometer Data-ver2' fields = [ DataField('date', 'Date', 'DATE'), DataField('time', 'Time', 'TIME'), DataField('aot675', 'AOT-675', 'DOUBLE'), DataField('aot675Corr', 'AOT_675-TRU', 'DOUBLE'), DataField('vFine', '1V1', 'DOUBLE'), DataField('vCoarse', '2V3', 'DOUBLE'), DataField('sphericity', 'T1Sphericity', 'DOUBLE'), DataField('bFine355', 'b1,0,355', 'DOUBLE'), DataField('bFine355P', 'b1,3,355', 'DOUBLE'), DataField('bFine355C', 'b1,2,355', 'DOUBLE'), DataField('bFine532', 'b1,0,532', 'DOUBLE'), DataField('bFine532P', 'b1,3,532', 'DOUBLE'), DataField('bFine532C', 'b1,2,532', 'DOUBLE'), DataField('bFine1064', 'b1,0,1064', 'DOUBLE'), DataField('bFine1064P', 'b1,3,1064', 'DOUBLE'), DataField('bFine1064C', 'b1,2,1064', 'DOUBLE'), DataField('bSpherical355', 'b2,0,355', 'DOUBLE'), DataField('bSpherical355P', 'b2,3,355', 'DOUBLE'), DataField('bSpherical355C', 'b2,2,355', 'DOUBLE'), DataField('bSpherical532', 'b2,0,532', 'DOUBLE'), DataField('bSpherical532P', 'b2,3,532', 'DOUBLE'), DataField('bSpherical532C', 'b2,2,532', 'DOUBLE'), DataField('bSpherical1064', 'b2,0,1064', 'DOUBLE'), DataField('bSpherical1064P', 'b2,3,1064', 'DOUBLE'), DataField('bSpherical1064C', 'b2,2,1064', 'DOUBLE'), DataField('bSpheroid355', 'b3,0,355', 'DOUBLE'), DataField('bSpheroid355P', 'b3,3,355', 'DOUBLE'), DataField('bSpheroid355C', 'b3,2,355', 'DOUBLE'), DataField('bSpheroid532', 'b3,0,532', 'DOUBLE'), DataField('bSpheroid532P', 'b3,3,532', 'DOUBLE'), DataField('bSpheroid532C', 'b3,2,532', 'DOUBLE'), DataField('bSpheroid1064', 'b3,0,1064', 'DOUBLE'), DataField('bSpheroid1064P', 'b3,3,1064', 'DOUBLE'), DataField('bSpheroid1064C', 'b3,2,1064', 'DOUBLE'), DataField('bCoarse355', 'b4,0,355', 'DOUBLE'), DataField('bCoarse355P', 'b4,3,355', 'DOUBLE'), DataField('bCoarse355C', 'b4,2,355', 'DOUBLE'), DataField('bCoarse532', 'b4,0,532', 'DOUBLE'), DataField('bCoarse532P', 'b4,3,532', 'DOUBLE'), DataField('bCoarse532C', 'b4,2,532', 'DOUBLE'), DataField('bCoarse1064', 'b4,0,1064', 'DOUBLE'), DataField('bCoarse1064P', 'b4,3,1064', 'DOUBLE'), DataField('bCoarse1064C', 'b4,2,1064', 'DOUBLE'), DataField('aFine355', 'a1,0,355', 'DOUBLE'), DataField('aFine532', 'a1,0,532', 'DOUBLE'), DataField('aFine1064', 'a1,0,1064', 'DOUBLE'), DataField('aSpherical355', 'a2,0,355', 'DOUBLE'), DataField('aSpherical532', 'a2,0,532', 'DOUBLE'), DataField('aSpherical1064', 'a2,0,1064', 'DOUBLE'), DataField('aSpheroid355', 'a3,0,355', 'DOUBLE'), DataField('aSpheroid532', 'a3,0,532', 'DOUBLE'), DataField('aSpheroid1064', 'a3,0,1064', 'DOUBLE'), DataField('aCoarse355', 'a4,0,355', 'DOUBLE'), DataField('aCoarse532', 'a4,0,532', 'DOUBLE'), DataField('aCoarse1064', 'a4,0,1064', 'DOUBLE'), DataField('wl355Um', 'T1Wave', 'DOUBLE'), DataField('re355', 'T1REFR', 'DOUBLE'), DataField('im355', 'T1REFI', 'DOUBLE'), DataField('ext355', 'T1Ext', 'DOUBLE'), DataField('ssa355', 'T1Albedo', 'DOUBLE'), DataField('f11_355', 'T1F11', 'DOUBLE'), DataField('f22by11_355', 'T1F22/F11', 'DOUBLE'), DataField('wl532Um', 'T2Wave', 'DOUBLE'), DataField('re532', 'T2REFR', 'DOUBLE'), DataField('im532', 'T2REFI', 'DOUBLE'), DataField('ext532', 'T2Ext', 'DOUBLE'), DataField('ssa532', 'T2Albedo', 'DOUBLE'), DataField('f11_532', 'T2F11', 'DOUBLE'), DataField('f22by11_532', 'T2F22/F11', 'DOUBLE'), DataField('wl1064Um', 'T3Wave', 'DOUBLE'), DataField('re1064', 'T3REFR', 'DOUBLE'), DataField('im1064', 'T3REFI', 'DOUBLE'), DataField('ext1064', 'T3Ext', 'DOUBLE'), DataField('ssa1064', 'T3Albedo', 'DOUBLE'), DataField('f11_1064', 'T3F11', 'DOUBLE'), DataField('f22by11_1064', 'T3F22/F11', 'DOUBLE'), # These fields are only used in Excel file export procedure. DataField('evRatio2Modes', 'EVRatio2', 'DOUBLE', required = False), DataField('evRatio3Modes', 'EVRatio3', 'DOUBLE', required = False), ] extraExcelHeaderRows = 1 openExcelFileErrorMessage = ( 'Photometer input is not a valid Excel file or may not be opened for ' 'reading') invalidExcelFormatErrorPrefix = ( 'Photometer input file format is invalid: ') exportDataToExcelErrorPrefix = ( 'Failed to write the data to the Excel file') # ---- Public methods ----------------------------------------------------- def onDataLoaded(self): # Convert micrometers to nanometers to improve readability and make # lidar and photometer data uniform. self.wl355 = self.getConvertedWavelength(self.wl355Um) self.wl532 = self.getConvertedWavelength(self.wl532Um) self.wl1064 = self.getConvertedWavelength(self.wl1064Um) # Create and initialize all the attributes for manually corrected # thicknesses and concentrations. self.updateCorrectedThicknesses() # Keep original value of 'aot675Corr' so that it would be possible to # say if the value has been modified interactively by the user (this is # used by 'PhotometerInputWidget'). self.aot675CorrOriginal = self.aot675Corr self.calcEigenvalueRatios() def getDateTime(self): if self.date is None or self.time is None: return None else: return datetime.datetime.combine(self.date, self.time) def updateCorrectedThicknesses(self): """Set values for all the attributes that depend on the manually corrected extinction aerosol optical thickness at 675 nm. Call this function whenever 'aot675Corr' attribute is modified.""" self.vFineCorr = self.getCorrectedThickness(self.vFine) self.vCoarseCorr = self.getCorrectedThickness(self.vCoarse) if self.vCoarseCorr is not None and self.sphericity is not None: self.vSphericalCorr = self.vCoarseCorr * (self.sphericity * 0.01) self.vSpheroidCorr = self.vCoarseCorr * ( 1.0 - self.sphericity * 0.01) else: self.vSphericalCorr = None self.vSpheroidCorr = None self.ext355Corr = self.getCorrectedThickness(self.ext355) self.ext532Corr = self.getCorrectedThickness(self.ext532) self.ext1064Corr = self.getCorrectedThickness(self.ext1064) def modifyPhaseFunction(self, modeSuffix, wlSuffix, f11Factor, f22By11Factor): """Modify the values of aerosol backscatter coefficients so that they would mimic the specified modifications of the aerosol scattering phase matrix for the given aerosol mode and at the given wavelength. Parameters: - 'modeSuffix', 'wlSuffix': strings identifying the aerosol mode and the wavelength for which the phase matrix is being modified. 'modeSuffix' may be 'Fine', 'Spherical', 'Spheroid', or 'Coarse'; 'wlSuffix' may be '355', '532', or '1064'. - 'f11Factor', 'f22By11Factor': ratios between the new and the old values of F11 and F22/F11 components of the scattering matrix, respectively. Warning: currently, this function does not make any attempts to modify 'f11_<...>' and 'f22by11_<...>' attributes of the class instance, thus effectively breaking its internal consistency. After this function has has been called, core aerosol parameters and utility attributes used to describe the aerosol as a whole won't agree any more.""" # Read current backscatter values for the given mode and wavelength. oldBscat = getattr(self, 'b' + modeSuffix + wlSuffix) oldBscatP = getattr(self, 'b' + modeSuffix + wlSuffix + 'P') oldBscatC = getattr(self, 'b' + modeSuffix + wlSuffix + 'C') # Derive the current phase matrix coefficient from the backscatters. oldF22By11 = (oldBscatP - oldBscatC) / oldBscat bscat = oldBscat * f11Factor f22By11 = oldF22By11 * f22By11Factor bscatC = bscat * (1 - f22By11) * 0.5 bscatP = bscat * (1 + f22By11) * 0.5 # Attenuation coefficients and optical thicknesses won't change when # the phase function is modified. setattr(self, 'b' + modeSuffix + wlSuffix, bscat) setattr(self, 'b' + modeSuffix + wlSuffix + 'P', bscatP) setattr(self, 'b' + modeSuffix + wlSuffix + 'C', bscatC) @staticmethod def readDataListFromAerlid(aerlidFilePath, aotCorrFilePath = None): """Same as 'readDataList', but read the data from an AERLID data file (i.e., from a Dubovik AERONET data file processed with AERLID tool). If 'aotCorrFilePath' is not 'None', it should be a path to the corresponding AOT correction data file (i.e., an auxiliary data file that holds manually corrected values of aerosol optical thicknesses at 675 nm).""" recordList = AerlidFileLoader().loadFile(aerlidFilePath) if aotCorrFilePath is not None: aot675CorrDict = AotCorrFileLoader().loadFile(aotCorrFilePath) for dataRecord in recordList: # Modify the 'dataRecord's 'aot675Corr' attribute if the value # for it was specified in the AOT correction data file. dateTimeValue = datetime.datetime.combine( dataRecord.date, dataRecord.time) if dateTimeValue in aot675CorrDict: dataRecord.aot675Corr = aot675CorrDict[dateTimeValue] # Update dependent instance attributes. dataRecord.updateCorrectedThicknesses() return recordList @staticmethod def updateAotCorrFileForDataList(aotCorrFilePath, recordList): """Write AOT correction data file for the given list of 'PhotometerInput' instances. Only those corrected AOT values that are different from the original ones are stored in the 'AOT corr' data file. Thus, if the original data file changes, its new values for 'aot675' won't get overwritten by the corrected values stored in the 'AOT corr' file. If 'recordList' contains no data records with user-modified AOT-675 values, 'AOT corr' data file is not created, and if there already exists one, it is deleted. Thus, an AOT correction data file will never exist unless it is really required.""" fileIsRequired = any(dataRecord.aot675Corr != dataRecord.aot675 for dataRecord in recordList) if not fileIsRequired: # Delete the existing file, if any, and do nothing more. fileExists = os.path.isfile(aotCorrFilePath) if fileExists: try: os.remove(aotCorrFilePath) except OSError: raise txt.Error( u'Couldn\u2019t delete the outdated AOT correction ' 'file (%s)' % txt.quotePath(aotCorrFilePath)) return fileDir = os.path.dirname(aotCorrFilePath) # Create destination directory, if it does not exist. try: # Using 'isdir' instead of 'exists' will make the code fail if the # required directory does not exist and may not be created due to a # naming conflict with an existing ordinary file. if not os.path.isdir(fileDir): os.makedirs(fileDir) except OSError: raise txt.Error( u'Couldn\u2019t create directory for the AOT correction file ' '(%s)' % txt.quotePath(aotCorrFilePath)) try: aotCorrFile = open(aotCorrFilePath, 'w') except IOError: raise txt.Error( 'AOT correction file %s may not be opened for writing' % txt.quotePath(aotCorrFilePath)) try: # Write the file header. Only date, time, and the corrected value # for AOT-675 are stored in the AOT correction data file. aotCorrFile.write('yyyy-mm-dd hh:mm:ss AOT_675-TRU\n') for dataRecord in recordList: # Skip values that are not different from the original ones. if dataRecord.aot675Corr != dataRecord.aot675: # Use the format that is used in AERLID data files. dateStr = dataRecord.date.strftime('%Y-%m-%d ') timeStr = dataRecord.time.strftime('%H:%M:%S ') # Use the same precision that is used in AERLID data files # (actually, the precision will be even higher, because the # most significant digit of AERLID's floating point values # is always zero, unlike the Python ones). aot675CorrStr = '% .5E' % dataRecord.aot675Corr aotCorrFile.write(dateStr + timeStr + aot675CorrStr + '\n') except IOError: raise txt.Error('Failed to write AOT correction file (%s)' % txt.quotePath(aotCorrFilePath)) finally: aotCorrFile.close() # ---- Private methods ---------------------------------------------------- def getConvertedWavelength(self, wavelenMicrometers): if wavelenMicrometers is None: # Conversion is not possible. return None return wavelenMicrometers * 1000.0 def getCorrectedThickness(self, originalThickness): """Correct an abstract 'originalThickness' quantity according to difference between original and manually corrected aerosol optical thicknesses at 675 nm. The quantity is assumed to be proportional to the optical thicnkess.""" if (self.aot675 is None or self.aot675 is None or originalThickness is None): # Correction is not possible. return originalThickness if self.aot675 == 0.0: # Correction is not possible. return originalThickness else: correctedValue = originalThickness * self.aot675Corr / self.aot675 if math.isinf(correctedValue): # Just to feel safe. return originalThickness return correctedValue @staticmethod def getEVRatio(matrix): eigenValues = numpy.linalg.eigvalsh(matrix.transpose() * matrix) numpy.sort(eigenValues) return eigenValues[-1] / eigenValues[0] def calcEigenvalueRatios(self): self.evRatio2Modes = self.getEVRatio(numpy.matrix([ [self.bFine355, self.bCoarse355], [self.bFine532, self.bCoarse532], [self.bFine1064, self.bCoarse1064]])) self.evRatio3Modes = self.getEVRatio(numpy.matrix([ [self.bFine355, self.bSpherical355, self.bSpheroid355], [self.bFine532, self.bSpherical532, self.bSpheroid532], [self.bFine1064, self.bSpherical1064, self.bSpheroid1064], [self.bFine532C, self.bSpherical532C, self.bSpheroid532C]])) # ***************************************************************************** class AerlidFileLoader: """Utility class used for parsing of AERLID photometer data files.""" # Private attributes: # - 'file': instance of the opened data file that is being loaded. # - 'errorPrefix': string to be used with error messages related to # invalid format of the data file. # - 'currentLine': zero-based index of the data file line being # processed. # ---- Public methods ----------------------------------------------------- def __init__(self): """Initialize auxiliary attributes used to describe the data format.""" # Data are stored in a text file with a large number of data columns. # Names of aerosol modes that should be present in the data file. self.aerosolModes = ('Total', 'Fine_Total', 'Fine_Spherical', 'Fine_Non-Spherical', 'Coarse_Total', 'Coarse_Spherical', 'Coarse_Non-Spherical') # Wavelength suffixes for 'PhotometerInput' instance attributes. self.wavelengthSuffixes = ('355', '532', '1064') # Corresponding wavelength values (used for data format checking). self.wavelengthsUm = (0.355, 0.532, 1.064) assert len(self.wavelengthSuffixes) == len(self.wavelengthsUm) # Headers of the initial data columns (not bound to any aerosol mode). self.prefixHeaders = ('yyyy-mm-dd', 'hh:mm:ss', 'Julian_day') # Headers of the final data columns (not bound to any aerosol mode). self.suffixHeaders = ('V-fine', 'V-coarse', 'AOT_675', 'AOT_675-TRU') # Headers of data columns that describe properties of a given aerosol # mode at a given wavelength. Names of aerosol modes are defined by # 'self.aerosolModes'; wavelengths are defined by 'self.wavelengthsUm'. self.specificHeaders = ('Wave', 'REFR', 'REFI', '%sphericity', 'ext', 'abs', 'sca', 'albedo', 'F11', 'F22/F11') # Indices of the required aerosol modes. self.totalModeIndex = self.aerosolModes.index('Total') self.fineModeIndex = self.aerosolModes.index('Fine_Total') self.sphericalModeIndex = self.aerosolModes.index('Coarse_Spherical') self.spheroidModeIndex = self.aerosolModes.index( 'Coarse_Non-Spherical') self.coarseModeIndex = self.aerosolModes.index('Coarse_Total') # Indices of the required generic data columns. self.dateCol = self.prefixHeaders.index('yyyy-mm-dd') self.timeCol = self.prefixHeaders.index('hh:mm:ss') # Indices of the final data columns are denoted with negative numbers. self.vFineCol = (self.suffixHeaders.index('V-fine') - len(self.suffixHeaders)) self.vCoarseCol = (self.suffixHeaders.index('V-coarse') - len(self.suffixHeaders)) self.aot675Col = (self.suffixHeaders.index('AOT_675') - len(self.suffixHeaders)) self.aot675CorrCol = (self.suffixHeaders.index('AOT_675-TRU') - len(self.suffixHeaders)) # Relative indices of data columns bound to a given aerosol mode and # wavelength. self.wlCol = self.specificHeaders.index('Wave') self.reCol = self.specificHeaders.index('REFR') self.imCol = self.specificHeaders.index('REFI') self.sphericityCol = self.specificHeaders.index('%sphericity') self.extCol = self.specificHeaders.index('ext') self.scatteringCol = self.specificHeaders.index('sca') self.ssaCol = self.specificHeaders.index('albedo') self.f11Col = self.specificHeaders.index('F11') self.f22by11Col = self.specificHeaders.index('F22/F11') def loadFile(self, filePath): """Load the specified AERLID data file. Return a list of 'PhotometerInput' instances. Raise 'txt.Error' if loading fails or file format is invalid.""" # The real loading is done in 'readFile'. Here, the file is opened, # and 'file' and 'errorPrefix' attributes are set. try: # If a Dubovik AERONET file contains no records, the corresponding # AERLID file will contain no data altogether, and its format # checking will fail. Treat this special case separately. if os.path.getsize(filePath) == 0: return [] else: self.file = open(filePath) # Include 'OSError' to catch exceptions from 'os.path.getsize'. except (IOError, OSError): raise txt.Error( 'AERLID data file %s may not be opened for reading' % txt.quotePath(filePath)) self.errorPrefix = ('AERLID data file format is invalid (%s): ' % txt.quotePath(filePath)) try: recordList = self.readFile() except IOError: raise txt.Error('Failed to read AERLID data file (%s)' % txt.quotePath(filePath)) finally: self.file.close() return recordList # ---- Private methods ---------------------------------------------------- def getColumnIndex(self, col, modeIndex = None, wlIndex = None): """Calculate the absolute index of a data column based on its relative index 'col'. If 'modeIndex' and 'wlIndex' are not 'None', they should be valid indices within 'self.aerosolModes' and 'self.wavelenghtsUm' arrays respectively, and 'col' is considered to be a relative data column index for a given aerosol mode and wavelength. Otherwise, 'col' is considered to refer to either initial ('col' >= 0) or final ('col' < 0) generic data column (thus, 'col' == 0 will refer to the first data column within the file, 'col' == -1 will refer to the last data column, 'col' == -2 will refer to the last but one column, and so on).""" if modeIndex is not None: # Specific headers are grouped by aerosol modes, and by wavelength # within a single mode. absIndex = (len(self.prefixHeaders) + len(self.specificHeaders) * (modeIndex * len(self.wavelengthsUm) + wlIndex) + col) elif col < 0: # Add the negative 'col' value to the total number of the columns. absIndex = (len(self.prefixHeaders) + len(self.specificHeaders) * (len(self.aerosolModes) * len(self.wavelengthsUm)) + len(self.suffixHeaders) + col) else: absIndex = col assert absIndex >= 0 return absIndex def readFile(self): """Load the opened data file referenced in 'self.file'.""" # This will be a list of 'PhotometerInput' instances. recordList = [] # Skip the AERONET site location header. self.file.readline() # Read the aerosol mode headers and check if all the modes are present. # Column headers are divided by whitespace in AERLID files. modeHeaders = self.file.readline().split() for modeIndex in range(len(self.aerosolModes)): if (modeIndex >= len(modeHeaders) or modeHeaders[modeIndex] != self.aerosolModes[modeIndex]): raise txt.Error(self.errorPrefix + '%s aerosol mode is missing' % txt.quote(self.aerosolModes[modeIndex])) # Read the data column headers and check them against the expected # values. Column headers are divided by whitespace in AERLID files. dataHeaders = self.file.readline().split() for col in range(len(self.prefixHeaders)): self.testColumnHeader(dataHeaders, self.prefixHeaders[col], col) for modeIndex in range(len(self.aerosolModes)): for wlIndex in range(len(self.wavelengthsUm)): for col in range(len(self.specificHeaders)): self.testColumnHeader(dataHeaders, self.specificHeaders[col], col, modeIndex, wlIndex) for col in range(-len(self.suffixHeaders), 0): self.testColumnHeader(dataHeaders, self.suffixHeaders[col], col) # Initialize the 'currentLine' attribute (there are 3 header lines). self.currentLine = 3 for line in self.file: # Read the data strings stored in the current line. # Data columns are divided by whitespace in AERLID files. dataStrings = line.split() dataRecord = self.readDataRecord(dataStrings) if dataRecord is not None: # Post-process the loaded data. dataRecord.onDataLoaded() recordList.append(dataRecord) # Update the 'currentLine' attribute. self.currentLine += 1 return recordList def readDataRecord(self, dataStrings): """Parse a data record that has been recently read from the data file. Return a valid 'PhotometerInput' instance or 'None' if some of the required data fields are invalid (i.e., contain the values of -999.0). 'dataList' should be the comlete list of data strings read from a single line of the file.""" dataRecord = PhotometerInput() # Read contents of the generic data columns. dateTimeValue = self.readDateTimeValue(dataStrings) dataRecord.date = dateTimeValue.date() dataRecord.time = dateTimeValue.time() dataRecord.aot675 = self.readFloatValue(dataStrings, self.aot675Col) dataRecord.aot675Corr = self.readFloatValue(dataStrings, self.aot675CorrCol) dataRecord.vFine = self.readFloatValue(dataStrings, self.vFineCol) dataRecord.vCoarse = self.readFloatValue(dataStrings, self.vCoarseCol) # Aerosol sphericity should be the same for all the wavelengths and # for all the 'Total' aerosol modes (for 'Spherical' modes it is # always 100, and for 'Non-Spherical' modes it is always 0). Read # the value stored in the most generic aerosol mode, using the # wavelength with zero index. dataRecord.sphericity = self.readFloatValue(dataStrings, self.sphericityCol, self.totalModeIndex, 0) # AERLID may produce incorrect results if some of AERONET input # fields are not available (i.e., if they have the value of 'N/A', # which is replaced with -999.0 during the preprocessing step). # Discars the incomplete (and apparently invalid) data records. if any(dataField == -999.0 for dataField in (dataRecord.aot675, dataRecord.aot675Corr, dataRecord.vFine, dataRecord.vCoarse, dataRecord.sphericity)): return None for wlIndex in range(len(self.wavelengthsUm)): # Read contents of data columns bound to the total aerosol mode # and to the given wavelength. wlSuffix = self.wavelengthSuffixes[wlIndex] wlValueUm = self.readFloatValue(dataStrings, self.wlCol, self.totalModeIndex, wlIndex) # Make sure that data stored in the file represent the required # wavelengths. if wlValueUm != self.wavelengthsUm[wlIndex]: raise txt.Error( 'wavelength value is %s instead of %s in line %s' % (txt.quoteNumber(wlValueUm), txt.quoteNumber(self.wavelengthsUm[wlIndex]), txt.quoteNumber(self.currentLine + 1))) setattr(dataRecord, 'wl' + wlSuffix + 'Um', wlValueUm) setattr(dataRecord, 're' + wlSuffix, self.readFloatValue( dataStrings, self.reCol, self.totalModeIndex, wlIndex)) setattr(dataRecord, 'im' + wlSuffix, self.readFloatValue( dataStrings, self.imCol, self.totalModeIndex, wlIndex)) setattr(dataRecord, 'ext' + wlSuffix, self.readFloatValue( dataStrings, self.extCol, self.totalModeIndex, wlIndex)) setattr(dataRecord, 'ssa' + wlSuffix, self.readFloatValue( dataStrings, self.ssaCol, self.totalModeIndex, wlIndex)) setattr(dataRecord, 'f11_' + wlSuffix, self.readFloatValue( dataStrings, self.f11Col, self.totalModeIndex, wlIndex)) setattr(dataRecord, 'f22by11_' + wlSuffix, self.readFloatValue( dataStrings, self.f22by11Col, self.totalModeIndex, wlIndex)) # Discars the data record if it contains invalid data. for attributeName in ( 're', 'im', 'ext', 'ssa', 'f11_', 'f22by11_'): if getattr(dataRecord, attributeName + wlSuffix) == -999.0: return None # This flag will be set to 'False' if some data required for # calculation of attenuation and backscatter coefficients are # missing (i.e., have the value of -999.0). isInputValid = True # Calculate attenuation and backscatter coefficients. # # Note: for 'Spherical' and 'Spheroid' aerosol modes, aerosol # optical thickness values specified in the AERONET data file # will correspond to total volumetric concentration of the 'Coarse' # aerosol mode (and therefore we have to specify 'vCoarse' as the # total volume for 'Spherical' and 'Spheroid' aerosol modes). This # seems to arise from the fact that AERONET only treats 'Fine' and # 'Coarse' aerosol modes as being "genuine" ones, whereas # sphericity is considered a mere tuning coefficient for these. isInputValid &= self.setAttBscattAttributes(dataRecord, dataStrings, 'Fine', dataRecord.vFine, self.fineModeIndex, wlIndex) isInputValid &= self.setAttBscattAttributes(dataRecord, dataStrings, 'Spherical', dataRecord.vCoarse, self.sphericalModeIndex, wlIndex) isInputValid &= self.setAttBscattAttributes(dataRecord, dataStrings, 'Spheroid', dataRecord.vCoarse, self.spheroidModeIndex, wlIndex) isInputValid &= self.setAttBscattAttributes(dataRecord, dataStrings, 'Coarse', dataRecord.vCoarse, self.coarseModeIndex, wlIndex) # Discars the data record if it contains invalid data. if not isInputValid: return None return dataRecord def testColumnHeader(self, headerList, requiredValue, col, modeIndex = None, wlIndex = None): """Check if the given data column header has the expected value. Raise 'txt.Error' on failure. 'headerList' should be the comlete list of header strings read from the file. 'col' should be the relative index of the data column (see 'getColumnIndex').""" columnIndex = self.getColumnIndex(col, modeIndex, wlIndex) if columnIndex >= len(headerList): raise txt.Error(self.errorPrefix + 'data column %s (%s) in missing' % (txt.quoteNumber(columnIndex + 1), txt.quote(requiredValue))) if headerList[columnIndex] != requiredValue: raise txt.Error(self.errorPrefix + 'header of column %s is %s instead of %s'% (txt.quoteNumber(columnIndex + 1), txt.quote(headerList[columnIndex]), txt.quote(requiredValue))) def readFloatValue(self, dataList, col, modeIndex = None, wlIndex = None): """Parse a floating-point value read from the data file. Raise 'txt.Error' on failure. 'dataList' should be the comlete list of data strings read from a single line of the data file. 'col' should be the relative index of the data column (see 'getColumnIndex').""" columnIndex = self.getColumnIndex(col, modeIndex, wlIndex) if columnIndex >= len(dataList): raise txt.Error(self.errorPrefix + 'missing data in column %s (line %s)' % (txt.quoteNumber(columnIndex + 1), txt.quoteNumber(self.currentLine + 1))) try: return float(dataList[columnIndex]) except ValueError: raise txt.Error(self.errorPrefix + 'invalid number in column %s (line %s)' % (txt.quoteNumber(columnIndex + 1), txt.quoteNumber(self.currentLine + 1))) def readDateTimeValue(self, dataList): """Parse the date and time values read from the AERLID the data file. Return a 'datetime.datetime' instance. Raise 'txt.Error' on failure. 'dataList' should be the comlete list of data strings read from a single line of the data file.""" dateColumn = self.getColumnIndex(self.dateCol) timeColumn = self.getColumnIndex(self.timeCol) if dateColumn >= len(dataList) or timeColumn >= len(dataList): raise txt.Error(self.errorPrefix + 'missing date or time in line %s' % txt.quoteNumber(self.currentLine + 1)) dateString = dataList[dateColumn] timeString = dataList[timeColumn] try: return datetime.datetime.strptime(dateString + timeString, '%Y-%m-%d%H:%M:%S') except ValueError: raise txt.Error(self.errorPrefix + 'invalid date or time in line %s' % txt.quoteNumber(self.currentLine + 1)) def setAttBscattAttributes(self, dataRecord, dataList, modeSuffix, volConcentration, modeIndex, wlIndex): """Set aerosol attenuation and backscatter attributes in the given 'PhotometerInput' instance 'dataRecord', for the given aerosol mode and wavelength. Return 'False' and do nothing more if any of the required data fields is invalid (i.e., has the value of -999.0). Otherwise, return 'True'. Raise 'txt.Error' on failure. Parameters: - 'dataList': the comlete list of data strings read from a single line of the data file. - 'modeSuffix': aerosol mode suffix for 'PhotometerInput' instance attributes. - 'volConcentration': total volumetric concentration of aerosol for the given aerosol mode, in micrometers, or 'None' if the value is not available for the given photometer measurement. - 'modeIndex', 'wlIndex': indices within 'self.aerosolModes' and 'self.wavelenghtsUm' arrays respectively that define aerosol mode and wavelength.""" wlSuffix = self.wavelengthSuffixes[wlIndex] extinction = self.readFloatValue(dataList, self.extCol, modeIndex, wlIndex) scattering = self.readFloatValue(dataList, self.scatteringCol, modeIndex, wlIndex) f11 = self.readFloatValue(dataList, self.f11Col, modeIndex, wlIndex) f22by11 = self.readFloatValue(dataList, self.f22by11Col, modeIndex, wlIndex) # Check if all the required data fields are valid. if any(dataField == -999.0 for dataField in (extinction, scattering, f11, f22by11)): return False # Aerosol attenuation and total aerosol backscatter coefficients. if volConcentration != 0.0: attenuation = extinction / volConcentration * 0.001 backscatter = (scattering / volConcentration * 0.001 * f11 / (4.0 * math.pi)) else: # Use arbitraty values if it's not possible to interpret the # experimental data. In this case, the values of the coefficients # won't influence the retrieval considerably. attenuation = self.getDefaultAttenuation(modeSuffix, wlSuffix) backscatter = self.getDefaultBackscatter(modeSuffix, wlSuffix) # Cross and parallel-polarized aerosol backscatter coefficients for a # linearly polarized lidar signal. backscatterC = backscatter * (1 - f22by11) * 0.5 backscatterP = backscatter * (1 + f22by11) * 0.5 setattr(dataRecord, 'a' + modeSuffix + wlSuffix, attenuation) setattr(dataRecord, 'b' + modeSuffix + wlSuffix, backscatter) setattr(dataRecord, 'b' + modeSuffix + wlSuffix + 'C', backscatterC) setattr(dataRecord, 'b' + modeSuffix + wlSuffix + 'P', backscatterP) return True @staticmethod def getDefaultAttenuation(modeSuffix, wlSuffix): """Return an arbitrary physically meaningful attenuation coefficient for the given aerosol mode and wavelength.""" # These are yearly average values for Minsk in 2010. attenuationMap = { 'Fine' : {'355' : 1.05e-2, '532': 5.92e-3, '1064' : 1.47e-3}, 'Spherical' : {'355' : 8.27e-4, '532': 8.67e-4, '1064' : 1.00e-3}, 'Spheroid' : {'355' : 9.31e-4, '532': 9.91e-4, '1064' : 1.11e-3}, 'Coarse' : {'355' : 8.55e-4, '532': 9.01e-4, '1064' : 1.03e-3}, } return attenuationMap[modeSuffix][wlSuffix] @staticmethod def getDefaultBackscatter(modeSuffix, wlSuffix): """Return an arbitrary physically meaningful backscatter coefficient for the given aerosol mode and wavelength.""" # These are yearly average values for Minsk in 2010. backscatterMap = { 'Fine' : {'355' : 1.40e-4, '532': 8.50e-5, '1064' : 3.92e-5}, 'Spherical' : {'355' : 2.40e-5, '532': 3.31e-5, '1064' : 4.09e-5}, 'Spheroid' : {'355' : 9.32e-6, '532': 1.25e-5, '1064' : 1.55e-5}, 'Coarse' : {'355' : 1.80e-5, '532': 2.43e-5, '1064' : 2.97e-5}, } return backscatterMap[modeSuffix][wlSuffix] # ***************************************************************************** class AotCorrFileLoader(AerlidFileLoader): """Utility class used for parsing of auxiliary data files that hold manually corrected values of aerosol optical thicknesses at 675 nm.""" # 'AOT corr' data files are in essence drastically reduced AERLID files, # with only date, time and aot675Corr data columns left. Thus, a modified # version of 'AerlidFileLoader' may be used to parse them. # ---- Public methods ----------------------------------------------------- def __init__(self): """Initialize auxiliary attributes used to describe the data format.""" # There are no data columns bound to aerosol modes in 'AOT corr' files. self.aerosolModes = () self.wavelengthsUm = () self.prefixHeaders = ('yyyy-mm-dd', 'hh:mm:ss', 'AOT_675-TRU') self.suffixHeaders = () self.specificHeaders = () # Indices of the required data columns. self.dateCol = self.prefixHeaders.index('yyyy-mm-dd') self.timeCol = self.prefixHeaders.index('hh:mm:ss') self.aot675CorrCol = self.prefixHeaders.index('AOT_675-TRU') def loadFile(self, filePath): """Load the specified AOT correction data file. Return a dictionary that maps 'datetime.datetime' values to floating point numbers representing corrected aerosol optical thicknesses. Raise 'txt.Error' if loading fails or file format is invalid.""" # The real loading is done in 'readFile'. Here, the file is opened, # and 'file' and 'errorPrefix' attributes are set. try: self.file = open(filePath) except IOError: raise txt.Error( 'AOT correction file %s may not be opened for reading' % txt.quotePath(filePath)) self.errorPrefix = ('AOT correction file format is invalid (%s): ' % txt.quotePath(filePath)) try: aot675CorrDict = self.readFile() except IOError: raise txt.Error('Failed to read AOT correction file (%s)' % txt.quotePath(filePath)) finally: self.file.close() return aot675CorrDict # ---- Private methods ---------------------------------------------------- def readFile(self): """Load the opened data file referenced in 'self.file'.""" # This will be a dictionary mapping 'datetime.datetime's to 'float's. aot675CorrDict = {} # Read the data column headers and check them against the expected # values. Column headers are divided by whitespace in 'AOT corr' files. dataHeaders = self.file.readline().split() for col in range(len(self.prefixHeaders)): self.testColumnHeader(dataHeaders, self.prefixHeaders[col], col) # Initialize the 'currentLine' attribute (there is only 1 header line). self.currentLine = 1 for line in self.file: # Read the data strings stored in the current line. # Data columns are divided by whitespace in 'AOT corr' files. dataStrings = line.split() dateTimeValue = self.readDateTimeValue(dataStrings) aot675CorrValue = self.readFloatValue(dataStrings, self.aot675CorrCol) # Store the loaded value in the dictionary. aot675CorrDict[dateTimeValue] = aot675CorrValue # Update the 'currentLine' attribute. self.currentLine += 1 return aot675CorrDict