diff --git a/cloudnetpy/instruments/ceilo.py b/cloudnetpy/instruments/ceilo.py index c9eb9c96..81d02542 100644 --- a/cloudnetpy/instruments/ceilo.py +++ b/cloudnetpy/instruments/ceilo.py @@ -1,8 +1,10 @@ """Module for reading and processing Vaisala / Lufft ceilometers.""" import linecache from typing import Union, Optional +import logging import numpy as np -from cloudnetpy.instruments.lufft import LufftCeilo +import netCDF4 +from cloudnetpy.instruments.lufft import LufftCeilo, Cl61d from cloudnetpy.instruments.vaisala import ClCeilo, Ct25k from cloudnetpy import utils, output, CloudnetArray from cloudnetpy.metadata import MetaData @@ -16,7 +18,7 @@ def ceilo2nc(full_path: str, date: Optional[str] = None) -> str: """Converts Vaisala / Lufft ceilometer data into Cloudnet Level 1b netCDF file. - This function reads raw Vaisala (CT25k, CL31, CL51) and Lufft (CHM15k) + This function reads raw Vaisala (CT25k, CL31, CL51, CL61-D) and Lufft (CHM15k) ceilometer files and writes the data into netCDF file. Three variants of the attenuated backscatter are saved in the file: @@ -24,6 +26,12 @@ def ceilo2nc(full_path: str, 2. Signal-to-noise screened backscatter, `beta` 3. SNR-screened backscatter with smoothed weak background, `beta_smooth` + With CL61-D three additional depolarisation parameters are saved: + + 1. Raw depolarisation, `depolarisation_raw` + 2. Signal-to-noise screened depolarisation, `depolarisation` + 3. SNR-screened depolarisation with smoothed weak background, `depolarisation_smooth` + Args: full_path: Ceilometer file name. For Vaisala it is a text file, for CHM15k it is a netCDF file. @@ -51,44 +59,58 @@ def ceilo2nc(full_path: str, >>> ceilo2nc('chm15k_raw.nc', 'chm15k.nc', site_meta) """ - ceilo = _initialize_ceilo(full_path, date) - ceilo.read_ceilometer_file(site_meta.get('calibration_factor', None)) - beta_variants = ceilo.calc_beta() - _append_data(ceilo, beta_variants) - _append_height(ceilo, site_meta['altitude']) - attributes = output.add_time_attribute(ATTRIBUTES, ceilo.date) - output.update_attributes(ceilo.data, attributes) - return _save_ceilo(ceilo, output_file, site_meta['name'], keep_uuid, uuid) + ceilo_obj = _initialize_ceilo(full_path, date) + logging.debug('reading daily file') + ceilo_obj.read_ceilometer_file(site_meta.get('calibration_factor', None)) + if 'cl61' in ceilo_obj.model.lower(): + depol_variants = ceilo_obj.calc_depol() + else: + depol_variants = None + beta_variants = ceilo_obj.calc_beta() + _append_data(ceilo_obj, beta_variants, depol_variants) + _append_height(ceilo_obj, site_meta['altitude']) + attributes = output.add_time_attribute(ATTRIBUTES, ceilo_obj.date) + output.update_attributes(ceilo_obj.data, attributes) + return _save_ceilo(ceilo_obj, output_file, site_meta['name'], keep_uuid, uuid) def _initialize_ceilo(full_path: str, - date: Union[str, None]) -> Union[ClCeilo, Ct25k, LufftCeilo]: + date: Optional[str] = None) -> Union[ClCeilo, Ct25k, LufftCeilo, Cl61d]: model = _find_ceilo_model(full_path) if model == 'cl31_or_cl51': return ClCeilo(full_path, date) if model == 'ct25k': return Ct25k(full_path) + if model == 'cl61d': + return Cl61d(full_path, date) return LufftCeilo(full_path, date) def _find_ceilo_model(full_path: str) -> str: - if full_path.lower().endswith('.nc'): + try: + nc = netCDF4.Dataset(full_path) + title = nc.title + nc.close() + for identifier in ['cl61d', 'cl61-d']: + if identifier in title.lower() or identifier in full_path.lower(): + return 'cl61d' return 'chm15k' - first_empty_line = utils.find_first_empty_line(full_path) - max_number_of_empty_lines = 10 - for n in range(1, max_number_of_empty_lines): - line = linecache.getline(full_path, first_empty_line + n) - if not utils.is_empty_line(line): - line = linecache.getline(full_path, first_empty_line + n + 1) - break - if 'CL' in line: - return 'cl31_or_cl51' - if 'CT' in line: - return 'ct25k' + except OSError: + first_empty_line = utils.find_first_empty_line(full_path) + max_number_of_empty_lines = 10 + for n in range(1, max_number_of_empty_lines): + line = linecache.getline(full_path, first_empty_line + n) + if not utils.is_empty_line(line): + line = linecache.getline(full_path, first_empty_line + n + 1) + break + if 'CL' in line: + return 'cl31_or_cl51' + if 'CT' in line: + return 'ct25k' raise RuntimeError('Error: Unknown ceilo model.') -def _append_height(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], +def _append_height(ceilo: Union[ClCeilo, Ct25k, LufftCeilo, Cl61d], site_altitude: float) -> None: """Finds height above mean sea level.""" tilt_angle = np.median(ceilo.metadata['tilt_angle']) @@ -98,11 +120,16 @@ def _append_height(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], ceilo.data['altitude'] = CloudnetArray(site_altitude, 'altitude') -def _append_data(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], - beta_variants: tuple): +def _append_data(ceilo: Union[ClCeilo, Ct25k, LufftCeilo, Cl61d], + beta_variants: tuple, + depol_variants: Optional[tuple] = None): """Adds data / metadata as CloudnetArrays to ceilo.data.""" for data, name in zip(beta_variants, ('beta_raw', 'beta', 'beta_smooth')): ceilo.data[name] = CloudnetArray(data, name) + if depol_variants is not None: + for data, name in zip(depol_variants, ('depolarisation', 'depolarisation_smooth')): + ceilo.data[name] = CloudnetArray(data, name) + del ceilo.data['beta_raw'] for field in ('range', 'time', 'wavelength', 'calibration_factor'): ceilo.data[field] = CloudnetArray(np.array(getattr(ceilo, field)), field) for field, data in ceilo.metadata.items(): @@ -111,16 +138,14 @@ def _append_data(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], ceilo.data[field] = CloudnetArray(np.array(ceilo.metadata[field], dtype=float), field) -def _save_ceilo(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], +def _save_ceilo(ceilo: Union[ClCeilo, Ct25k, LufftCeilo, Cl61d], output_file: str, location: str, keep_uuid: bool, uuid: Union[str, None]) -> str: """Saves the ceilometer netcdf-file.""" - dims = {'time': len(ceilo.time), 'range': len(ceilo.range)} - rootgrp = output.init_file(output_file, dims, ceilo.data, keep_uuid, uuid) uuid = rootgrp.file_uuid output.add_file_type(rootgrp, 'lidar') @@ -151,6 +176,31 @@ def _save_ceilo(ceilo: Union[ClCeilo, Ct25k, LufftCeilo], comment=('Range corrected, SNR screened backscatter coefficient.\n' 'Weak background is smoothed using Gaussian 2D-kernel.') ), + 'depolarisation': MetaData( + long_name='Depolarisation', + units='%', + comment='SNR screened lidar depolarisation' + ), + 'depolarisation_raw': MetaData( + long_name='Raw depolarisation', + units='%', + ), + 'depolarisation_smooth': MetaData( + long_name='Smoothed lidar depolarisation', + units='%', + comment=('SNR screened lidar depolarisation.\n' + 'Weak background is smoothed using Gaussian 2D-kernel.') + ), + 'p_pol': MetaData( + long_name='Raw attenuated backscatter coefficient - parallel component', + units='sr-1 m-1', + comment='Range-corrected, parallel-polarized component of attenuated backscatter.' + ), + 'x_pol': MetaData( + long_name='Raw attenuated backscatter coefficient - cross component', + units='sr-1 m-1', + comment='Range-corrected, cross-polarized component of attenuated backscatter.' + ), 'scale': MetaData( long_name='Scale', units='%', diff --git a/cloudnetpy/instruments/ceilometer.py b/cloudnetpy/instruments/ceilometer.py index 7c82f7ff..040e238e 100644 --- a/cloudnetpy/instruments/ceilometer.py +++ b/cloudnetpy/instruments/ceilometer.py @@ -10,81 +10,134 @@ class Ceilometer: def __init__(self, full_path: str): self.file_name = full_path - self.backscatter = np.array([]) + self.model = '' + self.processed_data = {} self.data = {} self.metadata = {} self.range = np.array([]) + self.range_squared = np.array([]) self.time = [] self.date = [] self.noise_params = (1, 1, 1, (1, 1)) - self.calibration_factor = None - self.model = None + self.calibration_factor = 1 self.wavelength = None def calc_beta(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Converts range-corrected raw beta to noise-screened beta.""" - - def _screen_beta(beta_in: np.ndarray, smooth: bool) -> np.ndarray: - beta_in = _calc_range_uncorrected_beta(beta_in, range_squared) - beta_in = self._screen_by_snr(beta_in, is_saturation, beta_is_smoothed=smooth) - return _calc_range_corrected_beta(beta_in, range_squared) - - range_squared = _get_range_squared(self.range) - is_saturation = self._find_saturated_profiles() - beta = _screen_beta(self.backscatter, False) - # smoothed version: - beta_smooth = ma.copy(self.backscatter) + snr_limit = 5 + noisy_data = NoisyData(*self._get_args(), snr_limit) + beta = noisy_data.screen_data(self.processed_data['backscatter']) + beta_smooth = self._calc_beta_smooth(beta) + beta_smooth = noisy_data.screen_data(beta_smooth, is_smoothed=True) + return self.processed_data['backscatter'], beta, beta_smooth + + def calc_depol(self) -> Tuple[np.ndarray, np.ndarray]: + """Converts raw depolarisation to noise-screened depolarisation.""" + snr_limit = 4 + noisy_data = NoisyData(*self._get_args(), snr_limit) + sigma = _calc_sigma_units(self.time, self.range) + x_pol = noisy_data.screen_data(self.processed_data['x_pol'], keep_negative=True) + depol = x_pol / self.processed_data['p_pol'] + p_pol_smooth = scipy.ndimage.filters.gaussian_filter(self.processed_data['p_pol'], sigma) + x_pol_smooth = scipy.ndimage.filters.gaussian_filter(self.processed_data['x_pol'], sigma) + x_pol_smooth = noisy_data.screen_data(x_pol_smooth, is_smoothed=True) + depol_smooth = x_pol_smooth / p_pol_smooth + return depol, depol_smooth + + def _calc_beta_smooth(self, beta: np.ndarray) -> np.ndarray: + beta_smooth = ma.copy(self.processed_data['backscatter']) cloud_ind, cloud_values, cloud_limit = _estimate_clouds_from_beta(beta) beta_smooth[cloud_ind] = cloud_limit sigma = _calc_sigma_units(self.time, self.range) beta_smooth = scipy.ndimage.filters.gaussian_filter(beta_smooth, sigma) beta_smooth[cloud_ind] = cloud_values - beta_smooth = _screen_beta(beta_smooth, True) - return self.backscatter, beta, beta_smooth + return beta_smooth + + def _get_range_squared(self) -> np.ndarray: + """Returns range (m), squared and converted to km.""" + m2km = 0.001 + return (self.range * m2km) ** 2 + + def _get_args(self): + return self.processed_data, self.range_squared, self.noise_params + + +class NoisyData: + def __init__(self, + data: dict, + range_squared: np.ndarray, + noise_params: tuple, + snr_limit: float): + self.data = data + self.range_squared = range_squared + self.noise_params = noise_params + self.snr_limit = snr_limit + self._is_saturation = self._find_saturated_profiles() + + def screen_data(self, + array: np.ndarray, + is_smoothed: Optional[bool] = False, + keep_negative: Optional[bool] = False) -> np.ndarray: + array = self._calc_range_uncorrected(array) + array = self._screen_by_snr(array, is_smoothed, keep_negative) + array = self._calc_range_corrected(array) + return array def _screen_by_snr(self, - beta_uncorrected: np.ndarray, - is_saturation: np.ndarray, - beta_is_smoothed: Optional[bool] = False) -> np.ndarray: - """Screens noise from ceilometer backscatter. - - Args: - beta_uncorrected: Range-uncorrected backscatter. - is_saturation: Boolean array denoting saturated profiles. - beta_is_smoothed: Should be true if input beta is smoothed. Default is False. - - """ - beta = ma.copy(beta_uncorrected) + array: np.ndarray, + is_smoothed: bool, + keep_negative: bool) -> np.ndarray: + """Screens noise from range-uncorrected lidar variable.""" n_gates, _, saturation_noise, noise_min = self.noise_params - noise_min = noise_min[0] if beta_is_smoothed else noise_min[1] - noise = _estimate_noise_from_top_gates(beta, n_gates, noise_min) - beta = _reset_low_values_above_saturation(beta, is_saturation, saturation_noise) - beta = _remove_noise(beta, noise) - return beta + noise_min = noise_min[0] if is_smoothed is True else noise_min[1] + noise = _estimate_noise_from_top_gates(array, n_gates, noise_min) + array = self._reset_low_values_above_saturation(array, saturation_noise) + array = self._remove_noise(array, noise, keep_negative) + return array def _find_saturated_profiles(self) -> np.ndarray: """Estimates saturated profiles using the variance of the top range gates.""" n_gates, var_lim, _, _ = self.noise_params - var = np.var(self.backscatter[:, -n_gates:], axis=1) + var = np.var(self.data['backscatter'][:, -n_gates:], axis=1) return var < var_lim - -def _remove_noise(beta_in: np.ndarray, noise: np.ndarray) -> np.ndarray: - beta = ma.copy(beta_in) - snr_limit = 5 - snr = (beta.T / noise) - beta[snr.T < snr_limit] = ma.masked - return beta - - -def _calc_sigma_units(time: list, range_los: np.ndarray) -> Tuple[float, float]: + def _reset_low_values_above_saturation(self, + array: np.ndarray, + saturation_noise: float) -> np.ndarray: + """Removes low values in saturated profiles above peak.""" + for saturated_profile in np.where(self._is_saturation)[0]: + profile = array[saturated_profile, :] + peak_ind = np.argmax(profile) + alt_ind = np.where(profile[peak_ind:] < saturation_noise)[0] + peak_ind + array[saturated_profile, alt_ind] = ma.masked + return array + + def _remove_noise(self, + array: np.ndarray, + noise: np.ndarray, + keep_negative: bool) -> np.ndarray: + snr = array / utils.transpose(noise) + if keep_negative is True: + array[np.abs(snr) < self.snr_limit] = ma.masked + else: + array[snr < self.snr_limit] = ma.masked + return array + + def _calc_range_uncorrected(self, array: np.ndarray) -> np.ndarray: + return array / self.range_squared + + def _calc_range_corrected(self, array: np.ndarray) -> np.ndarray: + return array * self.range_squared + + +def _calc_sigma_units(time_vector: list, range_los: np.ndarray) -> Tuple[float, float]: """Calculates Gaussian peak std parameters. The amount of smoothing is hard coded. This function calculates how many steps in time and height corresponds to this smoothing. Args: - time: 1D vector (fraction hour). + time_vector: 1D vector (fraction hour). range_los: 1D vector (m). Returns: @@ -95,7 +148,7 @@ def _calc_sigma_units(time: list, range_los: np.ndarray) -> Tuple[float, float]: minutes_in_hour = 60 sigma_minutes = 2 sigma_metres = 5 - time_step = utils.mdiff(time) * minutes_in_hour + time_step = utils.mdiff(time_vector) * minutes_in_hour alt_step = utils.mdiff(range_los) x_std = sigma_minutes / time_step y_std = sigma_metres / alt_step @@ -109,55 +162,8 @@ def _estimate_clouds_from_beta(beta: np.ndarray) -> Tuple[np.ndarray, np.ndarray return cloud_ind, beta[cloud_ind], cloud_limit -def _estimate_noise_from_top_gates(beta: np.ndarray, n_gates: int, noise_min: float) -> np.ndarray: - """Estimates backscatter noise from topmost range gates.""" - noise = ma.std(beta[:, -n_gates:], axis=1) +def _estimate_noise_from_top_gates(data: np.ndarray, n_gates: int, noise_min: float) -> np.ndarray: + """Estimates noise from topmost range gates.""" + noise = ma.std(data[:, -n_gates:], axis=1) noise[noise < noise_min] = noise_min return noise - - -def _reset_low_values_above_saturation(beta_in: np.ndarray, - is_saturation: np.ndarray, - saturation_noise: float) -> np.ndarray: - """Removes low values in saturated profiles above peak.""" - beta = ma.copy(beta_in) - for saturated_profile in np.where(is_saturation)[0]: - profile = beta[saturated_profile, :] - peak_ind = np.argmax(profile) - alt_ind = np.where(profile[peak_ind:] < saturation_noise)[0] + peak_ind - beta[saturated_profile, alt_ind] = ma.masked - return beta - - -def _calc_range_uncorrected_beta(beta: np.ndarray, range_squared: np.ndarray) -> np.ndarray: - """Calculates range uncorrected beta. - - Args: - beta: 2D attenuated backscatter. - range_squared: 1D altitude vector (km), squared. - - Returns: - ndarray: 2D range uncorrected beta. - - """ - return beta / range_squared - - -def _calc_range_corrected_beta(beta: np.ndarray, range_squared: np.ndarray) -> np.ndarray: - """Calculates range corrected beta. - - Args: - beta: 2D attenuated backscatter. - range_squared: 1D altitude vector (km), squared. - - Returns: - ndarray: 2D range corrected beta. - - """ - return beta * range_squared - - -def _get_range_squared(range_instru: np.ndarray) -> np.ndarray: - """Returns range (m), squared and converted to km.""" - m2km = 0.001 - return (range_instru*m2km)**2 diff --git a/cloudnetpy/instruments/lufft.py b/cloudnetpy/instruments/lufft.py index c2dfdca4..3524ae31 100644 --- a/cloudnetpy/instruments/lufft.py +++ b/cloudnetpy/instruments/lufft.py @@ -14,14 +14,14 @@ def __init__(self, file_name: str, date: Optional[str] = None): self._expected_date = date self.model = 'Lufft CHM15k' self.dataset = netCDF4.Dataset(self.file_name) - self.variables = self.dataset.variables self.noise_params = (70, 2e-14, 0.3e-6, (1e-9, 4e-9)) self.wavelength = 1064 def read_ceilometer_file(self, calibration_factor: Optional[float] = None) -> None: """Reads data and metadata from Jenoptik netCDF file.""" self.range = self._calc_range() - self.backscatter = self._calibrate_backscatter(calibration_factor) + self.range_squared = self._get_range_squared() + self.processed_data['backscatter'] = self._calibrate_backscatter(calibration_factor) self.time = self._fetch_time() self.date = self._read_date() self.metadata = self._read_metadata() @@ -31,33 +31,35 @@ def _calc_range(self) -> np.ndarray: ceilo_range = self._getvar('range') return ceilo_range - utils.mdiff(ceilo_range)/2 - def _calibrate_backscatter(self, calibration_factor: Union[float, None]) -> np.ndarray: + def _calibrate_backscatter(self, calibration_factor: Optional[float] = None) -> np.ndarray: beta_raw = self._getvar('beta_raw') overlap_function = _get_overlap(self.range) beta_raw /= overlap_function if calibration_factor is None: - logging.warning('Using default calibration factor for CHM15k') + logging.warning('Using default calibration factor') calibration_factor = 3e-12 self.calibration_factor = calibration_factor beta_raw *= calibration_factor return beta_raw def _fetch_time(self) -> np.ndarray: - time = self.variables['time'][:] + time = self.dataset.variables['time'][:] ind = time.argsort() time = time[ind] - self.backscatter = self.backscatter[ind, :] + for key in self.processed_data: + self.processed_data[key] = self.processed_data[key][ind, :] if self._expected_date is not None: - epoch = utils.get_epoch(self.variables['time'].units) + epoch = utils.get_epoch(self.dataset.variables['time'].units) valid_ind = [] for ind, timestamp in enumerate(time): date = '-'.join(utils.seconds2date(timestamp, epoch)[:3]) if date == self._expected_date: valid_ind.append(ind) if not valid_ind: - raise ValueError('Error: CHM15k date differs from expected.') + raise ValueError(f'Error: {self.model} date differs from expected.') time = time[valid_ind] - self.backscatter = self.backscatter[valid_ind, :] + for key in self.processed_data: + self.processed_data[key] = self.processed_data[key][valid_ind, :] return utils.seconds2hours(time) def _read_date(self) -> List[str]: @@ -67,13 +69,63 @@ def _read_date(self) -> List[str]: def _getvar(self, *args) -> Union[np.ndarray, float, None]: for arg in args: - if arg in self.variables: - var = self.variables[arg] + if arg in self.dataset.variables: + var = self.dataset.variables[arg] return var[0] if utils.isscalar(var) else var[:] return None def _read_metadata(self) -> dict: - return {'tilt_angle': self._getvar('zenith')} + tilt_angle = self._getvar('zenith') # 0 deg is vertical + if tilt_angle is None: + tilt_angle = 0 + logging.warning(f'Assuming {tilt_angle} deg tilt angle') + return {'tilt_angle': tilt_angle} + + +class Cl61d(LufftCeilo): + """Class for Vaisala CL61d ceilometer.""" + def __init__(self, file_name: str, date: Optional[str] = None): + super().__init__(file_name) + self._expected_date = date + self.model = 'Vaisala CL61d' + self.wavelength = 910.55 + self.noise_params = (100, 2e-14, 0.3e-6, (1e-9, 4e-9)) + + def read_ceilometer_file(self, calibration_factor: Optional[float] = None) -> None: + """Reads data and metadata from concatenated Vaisala CL61d netCDF file.""" + self.range = self._calc_range() + self.range_squared = self._get_range_squared() + self.processed_data['backscatter'] = self._calibrate_backscatter(calibration_factor) + for key in ('p_pol', 'x_pol', 'linear_depol_ratio'): + self.processed_data[key] = self._getvar(key) + self.time = self._fetch_time() + self.date = self._read_date() + self.metadata = self._read_metadata() + + def _calc_range(self) -> np.ndarray: + """Assumes 'range' means the lower limit of range gate.""" + ceilo_range = self._getvar('range') + return ceilo_range + utils.mdiff(ceilo_range)/2 + + def _read_date(self) -> List[str]: + if self._expected_date: + return self._expected_date.split('-') + time = self.dataset.variables['time'][:] + date_first = utils.seconds2date(time[0], epoch=(1970, 1, 1)) + date_last = utils.seconds2date(time[-1], epoch=(1970, 1, 1)) + date_middle = utils.seconds2date(time[round(len(time)/2)], epoch=(1970, 1, 1)) + if date_first != date_last: + logging.warning('No expected date given and different dates in CL61d timestamps.') + return date_middle[:3] + + def _calibrate_backscatter(self, calibration_factor: Optional[float] = None) -> np.ndarray: + beta_raw = self._getvar('beta_att') + if calibration_factor is None: + logging.warning('Using default calibration factor') + calibration_factor = 1 + self.calibration_factor = calibration_factor + beta_raw *= calibration_factor + return beta_raw def _get_overlap(range_ceilo: np.ndarray, diff --git a/cloudnetpy/instruments/vaisala.py b/cloudnetpy/instruments/vaisala.py index 0f640216..d9899f3c 100644 --- a/cloudnetpy/instruments/vaisala.py +++ b/cloudnetpy/instruments/vaisala.py @@ -167,7 +167,7 @@ def _read_header_line_2(lines: list) -> dict: def _range_correct_upper_part(self) -> None: altitude_limit = 2400 ind = np.where(self.range < altitude_limit) - self.backscatter[:, ind] *= (self.range[ind]*M2KM)**2 + self.processed_data['backscatter'][:, ind] *= (self.range[ind]*M2KM)**2 class ClCeilo(VaisalaCeilo): @@ -185,13 +185,14 @@ def read_ceilometer_file(self, calibration_factor: Optional[float] = None) -> No header.append(self._read_header_line_4(data_lines[-3])) self.metadata = self._handle_metadata(header) self.range = self._calc_range() - self.backscatter = self._read_backscatter(data_lines[-2]) + self.range_squared = self._get_range_squared() + self.processed_data['backscatter'] = self._read_backscatter(data_lines[-2]) self.calibration_factor = calibration_factor or 1 - self.backscatter *= self.calibration_factor + self.processed_data['backscatter'] *= self.calibration_factor self._store_ceilometer_info() def _store_ceilometer_info(self): - n_gates = self.backscatter.shape[1] + n_gates = self.processed_data['backscatter'].shape[1] if n_gates < 1000: self.model = 'cl31' self.wavelength = '910' @@ -236,10 +237,11 @@ def read_ceilometer_file(self, calibration_factor: Optional[float] = None) -> No header.append(self._read_header_line_3(data_lines[3])) self.metadata = self._handle_metadata(header) self.range = self._calc_range() + self.range_squared = self._get_range_squared() hex_profiles = self._parse_hex_profiles(data_lines[4:20]) - self.backscatter = self._read_backscatter(hex_profiles) + self.processed_data['backscatter'] = self._read_backscatter(hex_profiles) self.calibration_factor = calibration_factor or 1 - self.backscatter *= self.calibration_factor + self.processed_data['backscatter'] *= self.calibration_factor # TODO: should study the background noise to determine if the # next call is needed. It can be the case with cl31/51 also. # self._range_correct_upper_part() diff --git a/cloudnetpy/plotting/plot_meta.py b/cloudnetpy/plotting/plot_meta.py index 7e3ed754..5ef6a614 100644 --- a/cloudnetpy/plotting/plot_meta.py +++ b/cloudnetpy/plotting/plot_meta.py @@ -318,6 +318,30 @@ plot_scale=_LOG, plot_type='mesh' ), + 'depolarisation_raw': PlotMeta( + name='Raw depolarisation', + cbar='viridis', + clabel='', + plot_range=(0, 1), + plot_scale=_LIN, + plot_type='mesh' + ), + 'depolarisation': PlotMeta( + name='Lidar depolarisation', + cbar='viridis', + clabel='', + plot_range=(0, 1), + plot_scale=_LIN, + plot_type='mesh' + ), + 'depolarisation_smooth': PlotMeta( + name='Lidar depolarisation (smoothed)', + cbar='viridis', + clabel='', + plot_range=(0, 1), + plot_scale=_LIN, + plot_type='mesh' + ), 'Z': PlotMeta( name='Radar reflectivity factor', cbar='viridis', diff --git a/cloudnetpy/quality/data_quality_config.ini b/cloudnetpy/quality/data_quality_config.ini index 598e176b..cbe1461d 100644 --- a/cloudnetpy/quality/data_quality_config.ini +++ b/cloudnetpy/quality/data_quality_config.ini @@ -5,8 +5,8 @@ v = -10, 10 width = 0, 10 ldr = -40, 20 beta_raw = -1e-4, 1 -beta = 1e-13, 1 -beta_smooth = 1e-13, 1 +beta = 0, 1 +beta_smooth = 0, 1 v_sigma = 0, 5 insect_prob = 0, 1 is_rain = 0, 1 diff --git a/cloudnetpy/quality/metadata_config.ini b/cloudnetpy/quality/metadata_config.ini index 330b212a..9284c2cc 100644 --- a/cloudnetpy/quality/metadata_config.ini +++ b/cloudnetpy/quality/metadata_config.ini @@ -29,8 +29,7 @@ drizzle = Do, mu, S, beta_corr, drizzle_N, drizzle_lwc, drizzle_lwf, v_drizzle, v_drizzle_bias, drizzle_retrieval_status, height, time, altitude, latitude, longitude [required_global_attributes] -all = year, month, day, file_uuid, cloudnetpy_version, Conventions, location, history, title, - references +all = year, month, day, file_uuid, cloudnetpy_version, Conventions, location, history, title [attribute_limits] month = 1, 12 diff --git a/tests/unit/data/cl61d/live_20210829_000020.nc b/tests/unit/data/cl61d/live_20210829_000020.nc new file mode 100644 index 00000000..0af72ac6 Binary files /dev/null and b/tests/unit/data/cl61d/live_20210829_000020.nc differ diff --git a/tests/unit/data/cl61d/live_20210829_104420.nc b/tests/unit/data/cl61d/live_20210829_104420.nc new file mode 100644 index 00000000..a21055c7 Binary files /dev/null and b/tests/unit/data/cl61d/live_20210829_104420.nc differ diff --git a/tests/unit/data/cl61d/live_20210829_224520.nc b/tests/unit/data/cl61d/live_20210829_224520.nc new file mode 100644 index 00000000..338b91a8 Binary files /dev/null and b/tests/unit/data/cl61d/live_20210829_224520.nc differ diff --git a/tests/unit/data/cl61d/live_20210830_035020.nc b/tests/unit/data/cl61d/live_20210830_035020.nc new file mode 100644 index 00000000..3847cf13 Binary files /dev/null and b/tests/unit/data/cl61d/live_20210830_035020.nc differ diff --git a/tests/unit/test_ceilo.py b/tests/unit/test_ceilo.py index 36a54cb2..03fe8ec2 100644 --- a/tests/unit/test_ceilo.py +++ b/tests/unit/test_ceilo.py @@ -11,7 +11,13 @@ def test_find_ceilo_model_jenoptik(): - assert ceilo._find_ceilo_model('ceilo.nc') == 'chm15k' + file = f'{SCRIPT_PATH}/data/chm15k/00100_A202010220005_CHM170137.nc' + assert ceilo._find_ceilo_model(file) == 'chm15k' + + +def test_find_ceilo_model_cl61d(): + file = f'{SCRIPT_PATH}/data/cl61d/live_20210829_224520.nc' + assert ceilo._find_ceilo_model(file) == 'cl61d' @pytest.mark.parametrize("fix, result", [ diff --git a/tests/unit/test_ceilometer.py b/tests/unit/test_ceilometer.py index a050c204..edecbf43 100644 --- a/tests/unit/test_ceilometer.py +++ b/tests/unit/test_ceilometer.py @@ -2,49 +2,81 @@ import numpy.ma as ma from numpy.testing import assert_array_equal, assert_array_almost_equal from cloudnetpy.instruments import ceilometer - - -def test_remove_noise(): - noise = 2 - beta = ma.array([[1, 2, 3], - [1, 2, 3]], mask=False) - result = ma.array([[1, 2, 3], - [1, 2, 3]], - mask=[[1, 0, 0], - [1, 0, 0]]) - assert_array_equal(ceilometer._remove_noise(beta.mask, noise), result.mask) - - -def test_calc_sigma_units(): - time = np.linspace(0, 24, 721) # 2 min resolution - range_instru = np.arange(0, 1000, 5) # 5 m resolution - std_time, std_range = ceilometer._calc_sigma_units(time, range_instru) - assert_array_almost_equal(std_time, 1) - assert_array_almost_equal(std_range, 1) - - -def test_get_range_squared(): - range_instru = np.array([1000, 2000, 3000]) - result = np.array([1, 4, 9]) - assert_array_equal(ceilometer._get_range_squared(range_instru), result) - - -def test_calc_range_uncorrected_beta(): - beta = np.array([[1, 2, 3], - [1, 2, 3]]) - range_squared = np.array([1, 2, 3]) - result = np.ones((2, 3)) - assert_array_equal(ceilometer._calc_range_uncorrected_beta(beta, range_squared), - result) - - -def test_calc_range_corrected_beta(): - result = np.array([[1, 2, 3], - [1, 2, 3]]) - range_squared = np.array([1, 2, 3]) - beta = np.ones((2, 3)) - assert_array_equal(ceilometer._calc_range_corrected_beta(beta, range_squared), - result) +import pytest + + +class TestNoisyData: + + @pytest.fixture(autouse=True) + def run_before_and_after_tests(self): + data = {'backscatter': np.array([[1, 2, 3], + [1, 2, 3]])} + range_squared = np.array([1, 2, 3]) + snr_limit = 1 + noise_params = (1, 1, 1, (1, 1)) + self.noisy_data = ceilometer.NoisyData(data, range_squared, noise_params, snr_limit) + yield + + def test_remove_noise(self): + noise = np.array([1.1, -1.1]) + array = ma.array([[1, 2, 3], + [1, 2, 3]], mask=False) + expected = ma.array([[1, 2, 3], + [1, 2, 3]], mask=[[1, 0, 0], + [1, 0, 0]]) + screened_array = self.noisy_data._remove_noise(array, noise, True) + assert_array_equal(screened_array.mask, expected.mask) + + def test_calc_range_uncorrected(self): + beta = np.array([[1, 2, 3], + [1, 2, 3]]) + expected = np.ones((2, 3)) + assert_array_equal(self.noisy_data._calc_range_uncorrected(beta), expected) + + def test_calc_range_corrected(self): + beta = np.ones((2, 3)) + expected = np.array([[1, 2, 3], + [1, 2, 3]]) + assert_array_equal(self.noisy_data._calc_range_corrected(beta), expected) + + def test_reset_low_values_above_saturation(self): + beta = ma.array([[0, 10, 1e-6, 3], + [0, 0, 0, 0.1], + [0, 0.6, 1.2, 1e-8]]) + noise = 1e-3 + self.noisy_data._is_saturation = np.array([1, 0, 1]) + expected = ma.array([[0, 10, 1e-6, 3], + [0, 0, 0, 0.1], + [0, 0.6, 1.2, 1e-8]], mask=[[0, 0, 1, 0], + [0, 0, 0, 0], + [0, 0, 0, 1]]) + result = self.noisy_data._reset_low_values_above_saturation(beta, noise) + assert_array_equal(result.data, expected.data) + assert_array_equal(result.mask, expected.mask) + + def test_find_saturated_profiles(self): + self.noisy_data.noise_params = (2, 0.25, 1, (1, 1)) + self.noisy_data.data['backscatter'] = np.array([[0, 10, 1, 1.99], + [0, 10, 2.1, 1], + [0, 10, 1, 1]]) + result = [1, 0, 1] + assert_array_equal(self.noisy_data._find_saturated_profiles(), result) + + +class TestCeilometer: + + def test_calc_sigma_units(self): + time = np.linspace(0, 24, 721) # 2 min resolution + range_instru = np.arange(0, 1000, 5) # 5 m resolution + std_time, std_range = ceilometer._calc_sigma_units(time, range_instru) + assert_array_almost_equal(std_time, 1) + assert_array_almost_equal(std_range, 1) + + def test_get_range_squared(self): + obj = ceilometer.Ceilometer('/foo/bar') + obj.range = np.array([1000, 2000, 3000]) + result = np.array([1, 4, 9]) + assert_array_equal(obj._get_range_squared(), result) def test_estimate_clouds_from_beta(): @@ -74,30 +106,3 @@ def test_estimate_noise_from_top_gates(): result = np.array([np.std([0.4, 0.1, 0.2]), np.std([0, 0, 0.1]), noise_min]) assert_array_equal(ceilometer._estimate_noise_from_top_gates(beta, 3, noise_min), result) - - -def test_reset_low_values_above_saturation(): - beta = ma.array([[0, 10, 1e-6, 3], - [0, 0, 0, 0.1], - [0, 0.6, 1.2, 1e-8]]) - noise = 1e-3 - saturated = [1, 0, 1] - result = ma.array([[0, 10, 1e-6, 3], - [0, 0, 0, 0.1], - [0, 0.6, 1.2, 1e-8]], - mask=[[0, 0, 1, 0], - [0, 0, 0, 0], - [0, 0, 0, 1]]) - call = ceilometer._reset_low_values_above_saturation(beta, saturated, noise) - assert_array_equal(call.data, result.data) - assert_array_equal(call.mask, result.mask) - - -def test_find_saturated_profiles(): - obj = ceilometer.Ceilometer('ceilo.txt') - obj.noise_params = (2, 0.25, 1, (1, 1)) - obj.backscatter = np.array([[0, 10, 1, 1.99], - [0, 10, 2.1, 1], - [0, 10, 1, 1]]) - result = [1, 0, 1] - assert_array_equal(obj._find_saturated_profiles(), result) diff --git a/tests/unit/test_cl61d.py b/tests/unit/test_cl61d.py new file mode 100644 index 00000000..f34c5629 --- /dev/null +++ b/tests/unit/test_cl61d.py @@ -0,0 +1,58 @@ +import os +import glob +from cloudnetpy import concat_lib +from cloudnetpy.instruments import ceilo2nc +import pytest +import netCDF4 + +SCRIPT_PATH = os.path.dirname(os.path.realpath(__file__)) + + +class TestCl61d: + + site_meta = { + 'name': 'Hyytiälä', + 'altitude': 123, + 'calibration_factor': 2.0 + } + files = glob.glob(f'{SCRIPT_PATH}/data/cl61d/*.nc') + files.sort() + + @pytest.fixture(autouse=True) + def run_before_and_after_tests(self): + self.output = 'dummy_output_file.nc' + self.filename = 'dummy_daily_file.nc' + concat_lib.concatenate_files(self.files, self.filename, concat_dimension='profile') + yield + os.remove(self.filename) + os.remove(self.output) + + def test_variables(self): + ceilo2nc(self.filename, self.output, self.site_meta) + nc = netCDF4.Dataset(self.output) + for key in ('beta', 'depolarisation', 'beta_smooth', 'depolarisation_smooth'): + assert key in nc.variables + for key in ('beta_raw', 'depolarisation_raw'): + assert key not in nc.variables + for key in ('altitude', 'calibration_factor'): + assert nc.variables[key][:] == self.site_meta[key] + nc.close() + + def test_global_attributes(self): + uuid = ceilo2nc(self.filename, self.output, self.site_meta) + nc = netCDF4.Dataset(self.output) + assert nc.source == 'Vaisala CL61d' + assert nc.location == self.site_meta['name'] + assert nc.title == f'Ceilometer file from {self.site_meta["name"]}' + assert nc.file_uuid == uuid + assert nc.cloudnet_file_type == 'lidar' + nc.close() + + def test_date_argument(self): + ceilo2nc(self.filename, self.output, self.site_meta, date='2021-08-30') + nc = netCDF4.Dataset(self.output) + assert len(nc.variables['time']) == 12 + assert nc.year == '2021' + assert nc.month == '08' + assert nc.day == '30' + nc.close()