Source code for wradlib.zr

#!/usr/bin/env python
# Copyright (c) 2011-2023, wradlib developers.
# Distributed under the MIT License. See LICENSE.txt for more info.

"""
Z-R Conversions
^^^^^^^^^^^^^^^

Module zr takes care of transforming reflectivity
into rainfall rates and vice versa

.. currentmodule:: wradlib.zr

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

   {}
"""
__all__ = ["z_to_r", "r_to_z", "z_to_r_enhanced", "ZRMethods"]
__doc__ = __doc__.format("\n   ".join(__all__))

from functools import singledispatch

import numpy as np
import xarray as xr
from scipy import signal
from xradar.model import sweep_vars_mapping

from wradlib import trafo, util


[docs]@singledispatch def z_to_r(z, *, a=200.0, b=1.6): """Conversion from reflectivities to rain rates. Calculates rain rates from radar reflectivities using a power law Z/R relationship Z = a*R**b Parameters ---------- z : float or :class:`numpy:numpy.ndarray` Corresponds to reflectivity Z in mm**6/m**3 a : float Parameter ``a`` of the Z/R-relationship Standard value according to Marshall-Palmer is a=200., b=1.6 b : float Parameter ``b`` of the Z/R-relationship Standard value according to Marshall-Palmer is b=1.6 Note ---- The German Weather Service uses a=256 and b=1.42 instead of the Marshall-Palmer defaults. Returns ------- output : float or :class:`numpy:numpy.ndarray` rainfall intensity in mm/h """ return (z / a) ** (1.0 / b)
@z_to_r.register(xr.DataArray) def _z_to_r_xarray(obj, **kwargs): """Conversion from reflectivities to rain rates. Calculates rain rates from radar reflectivities using a power law Z/R relationship Z = a*R**b Parameters ---------- obj : :py:class:`xarray:xarray.DataArray` Corresponds to reflectivity Z in mm**6/m**3 Keyword Arguments ----------------- a : float Parameter `a` of the Z/R-relationship Standard value according to Marshall-Palmer is a=200., b=1.6 b : float Parameter `b` of the Z/R-relationship Standard value according to Marshall-Palmer is b=1.6 Note ---- The German Weather Service uses a=256 and b=1.42 instead of the Marshall-Palmer defaults. Returns ------- output : :py:class:`xarray:xarray.DataArray` rainfall intensity in mm/h """ dim0 = obj.wrl.util.dim0() out = xr.apply_ufunc( z_to_r, obj, input_core_dims=[[dim0, "range"]], output_core_dims=[[dim0, "range"]], kwargs=kwargs, dask="parallelized", dask_gufunc_kwargs=dict(allow_rechunk=True), ) out.attrs = sweep_vars_mapping["RATE"] out.name = "RATE" return out
[docs]@singledispatch def r_to_z(r, *, a=200.0, b=1.6): """Calculates reflectivity from rain rates using a power law Z/R relationship Z = a*R**b Parameters ---------- r : float or :class:`numpy:numpy.ndarray` Corresponds to rainfall intensity in mm/h a : float Parameter ``a`` of the Z/R-relationship Standard value according to Marshall-Palmer is a=200., b=1.6 b : float Parameter ``b`` of the Z/R-relationship Standard value according to Marshall-Palmer is b=1.6 Note ---- The German Weather Service uses a=256 and b=1.42 instead of the Marshall-Palmer defaults. Returns ------- output : float or :class:`numpy:numpy.ndarray` reflectivity in mm**6/m**3 """ return a * r**b
@r_to_z.register(xr.DataArray) def _r_to_z_xarray(obj, **kwargs): """Calculates reflectivity from rain rates using a power law Z/R relationship Z = a*R**b Parameters ---------- obj : :py:class:`xarray:xarray.DataArray` Corresponds to rainfall intensity in mm/h a : float Parameter ``a`` of the Z/R-relationship Standard value according to Marshall-Palmer is a=200., b=1.6 b : float Parameter ``b`` of the Z/R-relationship Standard value according to Marshall-Palmer is b=1.6 Note ---- The German Weather Service uses a=256 and b=1.42 instead of the Marshall-Palmer defaults. Returns ------- output : :py:class:`xarray:xarray.DataArray` reflectivity in mm**6/m**3 """ dim0 = obj.wrl.util.dim0() out = xr.apply_ufunc( r_to_z, obj, input_core_dims=[[dim0, "range"]], output_core_dims=[[dim0, "range"]], kwargs=kwargs, dask="parallelized", dask_gufunc_kwargs=dict(allow_rechunk=True), ) out.attrs = sweep_vars_mapping["ZH"] out.name = "ZH" return out
[docs]@singledispatch def z_to_r_enhanced(z, *, polar=True, shower=True): """Calculates rainrates from radar reflectivities using the enhanced \ three-part Z-R-relationship used by the DWD (as of 2009) To be used with polar representations so that one dimension is cyclical. i.e. z should be of shape (nazimuths, nbins) --> the first dimension is the cyclical one. For DWD DX-Data z's shape is (360,128). Neighborhood-means are taken only for available data via fast convolution sums. Refer to the RADOLAN final report or the RADOLAN System handbook for details on the calculations. Basically, for low reflectivities an index called the shower index is calculated as the mean of the differences along both axis in a neighborhood of 3x3 pixels. This means: +------+-----------------+ | | x-direction --> | +------+-----+-----+-----+ | | y | 1 | 2 | 3 | | | l +-----+-----+-----+ | | d | 4 | 5 | 6 | | | i +-----+-----+-----+ | | r | 7 | 8 | 9 | +------+-----+-----+-----+ If 5 is the pixel in question, it's shower index is calculated as: .. math:: ( &|1-2| + |2-3| + |4-5| + |5-6| + |7-8| + |8-9| + \\\\ &|1-4| + |4-7| + |2-5| + |5-8| + |3-6| + |6-9| ) / 12. then, the upper line of the sum would be diffx (DIFFerences in X-direction), the lower line would be diffy (DIFFerences in Y-direction) in the code below. Parameters ---------- z : :class:`numpy:numpy.ndarray` Corresponds to reflectivity Z in mm**6/m**3 ND-array, at least 2D polar : bool defaults to True for polar data, False for cartesian data. shower : bool output shower index, defaults to True Returns ------- r : :class:`numpy:numpy.ndarray` r - array of shape z.shape - calculated rain rates si : :class:`numpy:numpy.ndarray` si - array of shape z.shape - calculated shower index for control purposes. May be omitted in later versions """ z = np.asanyarray(z, dtype=np.float64) shape = z.shape z = z.reshape((-1,) + shape[-2:]) if polar: z0 = np.concatenate([z[:, -1:, :], z, z[:, 0:1, :]], axis=-2) x_ymin, x_ymax = 2, -2 y_ymin, y_ymax = 1, -1 x_xmin, x_xmax = None, None y_xmin, y_xmax = 1, -1 else: z0 = z.copy() x_ymin, x_ymax = 1, -1 y_ymin, y_ymax = None, None x_xmin, x_xmax = None, None y_xmin, y_xmax = 1, -1 z0 = trafo.decibel(z0) # create shower index using differences and convolution sum diffx = np.abs(np.diff(z0, n=1, axis=-1)) diffy = np.abs(np.diff(z0, n=1, axis=-2)) xkernel = np.ones((1, 3, 2)) ykernel = np.ones((1, 2, 3)) resx = signal.convolve(diffx, xkernel, mode="full", method="direct")[ :, x_ymin:x_ymax, x_xmin:x_xmax ] resy = signal.convolve(diffy, ykernel, mode="full", method="direct")[ :, y_ymin:y_ymax, y_xmin:y_xmax ] si = resx + resy # edge cases divide by 7, everything else divide by 12 if polar: si[:, :, 0] /= 7.0 si[:, :, -1] /= 7.0 si[:, :, 1:-1] /= 12.0 else: si[:, 1:-1, 1:-1] /= 12.0 si[:, :, 0] /= 7 si[:, :, -1] /= 7.0 si[:, 0, :] /= 7 si[:, -1, :] /= 7.0 rr = np.zeros(z.shape) z0_ = z0[:, y_ymin:y_ymax, :] # get masks gt44 = z0_ > 44 bt3644 = (z0_ >= 36.5) & (z0_ <= 44.0) si[bt3644] = -1.0 si[gt44] = -1.0 mn75h = (si > 7.5) & (si < 36.5) mn35 = (si > -1) & (si < 3.5) mn75l = (si >= 3.5) & (si <= 7.5) # calculate rainrates according DWD rr[mn75l] = z_to_r(z[mn75l], a=200.0, b=1.6) rr[mn75h] = z_to_r(z[mn75h], a=320.0, b=1.4) rr[mn35] = z_to_r(z[mn35], a=125.0, b=1.4) rr[gt44] = z_to_r(z[gt44], a=77.0, b=1.9) rr[bt3644] = z_to_r(z[bt3644], a=200.0, b=1.6) rr.shape = shape si.shape = shape if shower: return rr, si else: return rr
@z_to_r_enhanced.register(xr.DataArray) def _z_to_r_enhanced_xarray(obj, **kwargs): """Calculates rainrates from radar reflectivities using the enhanced \ three-part Z-R-relationship used by the DWD (as of 2009) To be used with polar representations so that one dimension is cyclical. i.e. z should be of shape (nazimuths, nbins) --> the first dimension is the cyclical one. For DWD DX-Data z's shape is (360,128). Neighborhood-means are taken only for available data via fast convolution sums. Refer to the RADOLAN final report or the RADOLAN System handbook for details on the calculations. Basically, for low reflectivities an index called the shower index is calculated as the mean of the differences along both axis in a neighborhood of 3x3 pixels. This means: +------+-----------------+ | | x-direction --> | +------+-----+-----+-----+ | | y | 1 | 2 | 3 | | | l +-----+-----+-----+ | | d | 4 | 5 | 6 | | | i +-----+-----+-----+ | | r | 7 | 8 | 9 | +------+-----+-----+-----+ If 5 is the pixel in question, it's shower index is calculated as: .. math:: ( &|1-2| + |2-3| + |4-5| + |5-6| + |7-8| + |8-9| + \\\\ &|1-4| + |4-7| + |2-5| + |5-8| + |3-6| + |6-9| ) / 12. then, the upper line of the sum would be diffx (DIFFerences in X-direction), the lower line would be diffy (DIFFerences in Y-direction) in the code below. Parameters ---------- obj : :py:class:`xarray:xarray.DataArray` Corresponds to reflectivity Z in mm**6/m**3 ND-array, at least 2D Keyword Arguments ----------------- shower : bool output shower index, defaults to True Returns ------- out : :py:class:`xarray:xarray.DataArray` calculated rain rates si : :py:class:`xarray:xarray.DataArray` calculated shower index """ dim0 = obj.wrl.util.dim0() shower = kwargs.setdefault("shower", True) output_core_dims = [[dim0, "range"]] if shower is True: output_core_dims *= 2 out = xr.apply_ufunc( z_to_r_enhanced, obj, input_core_dims=[[dim0, "range"]], output_core_dims=[[dim0, "range"], [dim0, "range"]], dask="parallelized", kwargs=kwargs, dask_gufunc_kwargs=dict(allow_rechunk=True), ) if shower is True: out, si = out out.attrs = sweep_vars_mapping["RATE"] out.name = "RATE" return out, si
[docs]class ZRMethods(util.XarrayMethods): """wradlib xarray SubAccessor methods for zr."""
[docs] @util.docstring(_z_to_r_xarray) def z_to_r(self, *args, **kwargs): if not isinstance(self, ZRMethods): return z_to_r(self, *args, **kwargs) else: return z_to_r(self._obj, *args, **kwargs)
[docs] @util.docstring(_r_to_z_xarray) def r_to_z(self, *args, **kwargs): if not isinstance(self, ZRMethods): return r_to_z(self, *args, **kwargs) else: return r_to_z(self._obj, *args, **kwargs)
[docs] @util.docstring(_z_to_r_enhanced_xarray) def z_to_r_enhanced(self, *args, **kwargs): if not isinstance(self, ZRMethods): return z_to_r_enhanced(self, *args, **kwargs) else: return z_to_r_enhanced(self._obj, *args, **kwargs)
if __name__ == "__main__": print("wradlib: Calling module <zr> as main...")