Skip to content

Data Classes

ProfilesData

The ProfilesData class contains profile data information. Most of the information can be found in the data attribute, which is an xarray.Dataset. Detection and retrieval methods might add information as additional xarray.DataArray.

ProfilesData

Base class representing profiles data returned by (aprofiles.reader.ReadProfiles):.

Source code in aprofiles/profiles.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
class ProfilesData:
    """
    Base class representing profiles data returned by (aprofiles.reader.ReadProfiles):.
    """

    def __init__(self, data):
        self.data = data

    @property
    def data(self):
        """
        Data attribute (instance of (xarray.Dataset):)
        """
        return self._data

    @data.setter
    def data(self, data):
        if not isinstance(data, xr.Dataset):
            raise ValueError("Wrong data type: an xarray Dataset is expected.")
        self._data = data

    def _get_index_from_altitude_AGL(self, altitude):
        """
        Returns the closest index of the ProfilesData vertical dimension to a given AGL altitude

        Args:
            altitude (float): in m, altitude AGL to look for

        Returns:
            (int): Closest index of the vertical dimension to the given altitude AGL
        """
        altitude_asl = altitude + self.data.station_altitude.data
        return int(np.argmin(abs(self.data.altitude.data - altitude_asl)))

    def _get_resolution(self, which):
        """
        Returns the resolution of a given dimension. Supports 'altitude' and 'time'.
        The altitude resolution is given in meters, while the time resolution is given in seconds.

        Args:
            which (str): Dimension: must be ['altitude', 'time'].

        Returns:
            (float): resolution, in m (if which=='altitude') or in s (if which=='time')
        """
        if which == "altitude":
            return min(np.diff(self.data.altitude.data))
        elif which == "time":
            return (
                min(np.diff(self.data.time.data)).astype("timedelta64[s]").astype(int)
            )

    def _get_lowest_clouds(self):
        # returns an array of the altitude (in m, ASL) of the lowest cloud for each timestamp
        lowest_clouds = []
        for i in np.arange(len(self.data.time.data)):
            i_clouds = apro.utils.get_true_indexes(self.data.clouds_bases.data[i, :])
            if len(i_clouds) > 0:
                lowest_clouds.append(self.data.altitude.data[i_clouds[0]])
            else:
                lowest_clouds.append(np.nan)

        return lowest_clouds

    def _get_itime(self, datetime):
        """
        Returns the index of closest datetime available of the ProfilesData object.
        """
        time = self.data.time.data
        i_time = np.argmin(abs(time - datetime))
        return i_time

    def snr(self, var="attenuated_backscatter_0", step=4, verbose=False):
        """
        Method that calculates the Signal to Noise Ratio.

        Args:
            var (str, optional): Variable of the DataArray to calculate the SNR from.
            step (int, optional): Number of steps around we calculate the SNR for a given altitude.
            verbose (bool, optional): Verbose mode.

        Returns:
            (ProfilesData): object with additional (xarray.DataArray): `snr`.

        Example:
            ```python
            import aprofiles as apro
            # read example file
            path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
            reader = apro.reader.ReadProfiles(path)
            profiles = reader.read()
            # snr calculation
            profiles.snr()
            # snr image
            profiles.plot(var='snr',vmin=0, vmax=3, cmap='Greys_r')
            ```
        """

        # Get the dimensions of the data array
        time_len, altitude_len = self.data[var].shape

        # Preallocate the SNR array with zeros
        snr_array = np.zeros((time_len, altitude_len))

        # Loop over each time step
        for t in track(range(time_len), description="snr   ", disable=not verbose):
            # Extract 1D slice for current time step
            array = self.data[var].data[t, :]

            # Create a sliding window view for the rolling calculation
            sliding_windows = np.lib.stride_tricks.sliding_window_view(array, window_shape=2*step+1, axis=0)

            # Calculate mean and std across the window axis
            means = np.nanmean(sliding_windows, axis=1)
            stds = np.nanstd(sliding_windows, axis=1)

            # Handle the edges (where sliding window can't be applied due to boundary)
            means = np.pad(means, pad_width=(step,), mode='constant', constant_values=np.nan)
            stds = np.pad(stds, pad_width=(step,), mode='constant', constant_values=np.nan)

            # Avoid division by zero
            stds = np.where(stds == 0, np.nan, stds)

            # Compute SNR
            snr_array[t, :] = np.divide(means, stds, where=(stds != 0))

        # Add the SNR DataArray to the object's data attribute
        self.data["snr"] = (('time', 'altitude'), snr_array)
        self.data["snr"] = self.data.snr.assign_attrs({
            'long_name': 'Signal to Noise Ratio',
            'units': '',
            'step': step
        })

        return self

    def gaussian_filter(
        self, sigma=0.25, var="attenuated_backscatter_0", inplace=False
    ):
        """
        Applies a 2D gaussian filter in order to reduce high frequency noise.

        Args:
            sigma (scalar or sequence of scalars, optional): Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.
            var (str, optional): variable name of the Dataset to be processed.
            inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

        Returns:
            (ProfilesData): object with additional attributes `gaussian_filter` for the processed (xarray.DataArray):.

        Example:
            ```python
            import aprofiles as apro
            # read example file
            path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
            reader = apro.reader.ReadProfiles(path)
            profiles = reader.read()
            # apply gaussian filtering
            profiles.gaussian_filter(sigma=0.5, inplace=True)
            profiles.data.attenuated_backscatter_0.attrs.gaussian_filter
            0.50
            ```

            ![Before gaussian filtering](../../assets/images/attenuated_backscatter.png)
            ![After gaussian filtering (sigma=0.5)](../../assets/images/gaussian_filter.png)
        """
        import copy

        from scipy.ndimage import gaussian_filter

        # apply gaussian filter
        filtered_data = gaussian_filter(self.data[var].data, sigma=sigma)

        if inplace:
            self.data[var].data = filtered_data
            new_dataset = self
        else:
            copied_dataset = copy.deepcopy(self)
            copied_dataset.data[var].data = filtered_data
            new_dataset = copied_dataset
        # add attribute
        new_dataset.data[var].attrs["gaussian_filter"] = sigma
        return new_dataset

    def time_avg(self, minutes, var="attenuated_backscatter_0", inplace=False):
        """
        Rolling median in the time dimension.

        Args:
            minutes (float): Number of minutes to average over.
            var (str, optional): variable of the Dataset to be processed.
            inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

        Returns:
            (ProfilesData):
        """
        rcs = self.data[var].copy()
        # time conversion from minutes to seconds
        t_avg = minutes * 60
        # time resolution in profiles data in seconds
        dt_s = self._get_resolution("time")
        # number of timestamps to be to averaged
        nt_avg = max([1, round(t_avg / dt_s)])
        # rolling median
        filtered_data = (
            rcs.rolling(time=nt_avg, min_periods=1, center=True).median().data
        )

        if inplace:
            self.data[var].data = filtered_data
            new_dataset = self
        else:
            copied_dataset = copy.deepcopy(self)
            copied_dataset.data[var].data = filtered_data
            new_dataset = copied_dataset
        # add attribute
        new_dataset.data[var].attrs["time averaged (minutes)"] = minutes
        return new_dataset

    def extrapolate_below(
        self, var="attenuated_backscatter_0", z=150, method="cst", inplace=False
    ):
        """
        Method for extrapolating lowest layers below a certain altitude. This is of particular intrest for instruments subject to After Pulse effect, with saturated signal in the lowest layers.
        We recommend to use a value of zmin=150m due to random values often found below that altitude which perturbs the clouds detection.

        Args:
            var (str, optional): variable of the :class:`xarray.Dataset` to be processed.
            z (float, optional): Altitude (in m, AGL) below which the signal is extrapolated.
            method ({'cst', 'lin'}, optional): Method to be used for extrapolation of lowest layers.
            inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

        Returns:
            (ProfilesData): object with additional attributes:

                - `extrapolation_low_layers_altitude_agl`
                - `extrapolation_low_layers_method` for the processed (xarray.DataArray):.

        Example:
            ```python
            import aprofiles as apro
            # read example file
            path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
            reader = apro.reader.ReadProfiles(path)
            profiles = reader.read()
            # desaturation below 4000m
            profiles.extrapolate_below(z=150., inplace=True)
            profiles.data.attenuated_backscatter_0.extrapolation_low_layers_altitude_agl
            150
            ```

            ![Before extrapolation](../../assets/images/lowest.png)
            ![After extrapolation](../../assets/images/lowest_extrap.png)
        """

        # get index of z
        imax = self._get_index_from_altitude_AGL(z)

        nt = np.shape(self.data[var].data)[0]

        if method == "cst":
            # get values at imin
            data_zmax = self.data[var].data[:, imax]
            # generates ones matrice with time/altitude dimension to fill up bottom
            ones = np.ones((nt, imax))
            # replace values
            filling_matrice = np.transpose(np.multiply(np.transpose(ones), data_zmax))
        elif method == "lin":
            raise NotImplementedError("Linear extrapolation is not implemented yet")
        else:
            raise ValueError("Expected string: lin or cst")

        if inplace:
            self.data[var].data[:, 0:imax] = filling_matrice
            new_profiles_data = self
        else:
            copied_dataset = copy.deepcopy(self)
            copied_dataset.data[var].data[:, 0:imax] = filling_matrice
            new_profiles_data = copied_dataset

        # add attributes
        new_profiles_data.data[var].attrs["extrapolation_low_layers_altitude_agl"] = z
        new_profiles_data.data[var].attrs["extrapolation_low_layers_method"] = method
        return new_profiles_data

    def range_correction(self, var="attenuated_backscatter_0", inplace=False):
        """
        Method that corrects the solid angle effect (1/z2) which makes that the backscatter beam is more unlikely to be detected with the square of the altitude.

        Args:
            var (str, optional): variable of the Dataset to be processed.
            inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

        Returns:
            (ProfilesData):

        .. warning::
            Make sure that the range correction is not already applied to the selected variable.
        """

        # for the altitude correction, must one use the altitude above the ground level
        z_agl = self.data.altitude.data - self.data.station_altitude.data

        data = self.data[var].data
        range_corrected_data = np.multiply(data, z_agl)

        if inplace:
            self.data[var].data = range_corrected_data
            new_profiles_data = self
        else:
            copied_dataset = copy.deepcopy(self)
            copied_dataset.data[var].data = range_corrected_data
            new_profiles_data = copied_dataset

        # add attribute
        new_profiles_data.data[var].attrs["range correction"] = True
        # remove units
        new_profiles_data.data[var].attrs["units"] = None
        return new_profiles_data

    def desaturate_below(self, var="attenuated_backscatter_0", z=4000.0, inplace=False):
        """
        Remove saturation caused by clouds at low altitude which results in negative values above the maximum.
        The absolute value of the signal is returned below the prescribed altitude.

        Args:
            var (str, optional): variable of the :class:`xarray.Dataset` to be processed.
            z (float, optional): Altitude (in m, AGL) below which the signal is unsaturated.
            inplace (bool, optional): if True, replace the instance of the (ProfilesData):.

        Todo:
            Refine method to desaturate only saturated areas.

        Returns:
            (ProfilesData): object with additional attribute `desaturate` for the processed (xarray.DataArray):.

        .. warning::
            For now, absolute values are returned everywhere below the prescribed altitude.

        Example:
            ```python
            import aprofiles as apro
            # read example file
            path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
            reader = apro.reader.ReadProfiles(path)
            profiles = reader.read()
            # desaturation below 4000m
            profiles.desaturate_below(z=4000., inplace=True)
            profiles.data.attenuated_backscatter_0.desaturated
            True
            ```

            ![Before desaturation](../../assets/images/saturated.png)
            ![After desaturation](../../assets/images/desaturated.png)
        """

        imax = self._get_index_from_altitude_AGL(z)
        unsaturated_data = copy.deepcopy(self.data[var].data)
        for i in range(len(self.data.time.data)):
            unsaturated_data[i, :imax] = abs(unsaturated_data[i, :imax])

        if inplace:
            self.data[var].data = unsaturated_data
            new_profiles_data = self
        else:
            copied_dataset = copy.deepcopy(self)
            copied_dataset.data[var].data = unsaturated_data
            new_profiles_data = copied_dataset

        # add attribute
        new_profiles_data.data[var].attrs["desaturated"] = True
        return new_profiles_data

    def foc(self, method="cloud_base", var="attenuated_backscatter_0", z_snr=2000., min_snr=2., zmin_cloud=200.):
        """
        Calls :meth:`aprofiles.detection.foc.detect_foc()`.
        """
        return apro.detection.foc.detect_foc(self, method, var, z_snr, min_snr, zmin_cloud)

    def clouds(self, time_avg=1, zmin=0, thr_noise=5., thr_clouds=4., min_snr=0.0, verbose=False):
        """
        Calls :meth:`aprofiles.detection.clouds.detect_clouds()`.
        """
        return apro.detection.clouds.detect_clouds(self, time_avg, zmin, thr_noise, thr_clouds, min_snr, verbose)

    def pbl(self, time_avg=1, zmin=100., zmax=3000., wav_width=200., under_clouds=True, min_snr=2., verbose=False):
        """
        Calls :meth:`aprofiles.detection.pbl.detect_pbl()`.
        """
        return apro.detection.pbl.detect_pbl(self, time_avg, zmin, zmax, wav_width, under_clouds, min_snr, verbose)

    def inversion(self, time_avg=1, zmin=4000., zmax=6000., min_snr=0., under_clouds=False, method="forward", 
        apriori={"lr": 50.}, remove_outliers=False, mass_conc=True, mass_conc_method="mortier_2013", verbose=False,
    ):
        """
        Calls :meth:`aprofiles.retrieval.extinction.inversion()` to calculate extinction profiles.
        Calls :meth:`aprofiles.retrieval.mass_conc.mec()` to calculate Mass to Extinction coefficients if `mass_conc` is true (Default).
        """
        apro.retrieval.extinction.inversion(self, time_avg, zmin, zmax, min_snr, under_clouds, method, apriori, remove_outliers, verbose)
        if mass_conc:
            apro.retrieval.mass_conc.concentration_profiles(self, mass_conc_method)
        return apro

    def plot(
        self, var="attenuated_backscatter_0", datetime=None, zref="agl", zmin=None, zmax=None, vmin=None, vmax=None, log=False,
        show_foc=False, show_pbl=False, show_clouds=False, cmap="coolwarm", show_fig=True, save_fig=None, **kwargs
    ):
        """
        Plotting method.
        Depending on the variable selected, this method will plot an image, a single profile or a time series of the requested variable.
        See also :ref:`Plotting`.

        Args:
            var (str, optional): Variable to be plotted.
            datetime (:class:`numpy.datetime64`, optional): if provided, plot the profile for closest time. If not, plot an image constructed on all profiles.
            zref ({'agl', 'asl'}, optional): Base reference for the altitude axis.
            zmin (float, optional): Minimum altitude AGL (m).
            zmax (float, optional): Maximum altitude AGL (m).
            vmin (float, optional): Minimum value.
            vmax (float, optional): Maximum value.
            log (bool, optional): Use logarithmic scale.
            show_foc (bool, optional): Show fog or condensation retrievals.
            show_pbl (bool, optional): Show PBL height retrievals.
            show_clouds (bool, optional): Show clouds retrievals.
            cmap (str, optional): Matplotlib colormap.
            show_fig (bool, optional): Show Figure.
            save_fig (str, optional): Path of the saved figure.
        """

        # check if var is available
        if var not in list(self.data.data_vars):
            raise ValueError(
                "{} is not a valid variable. \n List of available variables: {}".format(
                    var, list(self.data.data_vars)
                )
            )

        # here, check the dimension. If the variable has only the time dimension, calls timeseries method
        if datetime is None:
            # check dimension of var
            if len(list(self.data[var].dims)) == 2:
                apro.plot.image.plot(self.data, var, zref, zmin, zmax, vmin, vmax, log, show_foc, show_pbl, show_clouds, cmap, show_fig, save_fig)
            else:
                apro.plot.timeseries.plot(self.data, var, show_fig, save_fig, **kwargs)
        else:
            apro.plot.profile.plot(self.data, datetime, var, zref, zmin, zmax, vmin, vmax, log, show_foc, show_pbl, show_clouds, show_fig, save_fig)

    def write(self, base_dir=Path('examples', 'data', 'V-Profiles'), verbose=False):
        """
        Calls :meth:`aprofiles.io.write_profiles.write()`.
        """
        apro.io.write_profiles.write(self, base_dir, verbose=verbose)

data property writable

Data attribute (instance of (xarray.Dataset):)

clouds(time_avg=1, zmin=0, thr_noise=5.0, thr_clouds=4.0, min_snr=0.0, verbose=False)

Calls :meth:aprofiles.detection.clouds.detect_clouds().

Source code in aprofiles/profiles.py 
393
394
395
396
397
def clouds(self, time_avg=1, zmin=0, thr_noise=5., thr_clouds=4., min_snr=0.0, verbose=False):
    """
    Calls :meth:`aprofiles.detection.clouds.detect_clouds()`.
    """
    return apro.detection.clouds.detect_clouds(self, time_avg, zmin, thr_noise, thr_clouds, min_snr, verbose)

desaturate_below(var='attenuated_backscatter_0', z=4000.0, inplace=False)

Remove saturation caused by clouds at low altitude which results in negative values above the maximum. The absolute value of the signal is returned below the prescribed altitude.

Parameters:

Name Type Description Default
var str

variable of the :class:xarray.Dataset to be processed.

 'attenuated_backscatter_0'
z float

Altitude (in m, AGL) below which the signal is unsaturated.

 4000.0
inplace bool

if True, replace the instance of the (ProfilesData):.

 False
Todo

Refine method to desaturate only saturated areas.

Returns:

Type Description
ProfilesData): object with additional attribute `desaturate` for the processed (xarray.DataArray

.

.. warning:: For now, absolute values are returned everywhere below the prescribed altitude.

Example 
import aprofiles as apro
# read example file
path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
reader = apro.reader.ReadProfiles(path)
profiles = reader.read()
# desaturation below 4000m
profiles.desaturate_below(z=4000., inplace=True)
profiles.data.attenuated_backscatter_0.desaturated
True

Before desaturation After desaturation

Source code in aprofiles/profiles.py 
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
def desaturate_below(self, var="attenuated_backscatter_0", z=4000.0, inplace=False):
    """
    Remove saturation caused by clouds at low altitude which results in negative values above the maximum.
    The absolute value of the signal is returned below the prescribed altitude.

    Args:
        var (str, optional): variable of the :class:`xarray.Dataset` to be processed.
        z (float, optional): Altitude (in m, AGL) below which the signal is unsaturated.
        inplace (bool, optional): if True, replace the instance of the (ProfilesData):.

    Todo:
        Refine method to desaturate only saturated areas.

    Returns:
        (ProfilesData): object with additional attribute `desaturate` for the processed (xarray.DataArray):.

    .. warning::
        For now, absolute values are returned everywhere below the prescribed altitude.

    Example:
        ```python
        import aprofiles as apro
        # read example file
        path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
        reader = apro.reader.ReadProfiles(path)
        profiles = reader.read()
        # desaturation below 4000m
        profiles.desaturate_below(z=4000., inplace=True)
        profiles.data.attenuated_backscatter_0.desaturated
        True
        ```

        ![Before desaturation](../../assets/images/saturated.png)
        ![After desaturation](../../assets/images/desaturated.png)
    """

    imax = self._get_index_from_altitude_AGL(z)
    unsaturated_data = copy.deepcopy(self.data[var].data)
    for i in range(len(self.data.time.data)):
        unsaturated_data[i, :imax] = abs(unsaturated_data[i, :imax])

    if inplace:
        self.data[var].data = unsaturated_data
        new_profiles_data = self
    else:
        copied_dataset = copy.deepcopy(self)
        copied_dataset.data[var].data = unsaturated_data
        new_profiles_data = copied_dataset

    # add attribute
    new_profiles_data.data[var].attrs["desaturated"] = True
    return new_profiles_data

extrapolate_below(var='attenuated_backscatter_0', z=150, method='cst', inplace=False)

Method for extrapolating lowest layers below a certain altitude. This is of particular intrest for instruments subject to After Pulse effect, with saturated signal in the lowest layers. We recommend to use a value of zmin=150m due to random values often found below that altitude which perturbs the clouds detection.

Parameters:

Name Type Description Default
var str

variable of the :class:xarray.Dataset to be processed.

 'attenuated_backscatter_0'
z float

Altitude (in m, AGL) below which the signal is extrapolated.

 150
method {cst, lin}

Method to be used for extrapolation of lowest layers.

 'cst'
inplace bool

if True, replace the instance of the (ProfilesData): class.

 False

Returns:

Type Description
ProfilesData

object with additional attributes:

  • extrapolation_low_layers_altitude_agl
  • extrapolation_low_layers_method for the processed (xarray.DataArray):.
Example 
import aprofiles as apro
# read example file
path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
reader = apro.reader.ReadProfiles(path)
profiles = reader.read()
# desaturation below 4000m
profiles.extrapolate_below(z=150., inplace=True)
profiles.data.attenuated_backscatter_0.extrapolation_low_layers_altitude_agl
150

Before extrapolation After extrapolation

Source code in aprofiles/profiles.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
def extrapolate_below(
    self, var="attenuated_backscatter_0", z=150, method="cst", inplace=False
):
    """
    Method for extrapolating lowest layers below a certain altitude. This is of particular intrest for instruments subject to After Pulse effect, with saturated signal in the lowest layers.
    We recommend to use a value of zmin=150m due to random values often found below that altitude which perturbs the clouds detection.

    Args:
        var (str, optional): variable of the :class:`xarray.Dataset` to be processed.
        z (float, optional): Altitude (in m, AGL) below which the signal is extrapolated.
        method ({'cst', 'lin'}, optional): Method to be used for extrapolation of lowest layers.
        inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

    Returns:
        (ProfilesData): object with additional attributes:

            - `extrapolation_low_layers_altitude_agl`
            - `extrapolation_low_layers_method` for the processed (xarray.DataArray):.

    Example:
        ```python
        import aprofiles as apro
        # read example file
        path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
        reader = apro.reader.ReadProfiles(path)
        profiles = reader.read()
        # desaturation below 4000m
        profiles.extrapolate_below(z=150., inplace=True)
        profiles.data.attenuated_backscatter_0.extrapolation_low_layers_altitude_agl
        150
        ```

        ![Before extrapolation](../../assets/images/lowest.png)
        ![After extrapolation](../../assets/images/lowest_extrap.png)
    """

    # get index of z
    imax = self._get_index_from_altitude_AGL(z)

    nt = np.shape(self.data[var].data)[0]

    if method == "cst":
        # get values at imin
        data_zmax = self.data[var].data[:, imax]
        # generates ones matrice with time/altitude dimension to fill up bottom
        ones = np.ones((nt, imax))
        # replace values
        filling_matrice = np.transpose(np.multiply(np.transpose(ones), data_zmax))
    elif method == "lin":
        raise NotImplementedError("Linear extrapolation is not implemented yet")
    else:
        raise ValueError("Expected string: lin or cst")

    if inplace:
        self.data[var].data[:, 0:imax] = filling_matrice
        new_profiles_data = self
    else:
        copied_dataset = copy.deepcopy(self)
        copied_dataset.data[var].data[:, 0:imax] = filling_matrice
        new_profiles_data = copied_dataset

    # add attributes
    new_profiles_data.data[var].attrs["extrapolation_low_layers_altitude_agl"] = z
    new_profiles_data.data[var].attrs["extrapolation_low_layers_method"] = method
    return new_profiles_data

foc(method='cloud_base', var='attenuated_backscatter_0', z_snr=2000.0, min_snr=2.0, zmin_cloud=200.0)

Calls :meth:aprofiles.detection.foc.detect_foc().

Source code in aprofiles/profiles.py 
387
388
389
390
391
def foc(self, method="cloud_base", var="attenuated_backscatter_0", z_snr=2000., min_snr=2., zmin_cloud=200.):
    """
    Calls :meth:`aprofiles.detection.foc.detect_foc()`.
    """
    return apro.detection.foc.detect_foc(self, method, var, z_snr, min_snr, zmin_cloud)

gaussian_filter(sigma=0.25, var='attenuated_backscatter_0', inplace=False)

Applies a 2D gaussian filter in order to reduce high frequency noise.

Parameters:

Name Type Description Default
sigma scalar or sequence of scalars

Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.

 0.25
var str

variable name of the Dataset to be processed.

 'attenuated_backscatter_0'
inplace bool

if True, replace the instance of the (ProfilesData): class.

 False

Returns:

Type Description
ProfilesData): object with additional attributes `gaussian_filter` for the processed (xarray.DataArray

.
Example 
import aprofiles as apro
# read example file
path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
reader = apro.reader.ReadProfiles(path)
profiles = reader.read()
# apply gaussian filtering
profiles.gaussian_filter(sigma=0.5, inplace=True)
profiles.data.attenuated_backscatter_0.attrs.gaussian_filter
0.50

Before gaussian filtering After gaussian filtering (sigma=0.5)

Source code in aprofiles/profiles.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def gaussian_filter(
    self, sigma=0.25, var="attenuated_backscatter_0", inplace=False
):
    """
    Applies a 2D gaussian filter in order to reduce high frequency noise.

    Args:
        sigma (scalar or sequence of scalars, optional): Standard deviation for Gaussian kernel. The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.
        var (str, optional): variable name of the Dataset to be processed.
        inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

    Returns:
        (ProfilesData): object with additional attributes `gaussian_filter` for the processed (xarray.DataArray):.

    Example:
        ```python
        import aprofiles as apro
        # read example file
        path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
        reader = apro.reader.ReadProfiles(path)
        profiles = reader.read()
        # apply gaussian filtering
        profiles.gaussian_filter(sigma=0.5, inplace=True)
        profiles.data.attenuated_backscatter_0.attrs.gaussian_filter
        0.50
        ```

        ![Before gaussian filtering](../../assets/images/attenuated_backscatter.png)
        ![After gaussian filtering (sigma=0.5)](../../assets/images/gaussian_filter.png)
    """
    import copy

    from scipy.ndimage import gaussian_filter

    # apply gaussian filter
    filtered_data = gaussian_filter(self.data[var].data, sigma=sigma)

    if inplace:
        self.data[var].data = filtered_data
        new_dataset = self
    else:
        copied_dataset = copy.deepcopy(self)
        copied_dataset.data[var].data = filtered_data
        new_dataset = copied_dataset
    # add attribute
    new_dataset.data[var].attrs["gaussian_filter"] = sigma
    return new_dataset

inversion(time_avg=1, zmin=4000.0, zmax=6000.0, min_snr=0.0, under_clouds=False, method='forward', apriori={'lr': 50.0}, remove_outliers=False, mass_conc=True, mass_conc_method='mortier_2013', verbose=False)

Calls :meth:aprofiles.retrieval.extinction.inversion() to calculate extinction profiles. Calls :meth:aprofiles.retrieval.mass_conc.mec() to calculate Mass to Extinction coefficients if mass_conc is true (Default).

Source code in aprofiles/profiles.py 
405
406
407
408
409
410
411
412
413
414
415
def inversion(self, time_avg=1, zmin=4000., zmax=6000., min_snr=0., under_clouds=False, method="forward", 
    apriori={"lr": 50.}, remove_outliers=False, mass_conc=True, mass_conc_method="mortier_2013", verbose=False,
):
    """
    Calls :meth:`aprofiles.retrieval.extinction.inversion()` to calculate extinction profiles.
    Calls :meth:`aprofiles.retrieval.mass_conc.mec()` to calculate Mass to Extinction coefficients if `mass_conc` is true (Default).
    """
    apro.retrieval.extinction.inversion(self, time_avg, zmin, zmax, min_snr, under_clouds, method, apriori, remove_outliers, verbose)
    if mass_conc:
        apro.retrieval.mass_conc.concentration_profiles(self, mass_conc_method)
    return apro

pbl(time_avg=1, zmin=100.0, zmax=3000.0, wav_width=200.0, under_clouds=True, min_snr=2.0, verbose=False)

Calls :meth:aprofiles.detection.pbl.detect_pbl().

Source code in aprofiles/profiles.py 
399
400
401
402
403
def pbl(self, time_avg=1, zmin=100., zmax=3000., wav_width=200., under_clouds=True, min_snr=2., verbose=False):
    """
    Calls :meth:`aprofiles.detection.pbl.detect_pbl()`.
    """
    return apro.detection.pbl.detect_pbl(self, time_avg, zmin, zmax, wav_width, under_clouds, min_snr, verbose)

plot(var='attenuated_backscatter_0', datetime=None, zref='agl', zmin=None, zmax=None, vmin=None, vmax=None, log=False, show_foc=False, show_pbl=False, show_clouds=False, cmap='coolwarm', show_fig=True, save_fig=None, **kwargs)

Plotting method. Depending on the variable selected, this method will plot an image, a single profile or a time series of the requested variable. See also :ref:Plotting.

Parameters:

Name Type Description Default
var str

Variable to be plotted.

 'attenuated_backscatter_0'
datetime

class:numpy.datetime64, optional): if provided, plot the profile for closest time. If not, plot an image constructed on all profiles.

 None
zref {agl, asl}

Base reference for the altitude axis.

 'agl'
zmin float

Minimum altitude AGL (m).

 None
zmax float

Maximum altitude AGL (m).

 None
vmin float

Minimum value.

 None
vmax float

Maximum value.

 None
log bool

Use logarithmic scale.

 False
show_foc bool

Show fog or condensation retrievals.

 False
show_pbl bool

Show PBL height retrievals.

 False
show_clouds bool

Show clouds retrievals.

 False
cmap str

Matplotlib colormap.

 'coolwarm'
show_fig bool

Show Figure.

 True
save_fig str

Path of the saved figure.

 None
Source code in aprofiles/profiles.py 
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
def plot(
    self, var="attenuated_backscatter_0", datetime=None, zref="agl", zmin=None, zmax=None, vmin=None, vmax=None, log=False,
    show_foc=False, show_pbl=False, show_clouds=False, cmap="coolwarm", show_fig=True, save_fig=None, **kwargs
):
    """
    Plotting method.
    Depending on the variable selected, this method will plot an image, a single profile or a time series of the requested variable.
    See also :ref:`Plotting`.

    Args:
        var (str, optional): Variable to be plotted.
        datetime (:class:`numpy.datetime64`, optional): if provided, plot the profile for closest time. If not, plot an image constructed on all profiles.
        zref ({'agl', 'asl'}, optional): Base reference for the altitude axis.
        zmin (float, optional): Minimum altitude AGL (m).
        zmax (float, optional): Maximum altitude AGL (m).
        vmin (float, optional): Minimum value.
        vmax (float, optional): Maximum value.
        log (bool, optional): Use logarithmic scale.
        show_foc (bool, optional): Show fog or condensation retrievals.
        show_pbl (bool, optional): Show PBL height retrievals.
        show_clouds (bool, optional): Show clouds retrievals.
        cmap (str, optional): Matplotlib colormap.
        show_fig (bool, optional): Show Figure.
        save_fig (str, optional): Path of the saved figure.
    """

    # check if var is available
    if var not in list(self.data.data_vars):
        raise ValueError(
            "{} is not a valid variable. \n List of available variables: {}".format(
                var, list(self.data.data_vars)
            )
        )

    # here, check the dimension. If the variable has only the time dimension, calls timeseries method
    if datetime is None:
        # check dimension of var
        if len(list(self.data[var].dims)) == 2:
            apro.plot.image.plot(self.data, var, zref, zmin, zmax, vmin, vmax, log, show_foc, show_pbl, show_clouds, cmap, show_fig, save_fig)
        else:
            apro.plot.timeseries.plot(self.data, var, show_fig, save_fig, **kwargs)
    else:
        apro.plot.profile.plot(self.data, datetime, var, zref, zmin, zmax, vmin, vmax, log, show_foc, show_pbl, show_clouds, show_fig, save_fig)

range_correction(var='attenuated_backscatter_0', inplace=False)

Method that corrects the solid angle effect (1/z2) which makes that the backscatter beam is more unlikely to be detected with the square of the altitude.

Parameters:

Name Type Description Default
var str

variable of the Dataset to be processed.

 'attenuated_backscatter_0'
inplace bool

if True, replace the instance of the (ProfilesData): class.

 False

Returns:

Type Description
ProfilesData

.. warning:: Make sure that the range correction is not already applied to the selected variable.

Source code in aprofiles/profiles.py 
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
def range_correction(self, var="attenuated_backscatter_0", inplace=False):
    """
    Method that corrects the solid angle effect (1/z2) which makes that the backscatter beam is more unlikely to be detected with the square of the altitude.

    Args:
        var (str, optional): variable of the Dataset to be processed.
        inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

    Returns:
        (ProfilesData):

    .. warning::
        Make sure that the range correction is not already applied to the selected variable.
    """

    # for the altitude correction, must one use the altitude above the ground level
    z_agl = self.data.altitude.data - self.data.station_altitude.data

    data = self.data[var].data
    range_corrected_data = np.multiply(data, z_agl)

    if inplace:
        self.data[var].data = range_corrected_data
        new_profiles_data = self
    else:
        copied_dataset = copy.deepcopy(self)
        copied_dataset.data[var].data = range_corrected_data
        new_profiles_data = copied_dataset

    # add attribute
    new_profiles_data.data[var].attrs["range correction"] = True
    # remove units
    new_profiles_data.data[var].attrs["units"] = None
    return new_profiles_data

snr(var='attenuated_backscatter_0', step=4, verbose=False)

Method that calculates the Signal to Noise Ratio.

Parameters:

Name Type Description Default
var str

Variable of the DataArray to calculate the SNR from.

 'attenuated_backscatter_0'
step int

Number of steps around we calculate the SNR for a given altitude.

 4
verbose bool

Verbose mode.

 False

Returns:

Type Description
ProfilesData): object with additional (xarray.DataArray

snr.
Example 
import aprofiles as apro
# read example file
path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
reader = apro.reader.ReadProfiles(path)
profiles = reader.read()
# snr calculation
profiles.snr()
# snr image
profiles.plot(var='snr',vmin=0, vmax=3, cmap='Greys_r')
Source code in aprofiles/profiles.py 
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def snr(self, var="attenuated_backscatter_0", step=4, verbose=False):
    """
    Method that calculates the Signal to Noise Ratio.

    Args:
        var (str, optional): Variable of the DataArray to calculate the SNR from.
        step (int, optional): Number of steps around we calculate the SNR for a given altitude.
        verbose (bool, optional): Verbose mode.

    Returns:
        (ProfilesData): object with additional (xarray.DataArray): `snr`.

    Example:
        ```python
        import aprofiles as apro
        # read example file
        path = "examples/data/E-PROFILE/L2_0-20000-001492_A20210909.nc"
        reader = apro.reader.ReadProfiles(path)
        profiles = reader.read()
        # snr calculation
        profiles.snr()
        # snr image
        profiles.plot(var='snr',vmin=0, vmax=3, cmap='Greys_r')
        ```
    """

    # Get the dimensions of the data array
    time_len, altitude_len = self.data[var].shape

    # Preallocate the SNR array with zeros
    snr_array = np.zeros((time_len, altitude_len))

    # Loop over each time step
    for t in track(range(time_len), description="snr   ", disable=not verbose):
        # Extract 1D slice for current time step
        array = self.data[var].data[t, :]

        # Create a sliding window view for the rolling calculation
        sliding_windows = np.lib.stride_tricks.sliding_window_view(array, window_shape=2*step+1, axis=0)

        # Calculate mean and std across the window axis
        means = np.nanmean(sliding_windows, axis=1)
        stds = np.nanstd(sliding_windows, axis=1)

        # Handle the edges (where sliding window can't be applied due to boundary)
        means = np.pad(means, pad_width=(step,), mode='constant', constant_values=np.nan)
        stds = np.pad(stds, pad_width=(step,), mode='constant', constant_values=np.nan)

        # Avoid division by zero
        stds = np.where(stds == 0, np.nan, stds)

        # Compute SNR
        snr_array[t, :] = np.divide(means, stds, where=(stds != 0))

    # Add the SNR DataArray to the object's data attribute
    self.data["snr"] = (('time', 'altitude'), snr_array)
    self.data["snr"] = self.data.snr.assign_attrs({
        'long_name': 'Signal to Noise Ratio',
        'units': '',
        'step': step
    })

    return self

time_avg(minutes, var='attenuated_backscatter_0', inplace=False)

Rolling median in the time dimension.

Parameters:

Name Type Description Default
minutes float

Number of minutes to average over.

 required
var str

variable of the Dataset to be processed.

 'attenuated_backscatter_0'
inplace bool

if True, replace the instance of the (ProfilesData): class.

 False

Returns:

Type Description
ProfilesData
Source code in aprofiles/profiles.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def time_avg(self, minutes, var="attenuated_backscatter_0", inplace=False):
    """
    Rolling median in the time dimension.

    Args:
        minutes (float): Number of minutes to average over.
        var (str, optional): variable of the Dataset to be processed.
        inplace (bool, optional): if True, replace the instance of the (ProfilesData): class.

    Returns:
        (ProfilesData):
    """
    rcs = self.data[var].copy()
    # time conversion from minutes to seconds
    t_avg = minutes * 60
    # time resolution in profiles data in seconds
    dt_s = self._get_resolution("time")
    # number of timestamps to be to averaged
    nt_avg = max([1, round(t_avg / dt_s)])
    # rolling median
    filtered_data = (
        rcs.rolling(time=nt_avg, min_periods=1, center=True).median().data
    )

    if inplace:
        self.data[var].data = filtered_data
        new_dataset = self
    else:
        copied_dataset = copy.deepcopy(self)
        copied_dataset.data[var].data = filtered_data
        new_dataset = copied_dataset
    # add attribute
    new_dataset.data[var].attrs["time averaged (minutes)"] = minutes
    return new_dataset

write(base_dir=Path('examples', 'data', 'V-Profiles'), verbose=False)

Calls :meth:aprofiles.io.write_profiles.write().

Source code in aprofiles/profiles.py 
461
462
463
464
465
def write(self, base_dir=Path('examples', 'data', 'V-Profiles'), verbose=False):
    """
    Calls :meth:`aprofiles.io.write_profiles.write()`.
    """
    apro.io.write_profiles.write(self, base_dir, verbose=verbose)

Aeronet

Not implemented yet.

AeronetData

Base class representing profiles data returned by (aprofiles.io.reader.ReadAeronet):.

Source code in aprofiles/aeronet.py 
 5
 6
 7
 8
 9
10
11
class AeronetData:
    """
    Base class representing profiles data returned by (aprofiles.io.reader.ReadAeronet):."""

    def __init__(self, data):
        self.data = data
        raise NotImplementedError("AeronetData is not implemented yet")

Rayleigh

The RayleighData class is used for producing Rayleigh profiles in a standard atmosphere.

RayleighData

Class for computing rayleigh profile in a standard atmosphere. This class calls the :func:get_optics_in_std_atmo() method, which produces profiles of backscatter and extinction coefficients.

Attributes:

Name Type Description
altitude array - like

array of altitude ASL to be used to compute the rayleigh profile, in m.
wavelength float

Wavelength of the Rayleigh profile to be computed, in nm.
T0 float

Temperature at ground level, in K.
P0 float

Pressure at ground level, in hPa.
Example 
# some imports
import aprofiles as apro
import numpy as np
# creates altitude array
altitude = np.arange(15,15000,15)
wavelength = 1064.
# produce rayleigh profile
rayleigh = apro.rayleigh.RayleighData(altitude, wavelength, T0=298, P0=1013);
# checkout the instance attributes
rayleigh.__dict__.keys()
dict_keys(['altitude', 'T0', 'P0', 'wavelength', 'cross_section', 'backscatter', 'extinction'])
Source code in aprofiles/rayleigh.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
class RayleighData:
    """
    Class for computing *rayleigh profile* in a standard atmosphere.
    This class calls the :func:`get_optics_in_std_atmo()` method, which produces profiles of `backscatter` and `extinction` coefficients.

    Attributes:
        altitude (array-like): array of altitude ASL to be used to compute the rayleigh profile, in m.
        wavelength (float): Wavelength of the Rayleigh profile to be computed, in nm.
        T0 (float): Temperature at ground level, in K.
        P0 (float): Pressure at ground level, in hPa.

    Example:
        ```python
        # some imports
        import aprofiles as apro
        import numpy as np
        # creates altitude array
        altitude = np.arange(15,15000,15)
        wavelength = 1064.
        # produce rayleigh profile
        rayleigh = apro.rayleigh.RayleighData(altitude, wavelength, T0=298, P0=1013);
        # checkout the instance attributes
        rayleigh.__dict__.keys()
        dict_keys(['altitude', 'T0', 'P0', 'wavelength', 'cross_section', 'backscatter', 'extinction'])
        ```
    """

    def __init__(self, altitude: list, wavelength: float, T0=298, P0=1013):
        self.altitude = altitude
        self.T0 = T0
        self.P0 = P0
        self.wavelength = wavelength

        # calls functions
        self.get_optics_in_std_atmo()

    def get_optics_in_std_atmo(self):
        """
        Function that returns *backscatter* and *extinction* profiles [#]_ for an instance of a (RayleighData): class.

        .. [#] Bucholtz, A. (1995). Rayleigh-scattering calculations for the terrestrial atmosphere. Applied optics, 34(15), 2765-2773.

        Returns:
            (aprofiles.rayleigh.RayleighData): object with additional attributes:

                - `extinction` (numpy.ndarray): extinction coefficient (m-1)
                - `backscatter` (numpy.ndarray): backscatter coefficient (m-1.sr-1)
                - `cross_section` (numpy.ndarray): cross section (cm-2)
        """

        def _refi_air(wavelength):
            """
            Function that calculates the refractive index of the air at a given wavelength in a standard atmosphere [#]_.

            .. [#] Peck, E. R., & Reeder, K. (1972). Dispersion of air. JOSA, 62(8), 958-962.

            Args:
                wavelength (float): wavelength, in µm.

            Returns:
                refractive index of the air.
            """
            # wav (float): wavelength in µm
            # returns refractive index of the air in standard atmosphere (Peck and Reeder)
            var = (
                8060.51
                + 2480990 / (132.274 - wavelength ** -2)
                + 17455.7 / (39.32957 - wavelength ** -2)
            )
            return var * 1e-8 + 1

        # standard values at sea level
        T0 = 298  # K
        P0 = 1013 * 1e2  # Pa
        p_He = 8

        # standard gradients & parameters
        atmo = {
            "tropo": {"zmin": 0, "zmax": 13, "dTdz": -6.5},
            "strato": {"zmin": 13, "zmax": 55, "dTdz": 1.4},
            "meso": {"zmin": 55, "zmax": 100, "dTdz": -2.4},
        }

        # convert altitude to km
        z = self.altitude / 1000

        # vertical resolution
        dz = min(z[1:] - z[0:-1])

        # temperature profile
        Tz = [T0]
        for layer in atmo.keys():
            ibottom = int(np.floor(atmo[layer]["zmin"] / dz))
            itop = int(np.floor(atmo[layer]["zmax"] / dz))
            for i in range(ibottom, itop):
                Tz.append(Tz[ibottom] + (i - ibottom) * atmo[layer]["dTdz"] * dz)

        # vertical coordinates for temperature
        z_Tz = np.arange(0, 100, dz)

        # pressure profile
        Pz = P0 * np.exp(-z_Tz / p_He)

        # molecules density and cross section
        Ns = (6.02214e23 / 22.4141) / 1000  # molecules.cm-3
        # density (cm-3)
        N_m = [Ns * (T0 / P0) * (Pz[i] / Tz[i]) for i in range(len(Pz))]

        # cross section, in cm2 (adapted from Bodhaine et al., 1999)
        king_factor = 1.05  # tomasi et al., 2005
        num = 24 * (np.pi ** 3) * ((_refi_air(self.wavelength * 1e-3) ** 2 - 1) ** 2)
        denum = (
            ((self.wavelength * 1e-7) ** 4)
            * (Ns ** 2)
            * ((_refi_air(self.wavelength * 1e-3) ** 2 + 2) ** 2)
        )
        section_m = (num / denum) * king_factor

        # extinction profile
        amol = N_m * np.array(section_m)
        # backscatter profile
        lr_mol = 8 * np.pi / 3
        bmol = amol / lr_mol

        # colocate vertically to input altitude
        imin = np.argmin(np.abs(np.array(z_Tz) - z[0]))
        imax = imin + len(z)

        # output
        self.cross_section = section_m  # in cm2
        self.backscatter = bmol[imin:imax] * 1e2  # from cm-1 to m-1
        self.extinction = amol[imin:imax] * 1e2  # from cm-1 to m-1
        self.tau = np.cumsum(self.extinction*(dz*1000))[-1]

        return self

    def plot(self, show_fig=True):
        """
        Plot extinction profile of an instance of the (RayleighData): class.

        Example:
            ```python
            # some imports
            import aprofiles as apro
            import numpy as np
            # creates altitude array
            altitude = np.arange(15,15000,15)
            wavelength = 1064.
            # produce rayleigh profile
            rayleigh = apro.rayleigh.RayleighData(altitude, wavelength, T0=298, P0=1013);
            # plot profile
            rayleigh.plot()
            ```

            .. figure:: ../../docs/assets/images/rayleigh.png
                :scale: 80 %
                :alt: rayleigh profile

                Rayleigh extinction profile for a standard atmosphere.
        """

        fig, ax = plt.subplots(1, 1, figsize=(6, 6))
        plt.plot(self.extinction, self.altitude)
        plt.text(
            0.97,
            0.94,
            f"$\sigma_m: {self.cross_section:.2e} cm2$",
            horizontalalignment="right",
            verticalalignment="center",
            transform=ax.transAxes,
        )
        plt.text(
            0.97,
            0.88,
            f"$\u03C4_m: {self.tau:.2e}$",
            horizontalalignment="right",
            verticalalignment="center",
            transform=ax.transAxes,
        )

        plt.title(
            f"Rayleigh Profile in a Standard Atmosphere ({self.P0}hPa, {self.T0}K)",
            weight="bold",
        )
        plt.xlabel(f"Extinction coefficient @ {self.wavelength} nm (m-1)")
        plt.ylabel("Altitude ASL (m)")
        plt.tight_layout()
        if show_fig:
            plt.show()

get_optics_in_std_atmo()

Function that returns backscatter and extinction profiles [#]_ for an instance of a (RayleighData): class.

.. [#] Bucholtz, A. (1995). Rayleigh-scattering calculations for the terrestrial atmosphere. Applied optics, 34(15), 2765-2773.

Returns:

Type Description
RayleighData

object with additional attributes:

  • extinction (numpy.ndarray): extinction coefficient (m-1)
  • backscatter (numpy.ndarray): backscatter coefficient (m-1.sr-1)
  • cross_section (numpy.ndarray): cross section (cm-2)
Source code in aprofiles/rayleigh.py 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
def get_optics_in_std_atmo(self):
    """
    Function that returns *backscatter* and *extinction* profiles [#]_ for an instance of a (RayleighData): class.

    .. [#] Bucholtz, A. (1995). Rayleigh-scattering calculations for the terrestrial atmosphere. Applied optics, 34(15), 2765-2773.

    Returns:
        (aprofiles.rayleigh.RayleighData): object with additional attributes:

            - `extinction` (numpy.ndarray): extinction coefficient (m-1)
            - `backscatter` (numpy.ndarray): backscatter coefficient (m-1.sr-1)
            - `cross_section` (numpy.ndarray): cross section (cm-2)
    """

    def _refi_air(wavelength):
        """
        Function that calculates the refractive index of the air at a given wavelength in a standard atmosphere [#]_.

        .. [#] Peck, E. R., & Reeder, K. (1972). Dispersion of air. JOSA, 62(8), 958-962.

        Args:
            wavelength (float): wavelength, in µm.

        Returns:
            refractive index of the air.
        """
        # wav (float): wavelength in µm
        # returns refractive index of the air in standard atmosphere (Peck and Reeder)
        var = (
            8060.51
            + 2480990 / (132.274 - wavelength ** -2)
            + 17455.7 / (39.32957 - wavelength ** -2)
        )
        return var * 1e-8 + 1

    # standard values at sea level
    T0 = 298  # K
    P0 = 1013 * 1e2  # Pa
    p_He = 8

    # standard gradients & parameters
    atmo = {
        "tropo": {"zmin": 0, "zmax": 13, "dTdz": -6.5},
        "strato": {"zmin": 13, "zmax": 55, "dTdz": 1.4},
        "meso": {"zmin": 55, "zmax": 100, "dTdz": -2.4},
    }

    # convert altitude to km
    z = self.altitude / 1000

    # vertical resolution
    dz = min(z[1:] - z[0:-1])

    # temperature profile
    Tz = [T0]
    for layer in atmo.keys():
        ibottom = int(np.floor(atmo[layer]["zmin"] / dz))
        itop = int(np.floor(atmo[layer]["zmax"] / dz))
        for i in range(ibottom, itop):
            Tz.append(Tz[ibottom] + (i - ibottom) * atmo[layer]["dTdz"] * dz)

    # vertical coordinates for temperature
    z_Tz = np.arange(0, 100, dz)

    # pressure profile
    Pz = P0 * np.exp(-z_Tz / p_He)

    # molecules density and cross section
    Ns = (6.02214e23 / 22.4141) / 1000  # molecules.cm-3
    # density (cm-3)
    N_m = [Ns * (T0 / P0) * (Pz[i] / Tz[i]) for i in range(len(Pz))]

    # cross section, in cm2 (adapted from Bodhaine et al., 1999)
    king_factor = 1.05  # tomasi et al., 2005
    num = 24 * (np.pi ** 3) * ((_refi_air(self.wavelength * 1e-3) ** 2 - 1) ** 2)
    denum = (
        ((self.wavelength * 1e-7) ** 4)
        * (Ns ** 2)
        * ((_refi_air(self.wavelength * 1e-3) ** 2 + 2) ** 2)
    )
    section_m = (num / denum) * king_factor

    # extinction profile
    amol = N_m * np.array(section_m)
    # backscatter profile
    lr_mol = 8 * np.pi / 3
    bmol = amol / lr_mol

    # colocate vertically to input altitude
    imin = np.argmin(np.abs(np.array(z_Tz) - z[0]))
    imax = imin + len(z)

    # output
    self.cross_section = section_m  # in cm2
    self.backscatter = bmol[imin:imax] * 1e2  # from cm-1 to m-1
    self.extinction = amol[imin:imax] * 1e2  # from cm-1 to m-1
    self.tau = np.cumsum(self.extinction*(dz*1000))[-1]

    return self

plot(show_fig=True)

Plot extinction profile of an instance of the (RayleighData): class.

Example 
# some imports
import aprofiles as apro
import numpy as np
# creates altitude array
altitude = np.arange(15,15000,15)
wavelength = 1064.
# produce rayleigh profile
rayleigh = apro.rayleigh.RayleighData(altitude, wavelength, T0=298, P0=1013);
# plot profile
rayleigh.plot()

.. figure:: ../../docs/assets/images/rayleigh.png :scale: 80 % :alt: rayleigh profile

Rayleigh extinction profile for a standard atmosphere.
Source code in aprofiles/rayleigh.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def plot(self, show_fig=True):
    """
    Plot extinction profile of an instance of the (RayleighData): class.

    Example:
        ```python
        # some imports
        import aprofiles as apro
        import numpy as np
        # creates altitude array
        altitude = np.arange(15,15000,15)
        wavelength = 1064.
        # produce rayleigh profile
        rayleigh = apro.rayleigh.RayleighData(altitude, wavelength, T0=298, P0=1013);
        # plot profile
        rayleigh.plot()
        ```

        .. figure:: ../../docs/assets/images/rayleigh.png
            :scale: 80 %
            :alt: rayleigh profile

            Rayleigh extinction profile for a standard atmosphere.
    """

    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
    plt.plot(self.extinction, self.altitude)
    plt.text(
        0.97,
        0.94,
        f"$\sigma_m: {self.cross_section:.2e} cm2$",
        horizontalalignment="right",
        verticalalignment="center",
        transform=ax.transAxes,
    )
    plt.text(
        0.97,
        0.88,
        f"$\u03C4_m: {self.tau:.2e}$",
        horizontalalignment="right",
        verticalalignment="center",
        transform=ax.transAxes,
    )

    plt.title(
        f"Rayleigh Profile in a Standard Atmosphere ({self.P0}hPa, {self.T0}K)",
        weight="bold",
    )
    plt.xlabel(f"Extinction coefficient @ {self.wavelength} nm (m-1)")
    plt.ylabel("Altitude ASL (m)")
    plt.tight_layout()
    if show_fig:
        plt.show()

Size Distribution

The size distribution module is used to produce volume and number size distributions of a population of particles for a given type. The values describing the size distribution for different aerosol types are taken from the literature:

  • dust, biomass_burning, and urban aerosols1
  • volcanic_ash2

Aerosol properties are defined in config/aer_properties.json

SizeDistributionData

Class for computing size distributions for a given aerosol type. This class calls the aprofiles.size_distribution.SizeDistributionData.get_sd: method, which calculates VSD (Volume Size Distribution) and NSD (Number Size Distribution).

Attributes:

Name Type Description
`aer_type` {dust, volcanic_ash, biomass_burning, urban}

aerosol type.
`aer_properties` dict

dictionnary describing the optical and microphophysical properties of the prescribed aerosol (read from aer_properties.json)
Example 
# some imports
import aprofiles as apro
sd = apro.size_distribution.SizeDistributionData('urban')
# checkout the instance attributes
apro.size_distribution.SizeDistributionData('dust').__dict__.keys()
dict_keys(['aer_type', 'aer_properties', 'radius', 'vsd', 'nsd'])
Source code in aprofiles/size_distribution.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
class SizeDistributionData:
    """
    Class for computing *size distributions* for a given aerosol type.
    This class calls the `aprofiles.size_distribution.SizeDistributionData.get_sd`: method, which calculates VSD (Volume Size Distribution) and NSD (Number Size Distribution).

    Attributes:
        `aer_type` ({'dust', 'volcanic_ash', 'biomass_burning','urban'}): aerosol type.
        `aer_properties` (dict): dictionnary describing the optical and microphophysical properties of the prescribed aerosol (read from *aer_properties.json*)

    Example:
        ```python
        # some imports
        import aprofiles as apro
        sd = apro.size_distribution.SizeDistributionData('urban')
        # checkout the instance attributes
        apro.size_distribution.SizeDistributionData('dust').__dict__.keys()
        dict_keys(['aer_type', 'aer_properties', 'radius', 'vsd', 'nsd'])
        ```
    """

    def __init__(self, aer_type):
        self.aer_type = aer_type

        # read aer_properties.json files
        f = open(Path(Path(__file__).parent,'config','aer_properties.json'))
        aer_properties = json.load(f)
        f.close()
        # check if the aer_type exist in the json file
        if aer_type not in aer_properties.keys():
            raise ValueError(
                f"{aer_type} not found in aer_properties.json. `aer_type` must be one of the follwowing: {list(aer_properties.keys())}"
            )
        else:
            self.aer_properties = aer_properties[self.aer_type]
            self.get_sd()

    def _vsd_to_nsd(self):
        """
        Transforms Volume Size Distribution to Number Size Distribution
        """
        self.nsd = [
            self.vsd[i] * 3 / (4 * np.pi * self.radius[i] ** 4)
            for i in range(len(self.radius))
        ]
        return self

    def get_sd(self):
        """
        Returns the Volume and Number Size Distributions arrays from an instance of the (SizeDistributionData): class .

        Returns:
            (SizeDistribData): object with additional attributes:

                - `radius` (numpy.ndarray): radius (µm)
                - `vsd` (numpy.ndarray): Volume Size Distribution (µm2.µm-3)
                - `nsd` (numpy.ndarray): Number Size Distribution (µm-3.µm-1)
        """

        aer_properties = self.aer_properties
        radius = np.arange(1e-2, 2e1, 1e-3)  # radius in µm
        vsd = np.zeros(len(radius))

        # we loop though all the keys defining the different modes
        for mode in aer_properties["vsd"].keys():
            vsd += utils.gaussian(
                np.log(radius),
                np.log(aer_properties["vsd"][mode]["reff"]),
                aer_properties["vsd"][mode]["rstd"],
                aer_properties["vsd"][mode]["conc"],
            )

        self.radius = radius
        self.vsd = vsd
        self._vsd_to_nsd()

        return self

    def plot(self, show_fig=True):
        """
        Plot Size Distributions of an instance of the (SizeDistributionData): class.

        Example:
            ```python
            # import aprofiles
            import aprofiles as apro
            # get size distribution for urban particles
            sd = apro.size_distribution.SizeDistribData('urban');
            # plot profile
            sd.plot()
            ```

            ![Size distributions for urban particles](../../assets/images/urban_sd.png)
        """
        fig, ax = plt.subplots(1, 1, figsize=(6, 6))

        # plot Volume Size Distribution in 1st axis
        print(self.vsd)
        ax.plot(self.radius, self.vsd, label="VSD")
        ax.set_ylabel("V(r) (µm2.µm-3)")

        # plot Number Size Distribution in 2nd axis
        if "nsd" in self.__dict__:
            # add secondary yaxis
            ax2 = ax.twinx()
            ax2.plot(self.radius, self.nsd, "orange", label="NSD")
            ax2.set_ylabel("N(r) (µm-3.µm-1)")
            # ax2.set_ylim([0,10])

        ax.set_xlabel("Radius (µm)")
        ax.set_xscale("log")
        fig.legend(
            loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax.transAxes
        )
        plt.title(
            f"Size Distribution for {self.aer_type.capitalize().replace('_', ' ')} Particles",
            weight="bold",
        )
        plt.tight_layout()
        if show_fig:
            plt.show()

get_sd()

Returns the Volume and Number Size Distributions arrays from an instance of the (SizeDistributionData): class .

Returns:

Type Description
SizeDistribData

object with additional attributes:

  • radius (numpy.ndarray): radius (µm)
  • vsd (numpy.ndarray): Volume Size Distribution (µm2.µm-3)
  • nsd (numpy.ndarray): Number Size Distribution (µm-3.µm-1)
Source code in aprofiles/size_distribution.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def get_sd(self):
    """
    Returns the Volume and Number Size Distributions arrays from an instance of the (SizeDistributionData): class .

    Returns:
        (SizeDistribData): object with additional attributes:

            - `radius` (numpy.ndarray): radius (µm)
            - `vsd` (numpy.ndarray): Volume Size Distribution (µm2.µm-3)
            - `nsd` (numpy.ndarray): Number Size Distribution (µm-3.µm-1)
    """

    aer_properties = self.aer_properties
    radius = np.arange(1e-2, 2e1, 1e-3)  # radius in µm
    vsd = np.zeros(len(radius))

    # we loop though all the keys defining the different modes
    for mode in aer_properties["vsd"].keys():
        vsd += utils.gaussian(
            np.log(radius),
            np.log(aer_properties["vsd"][mode]["reff"]),
            aer_properties["vsd"][mode]["rstd"],
            aer_properties["vsd"][mode]["conc"],
        )

    self.radius = radius
    self.vsd = vsd
    self._vsd_to_nsd()

    return self

plot(show_fig=True)

Plot Size Distributions of an instance of the (SizeDistributionData): class.

Example 
# import aprofiles
import aprofiles as apro
# get size distribution for urban particles
sd = apro.size_distribution.SizeDistribData('urban');
# plot profile
sd.plot()

Size distributions for urban particles

Source code in aprofiles/size_distribution.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def plot(self, show_fig=True):
    """
    Plot Size Distributions of an instance of the (SizeDistributionData): class.

    Example:
        ```python
        # import aprofiles
        import aprofiles as apro
        # get size distribution for urban particles
        sd = apro.size_distribution.SizeDistribData('urban');
        # plot profile
        sd.plot()
        ```

        ![Size distributions for urban particles](../../assets/images/urban_sd.png)
    """
    fig, ax = plt.subplots(1, 1, figsize=(6, 6))

    # plot Volume Size Distribution in 1st axis
    print(self.vsd)
    ax.plot(self.radius, self.vsd, label="VSD")
    ax.set_ylabel("V(r) (µm2.µm-3)")

    # plot Number Size Distribution in 2nd axis
    if "nsd" in self.__dict__:
        # add secondary yaxis
        ax2 = ax.twinx()
        ax2.plot(self.radius, self.nsd, "orange", label="NSD")
        ax2.set_ylabel("N(r) (µm-3.µm-1)")
        # ax2.set_ylim([0,10])

    ax.set_xlabel("Radius (µm)")
    ax.set_xscale("log")
    fig.legend(
        loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax.transAxes
    )
    plt.title(
        f"Size Distribution for {self.aer_type.capitalize().replace('_', ' ')} Particles",
        weight="bold",
    )
    plt.tight_layout()
    if show_fig:
        plt.show()

Extinction to Mass Coefficient

The EMCData class is used for computing an Extinction to Mass Coefficient for a given aerosol type.

EMCData

Class for computing the Extinction to Mass Coefficient for a given aerosol type. This class calls the get_emc() method.

Attributes:

Name Type Description
`aer_type` {dust, volcanic_ash, biomass_burning, urban}

aerosol type.
`wavelength` int or float

wavelength, in mm.
`method` {mortier_2013, literature}

method to retrieve or compute EMC.
`aer_properties` dict

dictionnary describing the optical and microphophysical properties of the prescribed aerosol (read from aer_properties.json)
Example 
#some imports
import aprofiles as apro
emc_data = EMCData('volcanic_ash', 532.)
emc_data.__dict__.keys()
dict_keys(['aer_type', 'wavelength', 'aer_properties', 'nsd', 'vsd', 'radius', 'qext', 'conv_factor', 'emc'])
print(f'{emc_data.conv_factor:.2e} m {emc_data.emc):.2e} m2.g-1')
6.21e-07 m 0.62 m2.g-1
Source code in aprofiles/emc.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
class EMCData:
    """
    Class for computing the *Extinction to Mass Coefficient* for a given aerosol type.
    This class calls the [`get_emc()`](#aprofiles.emc.EMCData.get_emc) method.

    Attributes:
       `aer_type` ({'dust', 'volcanic_ash', 'biomass_burning', 'urban'}): aerosol type.
       `wavelength` (int or float): wavelength, in mm.
       `method` ({'mortier_2013', 'literature'}): method to retrieve or compute `EMC`.
       `aer_properties` (dict): dictionnary describing the optical and microphophysical properties of the prescribed aerosol (read from *aer_properties.json*)

    Example:
        ```python
        #some imports
        import aprofiles as apro
        emc_data = EMCData('volcanic_ash', 532.)
        emc_data.__dict__.keys()
        dict_keys(['aer_type', 'wavelength', 'aer_properties', 'nsd', 'vsd', 'radius', 'qext', 'conv_factor', 'emc'])
        print(f'{emc_data.conv_factor:.2e} m {emc_data.emc):.2e} m2.g-1')
        6.21e-07 m 0.62 m2.g-1
        ```
    """

    def __init__(self, aer_type: str, wavelength: float, method: str = "mortier_2013"):

        # check parameters type
        if not isinstance(aer_type, str):
            raise TypeError(
                "`aer_type` is expected to be a string. Got {} instead.".format(
                    type(aer_type)
                )
            )
        if not isinstance(wavelength, (int, float)):
            raise TypeError(
                "`wavelength` is expected to be an int or a float. Got {} instead".format(
                    type(wavelength)
                )
            )
        available_methods = ["mortier_2013", "literature"]
        if method not in available_methods:
            raise ValueError(
                "Invalid `method`. AAvailable methods are {}".format(available_methods)
            )

        self.aer_type = aer_type
        self.wavelength = wavelength
        self.method = method

        if self.method == "mortier_2013":
            # read aer_properties.json files
            f = open(Path(Path(__file__).parent,'config','aer_properties.json'))
            aer_properties = json.load(f)
            f.close()
            # check if the aer_type exist in the json file
            if aer_type not in aer_properties.keys():
                raise ValueError(
                    "{} not found in aer_properties.json. `aer_type` must be one of the follwowing: {}".format(
                        aer_type, list(aer_properties.keys())
                    )
                )
            else:
                self.aer_properties = aer_properties[self.aer_type]
                self.get_emc()
        elif self.method == "literature":
            self.emc = 3.33  # CHECK V-PROFILES VALUES. CHECK CODE IN LAPTOP?
            self.conv_factor = -99

    def get_emc(self):
        """
        Calculates the Extinction to Mass Coefficient for a given type of particles, assuming a prescribed size distribution shape (with unknown amplitude), density, and using [Mie theory](https://miepython.readthedocs.io) to calculate the extinction efficiency.

        Returns:
            (EMCData): with additional attributes:

                - `nsd` (1D Array): Number Size Distribution
                - `vsd` (1D Array): Volume Size Distribution
                - `radius` (1D Array): Radius in µm
                - `x` (1D Array): Size parameter (unitless)
                - `conv_factor` (float): Conversion factor in m
                - `emc` (float): Extinction to Mass Coefficient in m².g⁻¹

       !!! note
            For a population of particles, the extinction coefficient $\sigma_{ext}$ (m⁻¹) can be written as follows:

            $$
            \sigma_{ext} = \int_{r_{min}}^{r_{max}} N(r) Q_{ext}(m, r, \lambda) \pi r^2 dr
            $$

            where $Q_{ext}$ is the `extinction efficiency` and $N(r)$ is the `Number Size Distribution` (NSD).

            $Q_{ext}$ varies with the refractive index, $m$, the wavelength, $\lambda$, and can be calculated for spherical particles using Mie theory.

            The total aerosol mass concentration $M_0$ (µg.m⁻³) can be expressed as:

            $$
            M_0 = \int_{r_{min}}^{r_{max}} M(r) dr
            $$

            where $M(r)$ is the mass size distribution (MSD).

            This can be rewritten in terms of NSD and MSD as:

            $$
            M_0 = \int_{r_{min}}^{r_{max}} \\frac{4\pi r^3}{3} \\rho N(r) dr
            $$

            where $\\rho$ is the particle density (kg.m⁻³).

            By normalizing the NSD with respect to the fine mode ($N(r) = N_0 N_1(r)$), we arrive at:

            $$
            M_0 = \sigma_{ext} \\rho \\frac{4}{3} \\frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr}
            $$

            We define the `conversion factor` (in m) as:

            $$
            c_v = \\frac{4}{3} \\frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr}
            $$

            Thus, the equation simplifies to:

            $$
            M_0 = \sigma_{ext} \\rho c_v
            $$

            Finally, the `Extinction to Mass Coefficient` (EMC, also called `mass extinction cross section`, usually in m².g⁻¹) is defined as:

            $$
            EMC = \\frac{\sigma_{ext}}{M_0} = \\frac{1}{\\rho c_v}
            $$

            with $\\rho$ expressed in g.m⁻³.


            <table>
                <thead>
                    <tr>
                        <th>Aerosol Type</th>
                        <th colspan=2>Conversion Factor (µm)</th>
                        <th colspan=2>EMC (m².g⁻¹) </th>
                    </tr>
                    <tr>
                        <th></th>
                        <th>532 nm</th>
                        <th>1064 nm</th>
                        <th>532 nm</th>
                        <th>1064 nm</th>
                    </tr>
                </thead>
                <tbody>
                    <tr>
                        <td>Urban</td>
                        <td>0.31</td>
                        <td>1.92</td>
                        <td>1.86</td>
                        <td>0.31</td>
                    </tr>
                    <tr>
                        <td>Desert dust</td>
                        <td>0.68</td>
                        <td>1.04</td>
                        <td>0.58</td>
                        <td>0.38</td>
                    </tr>
                    <tr>
                        <td>Biomass burning</td>
                        <td>0.26</td>
                        <td>1.28</td>
                        <td>3.30</td>
                        <td>0.68</td>
                    </tr>
                    <tr>
                        <td>Volcanic ash</td>
                        <td>0.62</td>
                        <td>0.56</td>
                        <td>0.62</td>
                        <td>0.68</td>
                    </tr>
                </tbody>
                <caption>Conversion Factors and EMC calculated for the main aerosol types</caption>
            </table>

        """


        def _compute_conv_factor(nsd, qext, radius):
            """Compute Conversion Factor for a given Size Distribution and Efficiency

            Args:
                nsd (1D Array): Number Size Distribution (µm-3.µm-1)
                qext (1D Array): Extinction Efficiency (unitless)
                radius (1D Array): Radius, in µm

            Returns:
                (float): Conversion factor (m)

            """
            # radius resolution
            dr = min(np.diff(radius))

            # integrals
            numerator = [nsd[i] * (radius[i] ** 3) for i in range(len(radius))]
            denominator = [
                nsd[i] * qext[i] * (radius[i] ** 2) for i in range(len(radius))
            ]
            int1 = np.nancumsum(np.asarray(numerator) * dr)[-1]
            int2 = np.nancumsum(np.asarray(denominator) * dr)[-1]

            conv_factor = (4 / 3) * (int1 / int2)

            # conversion form µm to m
            conv_factor = conv_factor * 1e-6
            return conv_factor

        # generate a size distribution for given aer_type
        sd = size_distribution.SizeDistributionData(self.aer_type)

        # calculate efficiency extinction qext
        # size parameter
        # as the radius is in µm and the wavelength is in nm, one must convert the wavelength to µm
        x = [2 * np.pi * r / (self.wavelength * 1e-3) for r in sd.radius]
        # refractive index
        m = complex(
            self.aer_properties["ref_index"]["real"],
            -abs(self.aer_properties["ref_index"]["imag"]),
        )
        # mie calculation
        qext, _qsca, _qback, _g = miepython.mie(m, x)

        # output
        self.nsd = sd.nsd
        self.vsd = sd.vsd
        self.radius = sd.radius
        self.x = x
        self.qext = qext
        self.conv_factor = _compute_conv_factor(sd.nsd, qext, sd.radius)
        self.emc = 1 / (
            self.conv_factor * self.aer_properties["density"] * 1e6
        )  # convert density from g.cm-3 to g.m-3
        return self

    def plot(self, show_fig=True):
        """
        Plot main information of an instance of the [`EMCData`](#aprofiles.emc.EMCData) class.

        Example:
            ```python
            #import aprofiles
            import aprofiles as apro
            #compute emc for biomas burning particles at 532 nm
            emc = apro.emc.EMCData('volcanic_ash', 532.);
            #plot information
            emc.plot()
            ```

            ![Volcanic Ash particles properties used for EMC calculation](../../assets/images/volcanic_ash-emc.png)
        """
        fig, ax = plt.subplots(1, 1, figsize=(6, 6))

        # plot Volume Size Distribution in 1st axis
        ax.plot(self.radius, self.vsd, label="VSD")
        ax.set_ylabel("V(r) (µm2.µm-3)")

        # plot Number Size Distribution in 2nd axis
        if "nsd" in self.__dict__:
            # add secondary yaxis
            ax2 = ax.twinx()
            ax2.plot(self.radius, self.qext, label="Qext", color="gray")
            ax2.set_ylabel("Qext (unitless)")
            # ax2.set_ylim([0,10])

        # add additional information
        plt.text(
            0.975,
            0.85,
            f"$at\ \lambda={self.wavelength:.0f}\ nm$",
            horizontalalignment="right",
            verticalalignment="center",
            transform=ax.transAxes,
        )
        plt.text(
            0.975,
            0.80,
            f"$c_v: {self.conv_factor * 1e6:.2f}\ \mu m$",
            horizontalalignment="right",
            verticalalignment="center",
            transform=ax.transAxes,
        )
        plt.text(
            0.975,
            0.75,
            f"$EMC: {self.emc:.2f}\ m^2/g$",
            horizontalalignment="right",
            verticalalignment="center",
            transform=ax.transAxes,
        )

        ax.set_xlabel("Radius (µm)")
        ax.set_xscale("log")
        fig.legend(
            loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax.transAxes
        )
        plt.title(
            f"{self.aer_type.capitalize().replace('_', ' ')} particles properties for EMC calculation",
            weight="bold",
        )
        plt.tight_layout()
        if show_fig:
            plt.show()

get_emc()

Calculates the Extinction to Mass Coefficient for a given type of particles, assuming a prescribed size distribution shape (with unknown amplitude), density, and using Mie theory to calculate the extinction efficiency.

Returns: (EMCData): with additional attributes:

     - `nsd` (1D Array): Number Size Distribution
     - `vsd` (1D Array): Volume Size Distribution
     - `radius` (1D Array): Radius in µm
     - `x` (1D Array): Size parameter (unitless)
     - `conv_factor` (float): Conversion factor in m
     - `emc` (float): Extinction to Mass Coefficient in m².g⁻¹

Note

For a population of particles, the extinction coefficient \(\sigma_{ext}\) (m⁻¹) can be written as follows:

$$ \sigma_{ext} = \int_{r_{min}}^{r_{max}} N(r) Q_{ext}(m, r, \lambda) \pi r^2 dr $$

where \(Q_{ext}\) is the extinction efficiency and \(N(r)\) is the Number Size Distribution (NSD).

\(Q_{ext}\) varies with the refractive index, \(m\), the wavelength, \(\lambda\), and can be calculated for spherical particles using Mie theory.

The total aerosol mass concentration \(M_0\) (µg.m⁻³) can be expressed as:

$$ M_0 = \int_{r_{min}}^{r_{max}} M(r) dr $$

where \(M(r)\) is the mass size distribution (MSD).

This can be rewritten in terms of NSD and MSD as:

$$ M_0 = \int_{r_{min}}^{r_{max}} \frac{4\pi r^3}{3} \rho N(r) dr $$

where \(\rho\) is the particle density (kg.m⁻³).

By normalizing the NSD with respect to the fine mode (\(N(r) = N_0 N_1(r)\)), we arrive at:

$$ M_0 = \sigma_{ext} \rho \frac{4}{3} \frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr} $$

We define the conversion factor (in m) as:

$$ c_v = \frac{4}{3} \frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr} $$

Thus, the equation simplifies to:

$$ M_0 = \sigma_{ext} \rho c_v $$

Finally, the Extinction to Mass Coefficient (EMC, also called mass extinction cross section, usually in m².g⁻¹) is defined as:

$$ EMC = \frac{\sigma_{ext}}{M_0} = \frac{1}{\rho c_v} $$

with \(\rho\) expressed in g.m⁻³.

Aerosol Type Conversion Factor (µm) EMC (m².g⁻¹)
532 nm 1064 nm 532 nm 1064 nm
Urban 0.31 1.92 1.86 0.31
Desert dust 0.68 1.04 0.58 0.38
Biomass burning 0.26 1.28 3.30 0.68
Volcanic ash 0.62 0.56 0.62 0.68
Conversion Factors and EMC calculated for the main aerosol types

Source code in aprofiles/emc.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
def get_emc(self):
    """
    Calculates the Extinction to Mass Coefficient for a given type of particles, assuming a prescribed size distribution shape (with unknown amplitude), density, and using [Mie theory](https://miepython.readthedocs.io) to calculate the extinction efficiency.

    Returns:
        (EMCData): with additional attributes:

            - `nsd` (1D Array): Number Size Distribution
            - `vsd` (1D Array): Volume Size Distribution
            - `radius` (1D Array): Radius in µm
            - `x` (1D Array): Size parameter (unitless)
            - `conv_factor` (float): Conversion factor in m
            - `emc` (float): Extinction to Mass Coefficient in m².g⁻¹

   !!! note
        For a population of particles, the extinction coefficient $\sigma_{ext}$ (m⁻¹) can be written as follows:

        $$
        \sigma_{ext} = \int_{r_{min}}^{r_{max}} N(r) Q_{ext}(m, r, \lambda) \pi r^2 dr
        $$

        where $Q_{ext}$ is the `extinction efficiency` and $N(r)$ is the `Number Size Distribution` (NSD).

        $Q_{ext}$ varies with the refractive index, $m$, the wavelength, $\lambda$, and can be calculated for spherical particles using Mie theory.

        The total aerosol mass concentration $M_0$ (µg.m⁻³) can be expressed as:

        $$
        M_0 = \int_{r_{min}}^{r_{max}} M(r) dr
        $$

        where $M(r)$ is the mass size distribution (MSD).

        This can be rewritten in terms of NSD and MSD as:

        $$
        M_0 = \int_{r_{min}}^{r_{max}} \\frac{4\pi r^3}{3} \\rho N(r) dr
        $$

        where $\\rho$ is the particle density (kg.m⁻³).

        By normalizing the NSD with respect to the fine mode ($N(r) = N_0 N_1(r)$), we arrive at:

        $$
        M_0 = \sigma_{ext} \\rho \\frac{4}{3} \\frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr}
        $$

        We define the `conversion factor` (in m) as:

        $$
        c_v = \\frac{4}{3} \\frac{\int_{r_{min}}^{r_{max}} N_1(r) r^3 dr}{\int_{r_{min}}^{r_{max}} N_1(r) Q_{ext}(m, r, \lambda) r^2 dr}
        $$

        Thus, the equation simplifies to:

        $$
        M_0 = \sigma_{ext} \\rho c_v
        $$

        Finally, the `Extinction to Mass Coefficient` (EMC, also called `mass extinction cross section`, usually in m².g⁻¹) is defined as:

        $$
        EMC = \\frac{\sigma_{ext}}{M_0} = \\frac{1}{\\rho c_v}
        $$

        with $\\rho$ expressed in g.m⁻³.


        <table>
            <thead>
                <tr>
                    <th>Aerosol Type</th>
                    <th colspan=2>Conversion Factor (µm)</th>
                    <th colspan=2>EMC (m².g⁻¹) </th>
                </tr>
                <tr>
                    <th></th>
                    <th>532 nm</th>
                    <th>1064 nm</th>
                    <th>532 nm</th>
                    <th>1064 nm</th>
                </tr>
            </thead>
            <tbody>
                <tr>
                    <td>Urban</td>
                    <td>0.31</td>
                    <td>1.92</td>
                    <td>1.86</td>
                    <td>0.31</td>
                </tr>
                <tr>
                    <td>Desert dust</td>
                    <td>0.68</td>
                    <td>1.04</td>
                    <td>0.58</td>
                    <td>0.38</td>
                </tr>
                <tr>
                    <td>Biomass burning</td>
                    <td>0.26</td>
                    <td>1.28</td>
                    <td>3.30</td>
                    <td>0.68</td>
                </tr>
                <tr>
                    <td>Volcanic ash</td>
                    <td>0.62</td>
                    <td>0.56</td>
                    <td>0.62</td>
                    <td>0.68</td>
                </tr>
            </tbody>
            <caption>Conversion Factors and EMC calculated for the main aerosol types</caption>
        </table>

    """


    def _compute_conv_factor(nsd, qext, radius):
        """Compute Conversion Factor for a given Size Distribution and Efficiency

        Args:
            nsd (1D Array): Number Size Distribution (µm-3.µm-1)
            qext (1D Array): Extinction Efficiency (unitless)
            radius (1D Array): Radius, in µm

        Returns:
            (float): Conversion factor (m)

        """
        # radius resolution
        dr = min(np.diff(radius))

        # integrals
        numerator = [nsd[i] * (radius[i] ** 3) for i in range(len(radius))]
        denominator = [
            nsd[i] * qext[i] * (radius[i] ** 2) for i in range(len(radius))
        ]
        int1 = np.nancumsum(np.asarray(numerator) * dr)[-1]
        int2 = np.nancumsum(np.asarray(denominator) * dr)[-1]

        conv_factor = (4 / 3) * (int1 / int2)

        # conversion form µm to m
        conv_factor = conv_factor * 1e-6
        return conv_factor

    # generate a size distribution for given aer_type
    sd = size_distribution.SizeDistributionData(self.aer_type)

    # calculate efficiency extinction qext
    # size parameter
    # as the radius is in µm and the wavelength is in nm, one must convert the wavelength to µm
    x = [2 * np.pi * r / (self.wavelength * 1e-3) for r in sd.radius]
    # refractive index
    m = complex(
        self.aer_properties["ref_index"]["real"],
        -abs(self.aer_properties["ref_index"]["imag"]),
    )
    # mie calculation
    qext, _qsca, _qback, _g = miepython.mie(m, x)

    # output
    self.nsd = sd.nsd
    self.vsd = sd.vsd
    self.radius = sd.radius
    self.x = x
    self.qext = qext
    self.conv_factor = _compute_conv_factor(sd.nsd, qext, sd.radius)
    self.emc = 1 / (
        self.conv_factor * self.aer_properties["density"] * 1e6
    )  # convert density from g.cm-3 to g.m-3
    return self

plot(show_fig=True)

Plot main information of an instance of the EMCData class.

Example 
#import aprofiles
import aprofiles as apro
#compute emc for biomas burning particles at 532 nm
emc = apro.emc.EMCData('volcanic_ash', 532.);
#plot information
emc.plot()

Volcanic Ash particles properties used for EMC calculation

Source code in aprofiles/emc.py 
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
def plot(self, show_fig=True):
    """
    Plot main information of an instance of the [`EMCData`](#aprofiles.emc.EMCData) class.

    Example:
        ```python
        #import aprofiles
        import aprofiles as apro
        #compute emc for biomas burning particles at 532 nm
        emc = apro.emc.EMCData('volcanic_ash', 532.);
        #plot information
        emc.plot()
        ```

        ![Volcanic Ash particles properties used for EMC calculation](../../assets/images/volcanic_ash-emc.png)
    """
    fig, ax = plt.subplots(1, 1, figsize=(6, 6))

    # plot Volume Size Distribution in 1st axis
    ax.plot(self.radius, self.vsd, label="VSD")
    ax.set_ylabel("V(r) (µm2.µm-3)")

    # plot Number Size Distribution in 2nd axis
    if "nsd" in self.__dict__:
        # add secondary yaxis
        ax2 = ax.twinx()
        ax2.plot(self.radius, self.qext, label="Qext", color="gray")
        ax2.set_ylabel("Qext (unitless)")
        # ax2.set_ylim([0,10])

    # add additional information
    plt.text(
        0.975,
        0.85,
        f"$at\ \lambda={self.wavelength:.0f}\ nm$",
        horizontalalignment="right",
        verticalalignment="center",
        transform=ax.transAxes,
    )
    plt.text(
        0.975,
        0.80,
        f"$c_v: {self.conv_factor * 1e6:.2f}\ \mu m$",
        horizontalalignment="right",
        verticalalignment="center",
        transform=ax.transAxes,
    )
    plt.text(
        0.975,
        0.75,
        f"$EMC: {self.emc:.2f}\ m^2/g$",
        horizontalalignment="right",
        verticalalignment="center",
        transform=ax.transAxes,
    )

    ax.set_xlabel("Radius (µm)")
    ax.set_xscale("log")
    fig.legend(
        loc="upper right", bbox_to_anchor=(1, 1), bbox_transform=ax.transAxes
    )
    plt.title(
        f"{self.aer_type.capitalize().replace('_', ' ')} particles properties for EMC calculation",
        weight="bold",
    )
    plt.tight_layout()
    if show_fig:
        plt.show()

  1. Dubovik, O., Holben, B., Eck, T. F., Smirnov, A., Kaufman, Y. J., King, M. D., ... & Slutsker, I. (2002). Variability of absorption and optical properties of key aerosol types observed in worldwide locations. Journal of the Atmospheric Sciences, 59(3), 590-608. ↩

  2. Mortier, A., Goloub, P., Podvin, T., Deroo, C., Chaikovsky, A., Ajtai, N., ... & Derimian, Y. (2013). Detection and characterization of volcanic ash plumes over Lille during the Eyjafjallajökull eruption. Atmospheric Chemistry and Physics, 13(7), 3705-3720. ↩