<commit_before># -*- coding: UTF-8 -*-
#-------------------------------------------------------------------------------
# Name:         atten
# Purpose:
#
# Authors:      <NAME>, <NAME> and <NAME>
#
# Created:      26.10.2011
# Copyright:    (c) <NAME>, <NAME> and <NAME> 2011
# Licence:      The MIT License
#-------------------------------------------------------------------------------
#!/usr/bin/env python

"""
Attenuation Correction
^^^^^^^^^^^^^^^^^^^^^^

.. autosummary::
   :nosignatures:
   :toctree: generated/

    correctAttenuationHB

"""
import numpy as np
import os, sys
import math
import logging

logger = logging.getLogger('attcorr')


class AttenuationOverflowError(Exception):
    pass


def correctAttenuationHB(gateset, coefficients=None, mode='', thrs=59.0):
    """Gate-by-Gate attenuation correction according to Hitschfeld & Bordan
    [Hitschfeld1954]_



    Parameters
    ----------
    gateset : array
        multidimensional array. The range gates (over which iteration has to
        be performed) are supposed to vary along
        the *last* dimension so, e.g., for a set of `l` radar images stored in
        polar form with `m` azimuths and `n` range-bins the input array's
        shape can be either (l,m,n) or (m,l,n)
        data havs to be provided in decibel representation of reflectivity (dBZ)
    coefficients : dictionary
        correction coefficients
        a: proportionality factor of the k-Z relation ( :math:`k=a*Z^{b}` )
        b: exponent of the k-Z relation
        l: length of a range gate.
        if set to None the following default dictionary will be used
        {'a':1.67e-4, 'b':0.70, 'l':1.0}

    mode : string
        controls how the function reacts, if the sum of signal and attenuation
        exceeds the
        threshold ``thrs``
        Possible values:
        'warn' : emit a warning through the module's logger but continue
        execution
        'zero' : set offending gates to 0.0
        'nan' : set offending gates to nan
        Any other mode will raise an Exception

    thrs : float
        threshold, for the sum of attenuation and signal, which is deemed
        unplausible.

    Returns
    -------
    k : array
        Array with the same shape as ``gateset`` containing the calculated
        attenuation for each range gate

    Raises
    ------
    AttenuationOverflowError
        Exception, if attenuation exceeds ``thrs`` and no handling ``mode`` is
        set.

    References
    ----------

    .. [Hitschfeld1954] <NAME>. & <NAME>., 1954.
        Errors Inherent in the Radar Measurement of Rainfall at Attenuating
        Wavelengths. Journal of the Atmospheric Sciences, 11(1), p.58-67.
        DOI: 10.1175/1520-0469(1954)011<0058:EIITRM>2.0.CO;2

    .. comment
        _[1] Krämer2008 - <NAME> 2008: Quantitative Radardatenaufbereitung
        für die Niederschlagsvorhersage und die Siedlungsentwässerung,
        Mitteilungen Institut für Wasserwirtschaft, Hydrologie und
        Landwirtschaftlichen Wasserbau
        <NAME> Universität Hannover, Heft 92, ISSN 0343-8090

    """
    if coefficients is None:
        _coefficients = {'a':1.67e-4, 'b':0.7, 'l':1.0}
    else:
        _coefficients = coefficients

    a = _coefficients['a']
    b = _coefficients['b']
    l = _coefficients['l']


    k = np.empty(gateset.shape)
    k[...,0] = 0.
    ksum = 0.

    # multidimensional version
    # assumes that iteration is only along the last dimension (i.e. range gates)
    # all other dimensions are calculated simultaneously to gain some speed
    for gate in range(gateset.shape[-1]-1):
        # calculate k in dB/km from k-Z relation
        # c.f. Krämer2008(p. 147)
        kn = a * (10.0**((gateset[...,gate] + ksum)/10.0))**b  * 2.0 * l
        #kn = 10**(log10(a)+0.1*bin*b)
        #dBkn = 10*math.log10(a) + (bin+ksum)*b + 10*math.log10(2*l)
        ksum += kn

        k[...,gate+1] = ksum
        # stop-criterion, if corrected reflectivity is larger than 59 dBZ
        overflow = (gateset[...,gate] + ksum) > thrs
        if np.any(overflow):
            if mode == 'warn':
                logger.warning('dB-sum over threshold (%3.1f)'%thrs)
            if mode == 'nan':
                k[gate,overflow] = np.nan
            if mode == 'zero':
                k[gate,overflow] = 0.0
            else:
                raise AttenuationOverflowError

    return k


if __name__ == '__main__':
    print 'wradlib: Calling module <atten> as main...'<commit_msg>added attenuation correction according to Stefan Krämer as multidimensional version<commit_after># -*- coding: UTF-8 -*-
#-------------------------------------------------------------------------------
# Name:         atten
# Purpose:
#
# Authors:      <NAME>, <NAME> and <NAME>
#
# Created:      26.10.2011
# Copyright:    (c) <NAME>, <NAME> and <NAME> 2011
# Licence:      The MIT License
#-------------------------------------------------------------------------------
#!/usr/bin/env python

"""
Attenuation Correction
^^^^^^^^^^^^^^^^^^^^^^

.. autosummary::
   :nosignatures:
   :toctree: generated/

    correctAttenuationHB
    correctAttenuationKraemer

"""
import numpy as np
import os, sys
import math
import logging

logger = logging.getLogger('attcorr')


class AttenuationOverflowError(Exception):
    pass


def correctAttenuationHB(gateset, coefficients=None, mode='', thrs=59.0):
    """Gate-by-Gate attenuation correction according to Hitschfeld & Bordan
    [Hitschfeld1954]_



    Parameters
    ----------
    gateset : array
        multidimensional array. The range gates (over which iteration has to
        be performed) are supposed to vary along
        the *last* dimension so, e.g., for a set of `l` radar images stored in
        polar form with `m` azimuths and `n` range-bins the input array's
        shape can be either (l,m,n) or (m,l,n)
        data havs to be provided in decibel representation of reflectivity (dBZ)
    coefficients : dictionary
        correction coefficients
        a: proportionality factor of the k-Z relation ( :math:`k=a*Z^{b}` )
        b: exponent of the k-Z relation
        l: length of a range gate.
        if set to None the following default dictionary will be used
        {'a':1.67e-4, 'b':0.70, 'l':1.0}

    mode : string
        controls how the function reacts, if the sum of signal and attenuation
        exceeds the
        threshold ``thrs``
        Possible values:
        'warn' : emit a warning through the module's logger but continue
        execution
        'zero' : set offending gates to 0.0
        'nan' : set offending gates to nan
        Any other mode will raise an Exception

    thrs : float
        threshold, for the sum of attenuation and signal, which is deemed
        unplausible.

    Returns
    -------
    k : array
        Array with the same shape as ``gateset`` containing the calculated
        attenuation for each range gate.

    Raises
    ------
    AttenuationOverflowError
        Exception, if attenuation exceeds ``thrs`` and no handling ``mode`` is
        set.

    References
    ----------

    .. [Hitschfeld1954] <NAME>. & <NAME>., 1954.
        Errors Inherent in the Radar Measurement of Rainfall at Attenuating
        Wavelengths. Journal of the Atmospheric Sciences, 11(1), p.58-67.
        DOI: 10.1175/1520-0469(1954)011<0058:EIITRM>2.0.CO;2

    .. comment
        _[1] Krämer2008 - <NAME> 2008: Quantitative Radardatenaufbereitung
        für die Niederschlagsvorhersage und die Siedlungsentwässerung,
        Mitteilungen Institut für Wasserwirtschaft, Hydrologie und
        Landwirtschaftlichen Wasserbau
        <NAME> Universität Hannover, Heft 92, ISSN 0343-8090.

    """
    if coefficients is None:
        _coefficients = {'a':1.67e-4, 'b':0.7, 'l':1.0}
    else:
        _coefficients = coefficients

    a = _coefficients['a']
    b = _coefficients['b']
    l = _coefficients['l']


    k = np.empty(gateset.shape)
    k[...,0] = 0.
    ksum = 0.

    # multidimensional version
    # assumes that iteration is only along the last dimension (i.e. range gates)
    # all other dimensions are calculated simultaneously to gain some speed
    for gate in range(gateset.shape[-1]-1):
        # calculate k in dB/km from k-Z relation
        # c.f. Krämer2008(p. 147)
        kn = a * (10.0**((gateset[...,gate] + ksum)/10.0))**b  * 2.0 * l
        #kn = 10**(log10(a)+0.1*bin*b)
        #dBkn = 10*math.log10(a) + (bin+ksum)*b + 10*math.log10(2*l)
        ksum += kn

        k[...,gate+1] = ksum
        # stop-criterion, if corrected reflectivity is larger than 59 dBZ
        overflow = (gateset[...,gate] + ksum) > thrs
        if np.any(overflow):
            if mode == 'warn':
                logger.warning('dB-sum over threshold (%3.1f)'%thrs)
            if mode == 'nan':
                k[gate,overflow] = np.nan
            if mode == 'zero':
                k[gate,overflow] = 0.0
            else:
                raise AttenuationOverflowError

    return k


def correctAttenuationKraemer(gateset, coefficients = None, mode = 'zero',
                              thrs_dBZ = 59.0):
    """Gate-by-Gate attenuation correction according to <NAME>
    [Kraemer2008]_



    Parameters
    ----------
    gateset : array
        multidimensional array. The range gates (over which iteration has to
        be performed) are supposed to vary along
        the *last* dimension so, e.g., for a set of `l` radar images stored in
        polar form with `m` azimuths and `n` range-bins the input array's
        shape can be either (l,m,n) or (m,l,n)
        data havs to be provided in decibel representation of reflectivity (dBZ)
    coefficients : dictionary
        correction coefficients

        a_max: initial value for linear coefficient of the k-Z relation
        ( :math:`k=a*Z^{b}` )

        a_min: minimal allowed linear coefficient of the k-Z relation
        ( :math:`k=a*Z^{b}` ) in the downwards iteration of a in case of signal
        overflow (sum of signal and attenuation exceeds the threshold ``thrs``)

        b: exponent of the k-Z relation ( :math:`k=a*Z^{b}` )

        n: number of iterations from a_max to a_min

        l: length of a range gate.

        if set to None the following default dictionary will be used
        {'a_max':1.67e-4, 'a_min':2.33e-5, 'b':0.70, 'n':30, 'l':1.0}

    mode : string
        controls how the function reacts in case of signal overflow (sum of
        signal and attenuation exceeds the threshold ``thrs``)
        Possible values:

        'warn' : emit a warning through the module's logger but continue
        execution

        'zero' : set offending gates to 0.0

        'nan' : set offending gates to nan

        Any other mode will raise an Exception

    thrs_dBZ : float
        threshold, for the attenuation corrected signal in dBZ, which is deemed
        unplausible.

    Returns
    -------
    k : array
        Array with the same shape as ``gateset`` containing the calculated
        attenuation for each range gate.

    Raises
    ------
    AttenuationOverflowError
        Exception, if attenuation exceeds ``thrs`` even with smallest possible
        linear coefficient (a_min) and no handling ``mode`` is set.

    References
    ----------

    .. [Kraemer2008] <NAME> 2008: Quantitative Radardatenaufbereitung
        für die Niederschlagsvorhersage und die Siedlungsentwässerung,
        Mitteilungen Institut für Wasserwirtschaft, Hydrologie und
        Landwirtschaftlichen Wasserbau
        <NAME> Universität Hannover, Heft 92, ISSN 0343-8090.

    """
    if coefficients is None:
        _coefficients = {'a_max':1.67e-4,
                         'a_min':2.33e-5,
                         'b':0.7,
                         'n':30,
                         'l':1.0}
    else:
        _coefficients = coefficients

    a_max = _coefficients['a_max']
    a_min = _coefficients['a_min']
    b     = _coefficients['b']
    n     = _coefficients['n']
    l     = _coefficients['l']

    if np.max(np.isnan(gateset)): raise Exception('There are not processable NaN in the gateset!')

    da = (a_max - a_min) / (n - 1)
    ai = a_max + da
    k = np.zeros(gateset.shape)
    k[...,0] = 0.0
    # indexing all rows of last dimension (radarbeams)
    beams2correct = np.where(np.max(k, axis = k.ndim - 1) > (-1.))
    # iterate over possible a-parameters
    for i in range(n):
        ai = ai - da
        # subset of beams that have to be corrected and corresponding attenuations
        sub_gateset = gateset[beams2correct]
        sub_k = k[beams2correct]
        for gate in range(gateset.shape[-1] - 1):
            kn = ai * (10.0**((sub_gateset[...,gate] + sub_k[...,gate])/10.0))**b  * 2.0 * l
            sub_k[...,gate + 1] = sub_k[...,gate] + kn
        # integration of the calculated attenuation subset to the whole attenuation matrix
        k[beams2correct] = sub_k
        # indexing the rows of the last dimension (radarbeam), if any corrected values exceed the threshold
        beams2correct = np.where(np.max(gateset + k, axis = k.ndim - 1) > thrs_dBZ)
        # if there is no beam left for correction, the iteration can be interrupted prematurely
        if len(k[beams2correct]) == 0: break
    if len(k[beams2correct]) > 0:
        if mode == 'warn': logger.warning('dB-sum over threshold (%3.1f)'%thrs)
        elif mode == 'nan':  k[beams2correct] = np.nan
        else: k[beams2correct] = 0.0

    return k


def correctAttenuationHJ(gateset, coefficients = None, mode = 'zero',
                              thrs_dBZ = 59.0, max_PIA = 20.0):
    """Gate-by-Gate attenuation correction based on <NAME>
    [Kraemer2008]_, expanded by <NAME> and <NAME>
    [Jacobi2011]_



    Parameters
    ----------
    gateset : array
        multidimensional array. The range gates (over which iteration has to
        be performed) are supposed to vary along
        the *last* dimension so, e.g., for a set of `l` radar images stored in
        polar form with `m` azimuths and `n` range-bins the input array's
        shape can be either (l,m,n) or (m,l,n)
        data havs to be provided in decibel representation of reflectivity (dBZ)
    coefficients : dictionary
        correction coefficients

        a_max: initial value for linear coefficient of the k-Z relation
        ( :math:`k=a*Z^{b}` )

        a_min: minimal allowed linear coefficient of the k-Z relation
        ( :math:`k=a*Z^{b}` ) in the downwards iteration of a in case of signal
        overflow (sum of signal and attenuation exceeds the threshold ``thrs``)

        b: exponent of the k-Z relation ( :math:`k=a*Z^{b}` )

        n: number of iterations from a_max to a_min

        l: length of a range gate.

        if set to None the following default dictionary will be used
        {'a_max':1.67e-4, 'a_min':2.33e-5, 'b':0.70, 'n':30, 'l':1.0}

    mode : string
        controls how the function reacts in case of signal overflow (sum of
        signal and attenuation exceeds the threshold ``thrs``)
        Possible values:

        'warn' : emit a warning through the module's logger but continue
        execution

        'zero' : set offending gates to 0.0

        'nan' : set offending gates to nan

        Any other mode will raise an Exception

    thrs_dBZ : float
        threshold, for the attenuation corrected signal in dBZ, which is deemed
        unplausible.

    max_PIA : float
        threshold, for the attenuation corrected signal in dBZ, which is deemed
        unplausible.

    Returns
    -------
    k : array
        Array with the same shape as ``gateset`` containing the calculated
        attenuation for each range gate

    Raises
    ------
    AttenuationOverflowError
        Exception, if attenuation exceeds ``thrs`` even with smallest possible
        linear coefficient (a_min) and no handling ``mode`` is set.

    References
    ----------

    .. [Kraemer2008] Krämer, Stefan 2008: Quantitative Radardatenaufbereitung
        für die Niederschlagsvorhersage und die Siedlungsentwässerung,
        Mitteilungen Institut für Wasserwirtschaft, Hydrologie und
        Landwirtschaftlichen Wasserbau
        <NAME> Universität Hannover, Heft 92, ISSN 0343-8090

    """
    if coefficients is None:
        _coefficients = {'a_max':1.67e-4,
                         'a_min':2.33e-5,
                         'b':0.7,
                         'n':30,
                         'l':1.0}
    else:
        _coefficients = coefficients

    a_max = _coefficients['a_max']
    a_min = _coefficients['a_min']
    b     = _coefficients['b']
    n     = _coefficients['n']
    l     = _coefficients['l']

    if np.max(np.isnan(gateset)): raise Exception('There are not processable NaN in the gateset!')

    da = (a_max - a_min) / (n - 1)
    ai = a_max + da
    k = np.zeros(gateset.shape)
    k[...,0] = 0.0
    # indexing all rows of last dimension (radarbeams)
    beams2correct = np.where(np.max(k, axis = k.ndim - 1) > (-1.))
    # iterate over possible a-parameters
    for i in range(n):
        ai = ai - da
        # subset of beams that have to be corrected and corresponding attenuations
        sub_gateset = gateset[beams2correct]
        sub_k = k[beams2correct]
        for gate in range(gateset.shape[-1] - 1):
            kn = ai * (10.0**((sub_gateset[...,gate] + sub_k[...,gate])/10.0))**b  * 2.0 * l
            sub_k[...,gate + 1] = sub_k[...,gate] + kn
        # integration of the calculated attenuation subset to the whole attenuation matrix
        k[beams2correct] = sub_k
        # indexing the rows of the last dimension (radarbeam), if any corrected values exceed the threshold
        beams2correct = np.where(np.max(gateset + k, axis = k.ndim - 1) > thrs_dBZ)
        # if there is no beam left for correction, the iteration can be interrupted prematurely
        if len(k[beams2correct]) == 0: break
    if len(k[beams2correct]) > 0:
        if mode == 'warn': logger.warning('dB-sum over threshold (%3.1f)'%thrs)
        elif mode == 'nan':  k[beams2correct] = np.nan
        else: k[beams2correct] = 0.0

    return k

if __name__ == '__main__':
    print 'wradlib: Calling module <atten> as main...'
{'fraction_non_alphanumeric': 0.12452664507985292, 'fraction_numerical': 0.025849294769771144, 'mean_word_length': 2.13200412512891, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7414521', 'n_tokens_mistral': 5919, 'n_tokens_neox': 5403, 'n_words': 1832}