Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH: Rework ZR #408

Merged
merged 3 commits into from Feb 20, 2020
Merged

ENH: Rework ZR #408

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d07e365
ENH: Rework ZR
kmuehlbauer Feb 18, 2020
126c769
MIN: fix minor style and coding issues
kmuehlbauer Feb 19, 2020
ec0ac08
ENH: make shower index output triggered by keyword
kmuehlbauer Feb 19, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
19 changes: 9 additions & 10 deletions wradlib/tests/test_zr.py
View file
Open in desktop
Expand Up @@ -18,6 +18,7 @@ def setUp(self):
img[2:3, 6:9] = 2.0
img[3:4, 1:7] = 5.0

img = np.stack([img, img * 1.1, img * 1.2], axis=0)
self.img = img

def test_z_to_r(self):
Expand All @@ -27,12 +28,10 @@ def test_r_to_z(self):
self.assertEqual(zr.r_to_z(0.15), 9.611164492610417)

def test_z_to_r_enhanced(self):
res_rr, res_si = zr.z_to_r_enhanced(trafo.idecibel(self.img))
res_rr2, res_si2 = zr.z_to_r_enhanced(trafo.idecibel(self.img),
algo='mdfilt', mode='mirror')
res_rr3, res_si3 = zr.z_to_r_enhanced(trafo.idecibel(self.img),
algo='mdcorr', xmode='mirror',
ymode='mirror')
res_rr, res_si = zr.z_to_r_enhanced(trafo.idecibel(self.img),
polar=True)
res_rr2, res_si2 = zr.z_to_r_enhanced(trafo.idecibel(self.img[0]),
polar=True)

rr = np.array([[3.64633237e-02, 1.77564547e-01, 1.77564547e-01,
3.17838962e-02, 3.17838962e-02, 1.62407903e-02,
Expand Down Expand Up @@ -69,11 +68,11 @@ def test_z_to_r_enhanced(self):
[4.57142861, 4.00000004, 4.00000004, 3.08333336,
1.25000001, 8.75000003, 12.50000004, 12.08333337,
11.25000003, 7.50000002, 0.]])
np.testing.assert_allclose(rr, res_rr)
np.testing.assert_allclose(rr, res_rr2)
np.testing.assert_allclose(rr, res_rr3)

self.assertTrue(np.allclose(si, res_si))
np.testing.assert_almost_equal(si, res_si[0], decimal=6)
np.testing.assert_array_almost_equal(rr, res_rr[0], decimal=6)
np.testing.assert_almost_equal(si, res_si2, decimal=6)
np.testing.assert_array_almost_equal(rr, res_rr2, decimal=6)


if __name__ == '__main__':
Expand Down
309 changes: 84 additions & 225 deletions wradlib/zr.py
View file
Open in desktop
Expand Up @@ -21,7 +21,7 @@
__doc__ = __doc__.format('\n '.join(__all__))

import numpy as np
import scipy.ndimage.filters as filters
from scipy import signal

from wradlib import trafo

Expand Down Expand Up @@ -88,13 +88,16 @@ def r_to_z(r, a=200., b=1.6):
return a * r ** b


def _z_to_r_enhanced(z):
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).
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).

This function does the actual calculations without any padding.
Neighborhood-means are taken only for available data, reducing the number
of elements used near the edges of the array.
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
Expand All @@ -116,210 +119,16 @@ def _z_to_r_enhanced(z):
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.
"""
# get the shape of the input
dimy = z.shape[0]
dimx = z.shape[1]

# calculate the decibel values from the input
db = trafo.decibel(z)

# set up our output arrays
r = np.zeros(z.shape)
si = np.zeros(z.shape)

# calculate difference fields in x and y direction
# mainly for performance reasons, so that we can use numpy's efficient
# array operations and avoid calculating a difference more than once
diffx = np.abs(db[:, :-1] - db[:, 1:])
diffy = np.abs(db[:-1, :] - db[1:, :])

# if the reflectivity is larger than 44dBZ, then there is no need to
# calculate the shower index
gt44 = np.where(db > 44.)
r[gt44] = z_to_r(z[gt44], a=77., b=1.9)
# the same is true for values between 36.5 and 44 dBZ
bt3644 = np.where(np.logical_and(db >= 36.5, db <= 44.))
r[bt3644] = z_to_r(z[bt3644], a=200., b=1.6)

# now iterate over the array and look for the remaining values
# TODO : this could be a starting point for further optimization, if we
# iterated only over the remaining pixels instead of all
for i in range(dimy):
for j in range(dimx):
# if the reflectivity is too high, we coped with it already
# so we can skip that one
if db[i, j] >= 36.5:
# just set the shower index to some impossible value so that
# we know that there was no calculation done here
si[i, j] = -1
# continue with the next iteration
continue
else:
# calculate the bounds of the region where we have to consider
# the respective difference
xmin = max(0, j - 1)
xmax = min(dimx, j + 1)
ymin = max(0, i - 1)
ymax = min(dimy, i + 1)
# in fact python is quite forgiving with upper indices
# ours might go one index too far, so don't try to port this
# to another programming language straigt away!
diffxcut = diffx[ymin:ymax + 1, xmin:xmax]
diffycut = diffy[ymin:ymax, xmin:xmax + 1]
# calculate the mean for the current pixel
mn = (np.sum(diffxcut) + np.sum(diffycut)) / \
(diffxcut.size + diffycut.size)
# apply the three different Z/R relations
if mn < 3.5:
r[i, j] = z_to_r(z[i, j], a=125., b=1.4)
elif mn <= 7.5:
r[i, j] = z_to_r(z[i, j], a=200., b=1.6)
else:
r[i, j] = z_to_r(z[i, j], a=320., b=1.4)
# save the shower index
si[i, j] = mn
# return the results
return r, si


def z_to_r_esifilter(data):
"""calculates the shower index for the enhanced z-r relation

to be used as the callable for
:func:`scipy:scipy.ndimage.filters.generic_filter`
"""
if data[4] < 36.5:
tdata = data.reshape((3, 3))
# calculate difference fields in x and y direction
# mainly for performance reasons, so that we can use numpy's efficient
# array operations and avoid calculating a difference more than once
diffx = np.abs(tdata[:, :-1] - tdata[:, 1:])
diffy = np.abs(tdata[:-1, :] - tdata[1:, :])
return np.concatenate([diffx.ravel(), diffy.ravel()]).mean()
else:
return -1.


def _z_to_r_enhanced_mdfilt(z, mode='mirror'):
"""multidimensional version

assuming the two last dimensions represent a 2-D image
Uses :func:`scipy:scipy.ndimage.filters.generic_filter` to reduce the
number of for-loops even more.
"""
# get the shape of the input
# dimy = z.shape[-2]
# dimx = z.shape[-1]

# calculate the decibel values from the input
db = trafo.decibel(z)

# set up our output arrays
r = np.zeros(z.shape)
size = list(z.shape)
size[-2:] = [3, 3]
size[:-2] = [1] * len(size[:-2])
size = tuple(size)
si = filters.generic_filter(db, z_to_r_esifilter, size=size, mode=mode)

gt44 = db > 44.
r[gt44] = z_to_r(z[gt44], a=77, b=1.9)
si[gt44] = -1.
# the same is true for values between 36.5 and 44 dBZ
bt3644 = (db >= 36.5) & (db <= 44.)
r[bt3644] = z_to_r(z[bt3644], a=200, b=1.6)
si[bt3644] = -2.

si1 = (si >= 0.)
si2 = si1 & (si < 3.5)
si3 = si1 & ~si2 & (si <= 7.5)
si4 = si > 7.5

r[si2] = z_to_r(z[si2], a=125, b=1.4)
r[si3] = z_to_r(z[si3], a=200, b=1.6)
r[si4] = z_to_r(z[si4], a=320, b=1.4)

return r, si


def _z_to_r_enhanced_mdcorr(z, xmode='reflect', ymode='wrap'):
"""multidimensional version

assuming the two last dimensions represent a 2-D image
Uses :func:`scipy:scipy.ndimage.filters.correlate` to reduce the number of
for-loops even more.
"""
# get the shape of the input
# dimy = z.shape[-2]
# dimx = z.shape[-1]

# calculate the decibel values from the input
db = trafo.decibel(z)
# calculate the shower differences by 1-d correlation with a differencing
# kernel
db_diffx = np.abs(filters.correlate1d(db, [1, -1], axis=-1,
mode=xmode, origin=-1))
db_diffy = np.abs(filters.correlate1d(db, [1, -1], axis=-2,
mode=ymode, origin=-1))
diffxmode = 'wrap' if xmode == 'wrap' else 'constant'
diffymode = 'wrap' if ymode == 'wrap' else 'constant'
diffx_sum1 = filters.correlate1d(db_diffx, [1, 1, 1],
axis=-2, mode=diffymode)
diffxsum = filters.correlate1d(diffx_sum1, [1, 1, 0],
axis=-1, mode=diffxmode)
diffy_sum1 = filters.correlate1d(db_diffy, [1, 1, 1],
axis=-1, mode=diffxmode)
diffysum = filters.correlate1d(diffy_sum1, [1, 1, 0],
axis=-2, mode=diffymode)

divider = np.ones(db.shape) * 12.
if xmode != 'wrap':
divider[..., [0, -1]] = np.rint((divider[..., [0, -1]] + 1) /
1.618) - 1
if ymode != 'wrap':
divider[..., [0, -1], :] = np.rint((divider[..., [0, -1], :] + 1) /
1.618) - 1

# the shower index is the sum of the x- and y-differences
si = (diffxsum + diffysum) / divider

# set up our rainfall output array
r = np.zeros(z.shape)

gt44 = db > 44.
r[gt44] = z_to_r(z[gt44], a=77, b=1.9)
si[gt44] = -1.
# the same is true for values between 36.5 and 44 dBZ
bt3644 = (db >= 36.5) & (db <= 44.)
r[bt3644] = z_to_r(z[bt3644], a=200, b=1.6)
si[bt3644] = -2.

si1 = (si >= 0.)
si2 = si1 & (si < 3.5)
si3 = si1 & ~si2 & (si <= 7.5)
si4 = si > 7.5

r[si2] = z_to_r(z[si2], a=125, b=1.4)
r[si3] = z_to_r(z[si3], a=200, b=1.6)
r[si4] = z_to_r(z[si4], a=320, b=1.4)

return r, si


def z_to_r_enhanced(z, algo='plain', **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).

Parameters
----------
z : a float or an array of floats
z : ndarray a floats
Corresponds to reflectivity Z in mm**6/m**3
**must** be a 2-D array
ND-array, at least 2D
polar : bool
defaults to to True (polar data), False for cartesian data.
shower : bool
output shower index, defaults to True

Returns
-------
Expand All @@ -330,26 +139,76 @@ def z_to_r_enhanced(z, algo='plain', **kwargs):
for control purposes. May be omitted in later versions

"""
# create a padded version of the input array
padz = np.zeros((z.shape[0] + 2, z.shape[1]))
# fill the center with the original data
padz[1:-1, :] = z
# add the last beam before the first one
padz[0, :] = z[-1, :]
# add the first beam after the last one
padz[-1, :] = z[0, :]

# do the actual calculation
if algo == 'mdfilt':
padr, padsi = _z_to_r_enhanced_mdfilt(padz, **kwargs)
elif algo == 'mdcorr':
padr, padsi = _z_to_r_enhanced_mdcorr(padz, **kwargs)
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:
# fallback
padr, padsi = _z_to_r_enhanced(padz)

# return the unpadded field
return padr[1:-1, :], padsi[1:-1, :]
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.
si[:, :, -1] /= 7.
si[:, :, 1:-1] /= 12.
else:
si[:, 1:-1, 1:-1] /= 12.
si[:, :, 0] /= 7
si[:, :, -1] /= 7.
si[:, 0, :] /= 7
si[:, -1, :] /= 7.

rr = np.zeros(z.shape)
z0_ = z0[:, y_ymin:y_ymax, :]

# get masks
gt44 = (z0_ > 44)
bt3644 = (z0_ >= 36.5) & (z0_ <= 44.)
si[bt3644] = -1.
si[gt44] = -1.
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., b=1.6)
rr[mn75h] = z_to_r(z[mn75h], a=320., b=1.4)
rr[mn35] = z_to_r(z[mn35], a=125., b=1.4)
rr[gt44] = z_to_r(z[gt44], a=77., b=1.9)
rr[bt3644] = z_to_r(z[bt3644], a=200., b=1.6)

rr.shape = shape
si.shape = shape

if shower:
return rr, si
else:
return rr


if __name__ == '__main__':
Expand Down