Allow xarray.DataArray as input in pycomlink.processing.k_R_relation.a_b()

[maxmargraf]

Until now a_b() allows int, float or np.array.
In pycomlink.processing.nearby_rain_retrival.nearby_rainfall_retrival() the use of xarray.DataArrays is allowed addtionally to int, float and np.array but it uses a for loop to get the a and b values from a_b().

Hence, updating a_b() to digest individual or iterate-able int, float, np.arrray or xr.DataArrays would be good.

[cchwala]

Some comments on this:

This already works

pycml.processing.k_R_relation.a_b(
    f_GHz=ds_cmls.frequency/1e9,
    pol='V',
)

But using any type of array for pol does not work because we have the if-clause to distinguish V and H. So this does not work

pycml.processing.k_R_relation.a_b(
    f_GHz=ds_cmls.frequency/1e9,
    pol=ds_cmls.polarization,
)

An additional limitation of the current code is that the above example that works, does not return an xr.DataArray but a np.ndarray. It would be best if a xr.DataArray is returned in case of passing in one. Maybe the solution would be to always cast all input to an xr.DataArray and then use .where to calculate and set a and b for the two polarization separately.

[cchwala]

To be able to do the calculation of a and b in a vectorized form for the two polarizations, we can just build separate interpolators, e.g. instead of

        if pol == "V" or pol == "v" or pol == 'Vertical' or pol == "vertical":
            f_a = interp1d(ITU_table[0, :], ITU_table[2, :], kind="cubic")
            f_b = interp1d(ITU_table[0, :], ITU_table[4, :], kind="cubic")
        elif pol == "H" or pol == "h" or pol == "Horizontal" or pol == "horizontal":
            f_a = interp1d(ITU_table[0, :], ITU_table[1, :], kind="cubic")
            f_b = interp1d(ITU_table[0, :], ITU_table[3, :], kind="cubic")

do

f_a_v = interp1d(ITU_table[0, :], ITU_table[2, :], kind="cubic")
f_b_v = interp1d(ITU_table[0, :], ITU_table[4, :], kind="cubic")
f_a_h = interp1d(ITU_table[0, :], ITU_table[1, :], kind="cubic")
f_b_h = interp1d(ITU_table[0, :], ITU_table[3, :], kind="cubic")

and then use this with xarray or numpy indexing for the two cases of polarization.

Maybe something like

a_v = f_a_v(f_GHz.where(pol=='V')) # have to add all other options to specify pol
a_h = f_a_h(f_GHz.where(pol=='H'))

a = a_v.where(pol='V', a_v, a_h) # not very readable...

[cchwala]

closed via #141