wradlib Phase Processing - System PhiDP - ZPHI-Method

Overview

Within this notebook, we will cover:

Reading sweep data into xarray based Dataset
Retrieval of system PhiDP
ZPHI Phase processing
Attenuation correction using specific Attenuation

Prerequisites

Concepts	Importance	Notes
Xarray Basics	Helpful	Basic Dataset/DataArray
Matplotlib Basics	Helpful	Basic Plotting
Intro to Cartopy	Helpful	Projections

Time to learn: 10 minutes

import datetime as dt
import glob
import os
import sys
import warnings

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from scipy.integrate import cumulative_trapezoid

import wradlib as wrl
import xradar as xd

Import data

As a quick example to show the algorithm, we use a file from Down Under. For the further processing we use x-band data from BoXPol research radar at the University of Bonn, Germany.

boxpol = "data/hdf5/boxpol/2014-11-16--03_45_00,00.mvol"
terrey = "data/hdf5/terrey_39.h5"

swp0 = xr.open_dataset(terrey, engine="odim")
swp0 = swp0.pipe(wrl.georef.georeference_dataset).compute()

(we use .compute() here to load the data out of a dask array, otherwise some functions below will complain)

Pre-Processing

System PHIDP aka Phase Offset

The following function returns phase offset as well as start and stop ranges of the region of interest (first precipitating bins).

def phase_offset(phioff, method=None, rng=3000.0, npix=None, **kwargs):
    """Calculate Phase offset.

    Parameter
    ---------
    phioff : xarray.DataArray
        differential phase DataArray

    Keyword Arguments
    -----------------
    method : str
        aggregation method, defaults to 'median'
    rng : float
        range in m to calculate system phase offset

    Return
    ------
    phidp_offset : xarray.Dataset
        Dataset with PhiDP offset and start/stop ranges
    """
    range_step = np.diff(phioff.range)[0]
    nprec = int(rng / range_step)
    if not nprec % 2:
        nprec += 1

    if npix is None:
        npix = nprec // 2 + 1

    # create binary array
    phib = xr.where(np.isnan(phioff), 0, 1)

    # take nprec range bins and calculate sum
    phib_sum = phib.rolling(range=nprec, **kwargs).sum(skipna=True)

    # find at least N pixels in
    # phib_sum_N = phib_sum.where(phib_sum >= npix)
    phib_sum_N = xr.where(phib_sum <= npix, phib_sum, npix)

    # get start range of first N consecutive precip bins
    start_range = (
        phib_sum_N.idxmax(dim="range") - nprec // 2 * np.diff(phib_sum.range)[0]
    )
    start_range = xr.where(start_range < 0, 0, start_range)

    # get stop range
    stop_range = start_range + rng
    # get phase values in specified range
    off = phioff.where(
        (phioff.range >= start_range) & (phioff.range <= stop_range), drop=False
    )
    # calculate nan median over range
    if method is None:
        method = "median"
    func = getattr(off, method)
    off_func = func(dim="range", skipna=True)

    return xr.Dataset(
        dict(
            PHIDP_OFFSET=off_func,
            start_range=start_range,
            stop_range=stop_range,
            phib_sum=phib_sum,
            phib=phib,
        )
    )

Example Showcase

dr_m = swp0.range.diff("range").median()
swp_msk = swp0.where((swp0.DBZH >= 0.0))
swp_msk = swp_msk.where(swp_msk.RHOHV > 0.8)
swp_msk = swp_msk.where(swp_msk.range > dr_m * 5)

phi_masked = swp_msk.PHIDP.copy()
off = phase_offset(
    phi_masked, method="median", rng=2000.0, npix=7, center=True, min_periods=4
)
phioff = off.PHIDP_OFFSET.median(dim="azimuth", skipna=True)

fig = plt.figure(figsize=(16, 7))
ax1 = plt.subplot(111, projection="polar")
# set the lable go clockwise and start from the top
ax1.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
ax1.plot(theta, off.PHIDP_OFFSET, color="b", linewidth=3)

ax1.plot(theta, np.ones_like(theta) * phioff.values, color="r", lw=2)
ti = off.time.values.astype("M8[s]")[0]
om = phioff.values
tx = ax1.set_title(f"{ti}\n" + r"$\phi_{DP}-Offset$ " + f"{om:.1f} (deg)")
tx.set_y(1.1)
xticks = ax1.set_xticks(np.pi / 180.0 * np.linspace(0, 360, 36, endpoint=False))
ax1.set_ylim(50, 150)

(50.0, 150.0)

../../_images/3b89a045897c4ead0d687b5ebc6d50a362220049670e7e4dabdb729dc548422d.png

fig = plt.figure(figsize=(18, 5))
swp_msk.DBZH.plot(x="azimuth")
off.start_range.plot(c="b", lw=2)
off.stop_range.plot(c="r", lw=2)
plt.gca().set_ylim(0, 25000)

(0.0, 25000.0)

../../_images/ea1b15c9399b008ec8d31f0a8d5e677d54349ea8a9a21c7ad2d348bad385f857.png

Process BoXPol data

swp = xr.open_dataset(boxpol, engine="gamic")
swp = swp.pipe(wrl.georef.georeference_dataset)

vol = wrl.io.open_gamic_dataset(boxpol)

/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/io/hdf.py:139: FutureWarning: `open_gamic_dataset` functionality has been moved to `xradar`-package and will be removed in 2.0. Use `open_gamic_datatree` from `xradar`-package.
  return open_radar_dataset(filename_or_obj, engine=GamicBackendEntrypoint, **kwargs)

display(vol)

<wradlib.RadarVolume>
Dimension(s): (sweep: 1)
Elevation(s): (1.5)

swp = vol[0].copy()
swp = swp.pipe(wrl.georef.georeference_dataset)

Create Plot

fig = plt.figure(figsize=(13, 5))

ax1 = fig.add_subplot(121)
im1 = swp.PHIDP.where(swp.RHOHV > 0.8).plot(x="x", y="y", ax=ax1, cmap="turbo")
t = plt.title(r"Uncorrected $\phi_{DP}$")
t.set_y(1.1)

ax2 = fig.add_subplot(122)
im2 = swp.DBZH.where(swp.RHOHV > 0.8).plot(
    x="x", y="y", ax=ax2, cmap="turbo", vmin=-10, vmax=50
)
t = plt.title(r"Uncorrected $Z_{H}$")
t.set_y(1.1)
fig.suptitle(swp.time.values, fontsize=14)
fig.subplots_adjust(wspace=0.25)

../../_images/07d2ebed0d608163db53fa08ab6a4228a76eccb325b53dbdcce41aeab2737578.png

Apply reasonable masking

dr_m = swp.range.diff("range").median()
swp_msk = swp.where((swp.DBZH >= 0.0))
swp_msk = swp_msk.where(swp_msk.RHOHV > 0.8)
swp_msk = swp_msk.where(swp_msk.range > dr_m * 2)


phi_masked = swp_msk.PHIDP.copy()
off = phase_offset(
    phi_masked, method="median", rng=2000.0, npix=7, center=True, min_periods=2
)
phioff = off.PHIDP_OFFSET.median(dim="azimuth", skipna=True)

Plot phase offset distribution

fig = plt.figure(figsize=(16, 7))
ax1 = plt.subplot(111, projection="polar")
# set the lable go clockwise and start from the top
ax1.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
ax1.plot(theta, off.PHIDP_OFFSET, color="b", linewidth=3)

ax1.plot(theta, np.ones_like(theta) * phioff.values, color="r", lw=2)
ti = off.time.values.astype("M8[s]")
om = phioff.values
tx = ax1.set_title(f"{ti}\n" + r"$\phi_{DP}-Offset$ " + f"{om:.1f} (deg)")
tx.set_y(1.1)
xticks = ax1.set_xticks(np.pi / 180.0 * np.linspace(0, 360, 36, endpoint=False))
ax1.set_ylim(-120, -70)

(-120.0, -70.0)

../../_images/534011c8c90bf35d249d06576206274a0f0a90b3a1dff9c1bfa576a6ff5eb0aa.png

fig = plt.figure(figsize=(18, 5))
swp_msk.DBZH.plot(x="azimuth")
off.start_range.plot(c="b", lw=2)
off.stop_range.plot(c="r", lw=2)
plt.gca().set_ylim(0, 10000)

(0.0, 10000.0)

../../_images/be96f1cf55928e623bd79a303b9c5994623e78771aeff15033b35fa0861fbee4.png

Pleaser refer to the ZPHI-Method section at the bottom of this notebook for references and equations.

Retrieving \(\Delta \phi_{DP}\)

We will use the simple method of finding the first and the last non NAN values per ray from \(\phi_{DP}^{corr}\).

This is the most simple and probably not very robust method.

def phase_zphi(phi, rng=1000.0, **kwargs):
    range_step = np.diff(phi.range)[0]

    nprec = int(rng / range_step)

    if nprec % 2:
        nprec += 1

    # create binary array
    phib = xr.where(np.isnan(phi), 0, 1)

    # take nprec range bins and calculate sum
    phib_sum = phib.rolling(range=nprec, **kwargs).sum(skipna=True)

    offset = nprec // 2 * np.diff(phib_sum.range)[0]
    offset_idx = nprec // 2

    start_range = phib_sum.idxmax(dim="range") - offset
    start_range_idx = phib_sum.argmax(dim="range") - offset_idx

    stop_range = phib_sum[:, ::-1].idxmax(dim="range") - offset
    stop_range_idx = (
        len(phib_sum.range) - (phib_sum[:, ::-1].argmax(dim="range") - offset_idx) - 2
    )

    # get phase values in specified range
    first = phi.where(
        (phi.range >= start_range) & (phi.range <= start_range + rng), drop=True
    ).quantile(0.15, dim="range", skipna=True)
    last = phi.where(
        (phi.range >= stop_range - rng) & (phi.range <= stop_range), drop=True
    ).quantile(0.95, dim="range", skipna=True)

    return xr.Dataset(
        dict(
            phib=phib_sum,
            offset=offset,
            offset_idx=offset_idx,
            start_range=start_range,
            stop_range=stop_range,
            first=first.drop("quantile"),
            first_idx=start_range_idx,
            last=last.drop("quantile"),
            last_idx=stop_range_idx,
        )
    )

Apply extraction of phase parameters.

cphase = phase_zphi(swp_msk.PHIDP, rng=2000.0, center=True, min_periods=7)

Apply azimuthal averaging.

cphase = (
    cphase.pad(pad_width={"azimuth": 2}, mode="wrap")
    .rolling(azimuth=5, center=True)
    .median(skipna=True)
    .isel(azimuth=slice(2, -2))
)

\(\Delta \phi_{DP}\) - Polar Plots

This visualizes first and last indizes including \(\Delta \phi_{DP}\).

dphi = cphase.last - cphase.first
dphi = dphi.where(dphi >= 0).fillna(0)

fig = plt.figure(figsize=(20, 9))
ax1 = plt.subplot(131, projection="polar")
ax2 = plt.subplot(132, projection="polar")
ax3 = plt.subplot(133, projection="polar")
# set the lable go clockwise and start from the top
ax1.set_theta_zero_location("N")
ax2.set_theta_zero_location("N")
ax3.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
ax2.set_theta_direction(-1)
ax3.set_theta_direction(-1)
theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
ax1.plot(theta, cphase.start_range, color="b", linewidth=2)
ax1.plot(theta, cphase.stop_range, color="r", linewidth=2)
_ = ax1.set_title("Start/Stop Range")

ax2.plot(theta, cphase.first, color="b", linewidth=2)
ax2.plot(theta, cphase.last, color="r", linewidth=2)
_ = ax2.set_title("Start/Stop PHIDP")
ax2.set_ylim(-110, -40)

ax3.plot(theta, dphi, color="g", linewidth=3)
# ax3.plot(theta, dphi_old, color="k", linewidth=1)
_ = ax3.set_title("Delta PHIDP")

../../_images/c70294343f5db1e243ad82372b3ccc34df0903ec15be3848530bc14364f06b22.png

Calculating Reflectivity Integrals/Sums

\[za(r) = \left[Z_a(r) \right ]^b\]

\[iza(r,r2) = 0.46 \cdot b \cdot \int_{r}^{r2} \left [Z_a(s) \right ]^b ds\]

We do not restrict (mask) the reflectivities for now, but switch between DBTH and DBZH to see the difference.

zhraw = swp.DBZH.where(
    (swp.range > cphase.start_range) & (swp.range < cphase.stop_range)
)
zhraw.plot(x="x", y="y", cmap="turbo", vmin=0, vmax=100)

<matplotlib.collections.QuadMesh at 0x7f58841fb310>

../../_images/f678564d123a96fc9d90a765dc04a299556c78ef495ab61065fc6e03112679f7.png

# calculate linear reflectivity and ^b
zax = zhraw.pipe(wrl.trafo.idecibel).fillna(0)
za = zax**bx
# set masked to zero for integration
za_zero = za.fillna(0)

Calculate cumulative integral, and subtract from maximum. That way we have the cumulative sum for every bin until the end of the ray.

def cumulative_trapezoid_xarray(da, dim, initial=0):
    """Intgration with the scipy.integrate.cumtrapz.

    Parameter
    ---------
    da : xarray.DataArray
        array with differential phase data
    dim : int
        size of window in range dimension

    Keyword Arguments
    -----------------
    initial : float
        minimum number of valid bins

    Return
    ------
    kdp : xarray.DataArray
        DataArray with specific differential phase values
    """
    x = da[dim]
    dx = x.diff(dim).median(dim).values
    if x.attrs["units"] == "meters":
        dx /= 1000.0
    return xr.apply_ufunc(
        cumulative_trapezoid,
        da,
        input_core_dims=[[dim]],
        output_core_dims=[[dim]],
        dask="parallelized",
        kwargs=dict(axis=da.get_axis_num(dim), initial=initial, dx=dx),
        dask_gufunc_kwargs=dict(allow_rechunk=True),
    )

iza_x = 0.46 * bx * za_zero.pipe(cumulative_trapezoid_xarray, "range", initial=0)
iza = iza_x.max("range") - iza_x

Calculating Attenuation \(A_{H}\) for whole domain

\[A_{H}(r) = \frac{\left [Z_a(r) \right ]^b \cdot f(\Delta \phi_{DP})}{0.46b \int_{r1}^{r2} \left [Z_a(s) \right ]^b ds + f(\Delta \phi_{DP}) \cdot 0.46b \int_{r}^{r2} \left [Z_a(s) \right ]^b ds}\]

We can reduce the number of operations by rearranging the equation like this:

\[A_{H}(r) = \frac{\left [Z_a(r) \right ]^b}{\frac{0.46b \int_{r1}^{r2} \left [Z_a(s) \right ]^b ds}{f(\Delta \phi_{DP})} + 0.46b \int_{r}^{r2} \left [Z_a(s) \right ]^b ds}\]

iza_fdphi = iza / fdphi
idx = cphase.first_idx.astype(int)
iza_first = iza_fdphi[:, idx]
ah = za / (iza_first + iza)

Give it a name!

ah.name = "AH"
ah.attrs["short_name"] = "specific_attenuation_h"
ah.attrs["long_name"] = "Specific attenuation H"
ah.attrs["units"] = "dB/km"

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)
ticks_ah = np.arange(0, 5, 0.2)
im = ah.plot(x="x", y="y", ax=ax, cmap="turbo", levels=np.arange(0, 0.5, 0.025))

../../_images/61be8740d7b9f5f768077ce4c23e72a1712ba92a3e30acd655367cf9ab427d31.png

Calculate \(\phi_{DP}^{cal}(r, \alpha)\) for whole domain

\[\phi_{DP}^{cal}(r_i, \alpha) = 2 \cdot \int_{r1}^{r2} \frac{A_H(s; \alpha)}{\alpha}ds\]

phical = 2 * (ah / alphax).pipe(cumulative_trapezoid_xarray, "range", initial=0)
phical.name = "PHICAL"
phical.attrs = xd.model.sweep_vars_mapping["PHIDP"]

phical.where(swp_msk.PHIDP).plot(x="x", y="y", vmin=0, vmax=50, cmap="turbo")

<matplotlib.collections.QuadMesh at 0x7f586ff7d050>

../../_images/4b099c304a1d6977cb9835357b15b6a68f6b6adef2034678eb81a05eb23c0b5a.png

Apply attenuation correction

print(alphax)

0.28

zhraw = swp.DBZH.copy()
zdrraw = swp.ZDR.copy()

with xr.set_options(keep_attrs=True):
    zhcorr = zhraw + alphax * (phical)
    zdiff = zhcorr - zhraw
    zdrcorr = zdrraw + betax * (phical)
    zdrdiff = zdrcorr - zdrraw

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    nrows=2,
    ncols=2,
    figsize=(15, 12),
    sharex=True,
    sharey=True,
    squeeze=True,
    constrained_layout=True,
)

scantime = zhraw.time.values.astype("<M8[s]")
t = fig.suptitle(scantime, fontsize=14)
t.set_y(1.05)

zhraw.plot(x="x", y="y", ax=ax1, cmap="turbo", levels=np.arange(0, 40, 2))
ax1.set_title(r"Uncorrected $Z_{H}$", fontsize=16)
zhcorr.plot(x="x", y="y", ax=ax2, cmap="turbo", levels=np.arange(0, 40, 2))
ax2.set_title(r"Corrected $Z_{H}$", fontsize=16)

zdrraw.plot(x="x", y="y", ax=ax3, cmap="turbo", levels=np.arange(-0.5, 3, 0.1))
ax3.set_title(r"Uncorrected $Z_{DR}$", fontsize=16)
zdrcorr.plot(x="x", y="y", ax=ax4, cmap="turbo", levels=np.arange(-0.5, 3, 0.1))
ax4.set_title(r"Corrected $Z_{DR}$", fontsize=16)

Text(0.5, 1.0, 'Corrected $Z_{DR}$')

../../_images/67500b39a0c7d980ad8f0b9b810616bdbe86a87e5853bbd7b24fc355b338cc70.png

\(K_{DP}\) from \(A_H\) vs. \(K_{DP}\) from \(\phi_{DP}\)

\(K_{DP} = \frac{A_H}{\alpha}\)
\(K_{DP} = \frac{1}{2}\frac{\mathrm{d}\phi_{DP}}{\mathrm{d}r}\)

What are the benefits of \(K_{DP}(A_H)\)?

no noise artefacts
no \(\delta\)
no negative \(K_{DP}\)
no spatial degradation

def kdp_from_phidp(da, winlen, min_periods=2):
    """Derive KDP from PHIDP (based on convolution filter).

    Parameter
    ---------
    da : xarray.DataArray
        array with differential phase data
    winlen : int
        size of window in range dimension

    Keyword Arguments
    -----------------
    min_periods : int
        minimum number of valid bins

    Return
    ------
    kdp : xarray.DataArray
        DataArray with specific differential phase values
    """
    dr = da.range.diff("range").median("range").values / 1000.0
    print("range res [km]:", dr)
    print("processing window [km]:", dr * winlen)
    return xr.apply_ufunc(
        wrl.dp.kdp_from_phidp,
        da,
        input_core_dims=[["range"]],
        output_core_dims=[["range"]],
        dask="parallelized",
        kwargs=dict(winlen=winlen, dr=dr, min_periods=min_periods),
        dask_gufunc_kwargs=dict(allow_rechunk=True),
    )


def kdp_phidp_vulpiani(da, winlen, min_periods=2):
    """Derive KDP from PHIDP (based on Vulpiani).

    ParameterRHOHV_NC
    ---------
    da : xarray.DataArray
        array with differential phase data
    winlen : int
        size of window in range dimension

    Keyword Arguments
    -----------------
    min_periods : int
        minimum number of valid bins

    Return
    ------
    kdp : xarray.DataArray
        DataArray with specific differential phase values
    """
    dr = da.range.diff("range").median("range").values / 1000.0
    print("range res [km]:", dr)
    print("processing window [km]:", dr * winlen)
    return xr.apply_ufunc(
        wrl.dp.process_raw_phidp_vulpiani,
        da,
        input_core_dims=[["range"]],
        output_core_dims=[["range"], ["range"]],
        dask="parallelized",
        kwargs=dict(winlen=winlen, dr=dr, min_periods=min_periods),
        dask_gufunc_kwargs=dict(allow_rechunk=True),
    )

%%time
kdp1 = kdp_from_phidp(swp_msk.PHIDP, winlen=31, min_periods=11)
kdp1.attrs = xd.model.sweep_vars_mapping["KDP"]

kdp2 = kdp_phidp_vulpiani(swp.PHIDP, winlen=71, min_periods=21)[1]
kdp2.attrs = xd.model.sweep_vars_mapping["KDP"]

kdp3 = xr.zeros_like(kdp1)
kdp3.attrs = xd.model.sweep_vars_mapping["KDP"]
kdp3.data = ah / alphax

range res [km]: 0.1
processing window [km]: 3.1
range res [km]: 0.1
processing window [km]: 7.1000000000000005
CPU times: user 383 ms, sys: 84.9 ms, total: 468 ms
Wall time: 475 ms

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    nrows=2, ncols=2, figsize=(12, 10), constrained_layout=True
)

swp.KDP.plot(
    x="x",
    y="y",
    ax=ax1,
    cmap="turbo",
    levels=np.arange(-0.5, 1, 0.1),
    cbar_kwargs=dict(shrink=0.64),
)
ax1.set_title(r"$K_{DP}$ - Signalprocessor", fontsize=16)
ax1.set_aspect("equal")
kdp1.plot(
    x="x",
    y="y",
    ax=ax2,
    cmap="turbo",
    levels=np.arange(-0.5, 1, 0.1),
    cbar_kwargs=dict(shrink=0.64),
)
ax2.set_title(r"$K_{DP}$ - Simple Derivative", fontsize=16)
ax2.set_aspect("equal")
kdp2.plot(
    x="x",
    y="y",
    ax=ax3,
    cmap="turbo",
    levels=np.arange(-0.5, 1, 0.1),
    cbar_kwargs=dict(shrink=0.64),
)
ax3.set_title(r"$K_{DP}$ - Vulpiani", fontsize=16)
ax3.set_aspect("equal")
kdp3.plot(
    x="x",
    y="y",
    ax=ax4,
    cmap="turbo",
    levels=np.arange(-0.5, 1, 0.1),
    cbar_kwargs=dict(shrink=0.64),
)
ax4.set_title(r"$K_{DP}$ - spec. Attenuation/ZPHI", fontsize=16)
ax4.set_aspect("equal")

../../_images/af9ab5654e68d79750f23664b23d65914704e68a73d9a823f57bfc7ce0d021d1.png

Summary

We’ve just learned how to derive System Phase Offset and Specific Attenuation AH using the ZPHI-Method. Different KDP derivation methods have been compared.

Resources and references

xarray
dask
wradlib xarray backends
OPERA ODIM_H5
WMO JET-OWR
Testud, J., Le Bouar, E., Obligis, E., & Ali-Mehenni, M. (2000). The Rain Profiling Algorithm Applied to Polarimetric Weather Radar, Journal of Atmospheric and Oceanic Technology, 17(3), 332-356. Retrieved Nov 24, 2021, from https://journals.ametsoc.org/view/journals/atot/17/3/1520-0426_2000_017_0332_trpaat_2_0_co_2.xml
Diederich, M., Ryzhkov, A., Simmer, C., Zhang, P., & Trömel, S. (2015). Use of Specific Attenuation for Rainfall Measurement at X-Band Radar Wavelengths.: Part I: Radar Calibration and Partial Beam Blockage Estimation. Journal of Hydrometeorology, 16(2), 487–502. http://www.jstor.org/stable/24914953

ZPHI-Method

see Testud et.al. (chapter 4. p. 339ff.), Diederich et.al. (chapter 3. p. 492 ff).

There is a equational difference in the two papers, which can be solved like this:

\(\begin{equation} f\Delta\phi_{DP} = 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1 \tag{1} \end{equation}\)

\(\begin{equation} C(b, PIA) = \exp[{0.23 \cdot b \cdot (PIA)}] - 1 \tag{2} \end{equation}\)

with

\(\begin{equation} PIA = \alpha \cdot \Delta\phi_{DP} \tag{3} \end{equation}\)

\(\begin{equation} C(b, PIA) = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{4} \end{equation}\)

Both expressions are used equivalently:

\(\begin{equation} 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{5} \end{equation}\)

Using logarithmic identities:

\(\begin{equation} \ln {u^r} = r \cdot \ln {u} \tag{6a} \end{equation}\)

\(\begin{equation} \exp {\ln x} = x \tag{6b} \end{equation}\)

the left hand side can be further expressed as:

\(\begin{equation} \exp [\ln {10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7a} \end{equation}\)

\(\begin{equation} \exp[0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP} \cdot \ln {10}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7b} \end{equation}\)

leading to equality

\(\begin{equation} \exp[0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7c} \end{equation}\)