A Typical Workflow For Radar-Based Rainfall Estimation#

Raw, unprocessed reflectivity products can already provide useful visual information about the spatial distribution of rainfall fields. However, in order to use weather radar observations for quantitative studies (e.g. in hydrological modelling or for assimilation into Numerical Weather Prediction models), the data has to be carefully processed in order to account for typical errors sources such as ground echoes (clutter), attenuation of the radar signal, or uncertainties in the Z/R relationship.

Moreover, it might be necessary to transfer the data from polar coordinates to cartesian grids, or to combine observations from different radar locations in overlapping areas on a common grid (composition). And in the end, you would typically like to visualise the spatial rainfall distribution on a map. Many users also need to quantify the potential error (uncertainty) of their data-based rainfall estimation.

These are just some steps that might be necessary in order to make radar data useful in a specific quantitative application environment. All steps together are typically referred to as a “radar data processing chain”. \(\omega radlib\) was designed to support you in establishing your own processing chain, suited to your specific requirements. In the following, we will provide an outline of a typical processing chain, step-by-step. You might not need all steps for your own workflow, or you might need steps which are not yet included here.


Consider this just as an example. We will not go into detail for each step in this section, but refer to more detailed tutorials (if available) or the corresponding entry in the library reference. Most of the steps have a corresponding \(\omega radlib\) module. In order to access the functions of \(\omega radlib\), you have to import \(\omega radlib\) in your Python environment:

import wradlib as wrl

If you have trouble with that import, please head back to the Installation section.


The data used in this tutorial can be found in the wradlib-data repository. Follow these instructions to install and use this data files.


Be aware that applying an algorithm for error correction does not guarantee that the error is totally removed. Error correction procedures are suceptible to errors, too. Not only might they fail to remove the error. They might also introduce new errors. The trade-off between costs (introduction of new errors) and benefits (error reduction) can turn out differently for different locations, different points in time, or different rainfall situations.

Reading the data#

The binary encoding of many radar products is a major obstacle for many potential radar users. Often, decoder software is not easily available. \(\omega radlib\) supports a couple of formats such as the ODIM_H5 implementation, NetCDF, and some formats used by the German Weather Service. With the porting of xarray-based code from \(\omega radlib\) to xradar the range of supported formats will surely increase in the future.

Ever since the basic data type used in \(\omega radlib\) is a multi-dimensional array, the numpy.ndarray. Such an array might e.g. represent a polar or cartesian grid, or a series of rain gage observations. Metadata are normally managed as Python dictionaries. In order to read the content of a data file into a numpy array, you would normally use the wradlib.io module. With the evolution in the scientific Python stack (eg. xarray) \(\omega radlib\) adopted DataTree, Dataset and DataArray to conveniently hold radar data together with coordinates and metadata (now available in xradar). In the following example, a local PPI from the German Weather Service, a DX file, is read and converted to an xarray.DataArray and plotted:

The metadata object can be inspected via keywords. The data object contains the actual data, in this case a polar grid with 360 azimuth angles and 128 range bins.

filename = wrl.util.get_wradlib_data_file("dx/raa00-dx_10908-0806021655-fbg---bin.gz")
data, metadata = wrl.io.read_dx(filename)
(360, 128)
dict_keys(['producttype', 'datetime', 'radarid', 'bytes', 'version', 'cluttermap', 'dopplerfilter', 'statfilter', 'elevprofile', 'message', 'elev', 'azim', 'clutter'])

Using xarray we can conveniently hold the reflectivity data as well as the needed coordinates in one convenient structure.

radar_location = (8.003611, 47.873611, 1516.0)
da = wrl.georef.create_xarray_dataarray(
    r=np.arange(500, data.shape[1] * 1000 + 500, 1000),

We can use xarray directly to create a nice plot.

<matplotlib.collections.QuadMesh at 0x7f34313a9e10>

Georeferencing and Projection#

In order to define the horizontal and vertical position of the radar bins, we need to retrieve the corresponding 3-dimensional coordinates in terms of longitude, latitude and altitude. This information is required e.g. if the positions should be plotted on a map. It is also required for constructing CAPPIs. The position of a radar bin in 3-dimensional space depends on the position of the radar device, the elevation and azimuth angle of the radar beam, the range of the bin, and the assumed influence of atmospheric refraction on the beam propagation. For the sample data used above, the position of the radar device is the Feldberg in Germany (8.005, 47.8744, 1517):

We have the ability to easily georeference all radar bins (here with an azimuthal equidistant projection). This adds the projection as crs_wkt coordinate to the Dataset/DataArray.

da = da.wrl.georef.georeference()
<xarray.DataArray (azimuth: 360, range: 128)>
array([[-10. , -18.5, -26.5, ..., -32.5, -32.5, -32.5],
       [-19. , -19.5, -14. , ..., -32.5, -32.5, -32.5],
       [-15. , -21. , -32.5, ..., -32.5, -32.5, -32.5],
       [ -0.5, -15.5, -32.5, ..., -32.5, -32.5, -32.5],
       [  1.5, -18.5, -19. , ..., -32.5, -32.5, -32.5],
       [ -2.5, -20.5, -32.5, ..., -32.5, -32.5, -32.5]])
Coordinates: (12/14)
  * range       (range) float64 500.0 1.5e+03 2.5e+03 ... 1.265e+05 1.275e+05
  * azimuth     (azimuth) float64 0.0 1.0 2.0 3.0 ... 356.0 357.0 358.0 359.0
    elevation   (azimuth) float64 0.2 0.2 0.1 0.1 0.1 ... 0.3 0.3 0.3 0.3 0.2
    longitude   float64 8.004
    latitude    float64 47.87
    altitude    float64 1.516e+03
    ...          ...
    y           (azimuth, range) float64 499.9 1.5e+03 ... 1.264e+05 1.274e+05
    z           (azimuth, range) float64 1.518e+03 1.521e+03 ... 2.918e+03
    gr          (azimuth, range) float64 499.9 1.5e+03 ... 1.265e+05 1.275e+05
    rays        (azimuth, range) float64 0.0 0.0 0.0 0.0 ... 359.0 359.0 359.0
    bins        (azimuth, range) float64 500.0 1.5e+03 ... 1.265e+05 1.275e+05
    crs_wkt     int64 0

\(\omega radlib\) supports the projection between geographical coordinates (lon/lat) and other reference systems. It uses GDAL/OSR Spatial References Objects as function parameters. Basically, you have to create the OSR-object by using GDAL-capabilities or one of the provided helper functions. We recommend the creation using EPSG numbers:

# UTM Zone 32, EPSG-Number 32632
utm = wrl.georef.epsg_to_osr(32632)
da_utm = da.wrl.georef.reproject(trg_crs=utm)
da_utm.plot(x="x", y="y")
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
            SPHEROID["WGS 84",6378137,298.257223563,
<matplotlib.collections.QuadMesh at 0x7f343123c990>

Second, you can provide a string which represents the projection - based on the PROJ library. You can look up projection strings, but for some projections, \(\omega radlib\) helps you to define a projection string. In the following example, the target projection is ‘dwd-radolan’:

radolan = wrl.georef.create_osr("dwd-radolan")
da_radolan = da.wrl.georef.reproject(trg_crs=radolan)
da_radolan.plot(x="x", y="y")
<matplotlib.collections.QuadMesh at 0x7f3428118390>


Get more info in the library reference section Georeferencing.

For plotting we also can use the \(\omega radlib\) plot-function, which detects the data as PPI and plots it.

pm = da.wrl.vis.plot(add_colorbar=True)


Get more info in the section Supported radar data formats and in the library reference section Raw Data I/O.

Clutter removal#

Clutter are non-meteorological echos. They are caused by the radar beam hitting objects on the earth’s surface (e.g. mountain or hill tops, houses, wind turbines) or in the air (e.g. airplanes, birds). These objects can potentially cause high reflectivities due large scattering cross sections. Static clutter, if not efficiently removed by Doppler filters, can cause permanent echos which could introduce severe bias in quantitative applications. Thus, an efficient identification and removal of clutter is mandatory e.g. for hydrological studies. Clutter removal can be based on static maps or dynamic filters. Normally, static clutter becomes visible more clearly in rainfall accumulation maps over periods of weeks or months. We recommend such accumulations to create static clutter maps which can in turn be used to remove the static clutter from an image and fill the resulting gaps by interpolation.

In the following example, the clutter filter published by Gabella et al., 2002) is applied to the single radar sweep of the above example:

clutter = da.wrl.classify.filter_gabella(tr1=12, n_p=6, tr2=1.1)
pm = clutter.wrl.vis.plot(cmap=plt.cm.gray)
plt.title("Clutter Map")
Text(0.5, 1.0, 'Clutter Map')

The resulting Boolean array clutter indicates the position of clutter. It can be used to interpolate the values at those positons from non-clutter values, as shown in the following lines:

data_no_clutter = da.wrl.ipol.interpolate_polar(clutter)
pm = data_no_clutter.wrl.vis.plot(add_colorbar=True)  # simple diagnostic plot

It is generally recommended to remove the clutter before e.g. gridding the data. Otherwise the clutter signal might be “smeared” over multiple grid cells, resulting into a decrease in detectability.


Get more info in the library reference section Echo Classification.

Attenuation correction#

Attenuation by wet radome and by heavy rainfall can cause serious underestimation of rainfall for C-Band and X-Band devices. For such radar devices, situations with heavy rainfall require a correction of attenuation effects. The general approach with single-polarized radars is to use a recursive gate-by-gate approach. See Kraemer et al., 2008 for an introduction to this concept. Basically, the specific attenuation k of the first range gate is computed via a so-called k-Z relationship. Based on k, the reflectivity of the second range gate is corrected and then used to compute the specific attenuation for the second range gate (and so on). The concept was first introduced by Hitschfeld et al., 1954). Its main drawback is its suceptibility to instable behaviour. \(\omega radlib\) provides different implementations which address this problem.

One example is the algorithm published by Jacobi and Heistermann, 2016:

pia = data_no_clutter.wrl.atten.correct_attenuation_constrained(
    constraints=[wrl.atten.constraint_dbz, wrl.atten.constraint_pia],
    constraint_args=[[59.0], [20.0]],
data_attcorr = data_no_clutter + pia

The first line computes the path integrated attenuation pia for each radar bin. The second line uses pia to correct the reflectivity values. Let’s inspect the effect of attenuation correction for an azimuth angle of 65 deg:

plt.figure(figsize=(10, 8))
data_no_clutter[65].plot(label="no attcorr")
<matplotlib.legend.Legend at 0x7f34311adf10>


Get more info in the library reference section Attenuation Correction. There you will learn to know the algorithms available for attenuation correction and how to manipulate their behaviour by using additonal keyword arguments.

Vertical Profile of Reflectivity#

Precipitation is 3-dimensional in space. The vertical distribution of precipitation (and thus reflectivity) is typically non-uniform. As the height of the radar beam increases with the distance from the radar location (beam elevation, earth curvature), one sweep samples from different heights. The effects of the non-uniform VPR and the different sampling heights need to be accounted for if we are interested in the precipiation near the ground or in defined altitudes.


Get more info in the library reference section Vertical Profile of Reflectivity (VPR). There you will learn how to reference polar volume data, to create CAPPIs and Pseudo CAPPIs, to inspect vertical profiles of reflectivity (UNDER DEVELOPMENT), and to use these for correction (UNDER DEVELOPMENT).

Conversion of Reflectivity into Rainfall#

Reflectivity (Z) and precipitation rate (R) can be related in form of a power law \(R=a*Z^b\). The parameters a and b depend on the type of precipitation in terms of drop size distribution and water temperature. Before applying the Z-R relationship, we need to convert from dBZ to Z:

z = data_attcorr.wrl.trafo.idecibel()
R = z.wrl.zr.z_to_r()  # wrl.zr.z_to_r(wrl.trafo.idecibel(data_attcorr))

The above line uses the default parameters a=200 and b=1.6 for the Z-R relationship. In order to compute a rainfall depth from rainfall intensity, we have to specify an integration interval in seconds. In this example, we choose five minutes (300 s), corresponding to the scan interval.

depths = R.wrl.trafo.r_to_depth(300)
<matplotlib.collections.QuadMesh at 0x7f3421a9e6d0>


Get more info in the section Converting reflectivity to rainfall and in the library reference sections Z-R Conversions and Data Transformation. Here you will learn about the effects of the Z-R parameters a and b.

Rainfall accumulation#

For many applications, accumulated rainfall depths over specific time intervals are required, e.g. hourly or daily accumulations. In the following example, we will use a synthetic time series of 5 minute intervals. Just imagine we have repeated the above procedure for one day of five-minute sweeps and combined the arrays of rainfall depth in a 3-dimensional array of shape (number of time steps, number of azimuth angles, number of range gates).

Now we want to compute the daily rainfall sum:

sweep_times = wrl.util.from_to("2012-10-26 00:00:00", "2012-10-27 00:00:00", 300)
depths_5min = np.random.uniform(size=(len(sweep_times) - 1, 360, 128))
depth_daily = depths_5min.sum(axis=0)

Check the shape the resulting array for plausibility:

(360, 128)


For more advanced operations on time series, we recommend the Pandas package.


Assume you would like to transfer the rainfall intensity from the above example from polar coordinates to a Cartesian grid, or to an arbitrary set of irregular points in space (e.g. centroids of sub-catchments). You already retrieved the Cartesian coordinates of the radar bins in the previous section Georeferencing and Projection. Now you only need to define the target coordinates (e.g. a grid) and apply the togrid function of the wradlib.comp module. In this example, we want our grid only to represent the South-West sector of our radar circle on a 100 x 100 grid. First, we define the target grid coordinates (these must be an array of 100x100 rows with one coordinate pair each):

xgrid = np.linspace(R.x.min(), R.x.mean(), 100)
ygrid = np.linspace(R.y.min(), R.y.mean(), 100)
cart = xr.Dataset(coords={"x": (["x"], xgrid), "y": (["y"], ygrid)})
Dimensions:  (x: 100, y: 100)
  * x        (x) float64 -1.275e+05 -1.262e+05 -1.249e+05 ... -1.287e+03 0.5076
  * y        (y) float64 -1.275e+05 -1.262e+05 -1.249e+05 ... -1.287e+03 0.7419
Data variables:

Now we transfer the polar data to the grid and mask out invalid values for plotting (values outside the radar circle receive NaN):

gridded = data_attcorr.wrl.comp.togrid(
    cart, radius=128000.0, center=(0, 0), interpol=wrl.ipol.Nearest
<xarray.DataArray 'None.togrid' (y: 100, x: 100)>
array([[         nan,          nan,          nan, ..., -32.48363104,
        -32.48363104, -32.48289347],
       [         nan,          nan,          nan, ..., -32.48363282,
        -32.48363282, -32.48289525],
       [         nan,          nan,          nan, ..., -32.48363637,
        -32.48363637, -32.48289881],
       [ 23.14862042,  23.63355105,  20.61339187, ..., -32.49999468,
        -32.49999645, -32.49999645],
       [ 23.14862042,  23.63355105,  20.61339187, ..., -32.49999645,
        -32.49999823, -32.49999823],
       [ 22.09449356,  21.58366419,  18.06968247, ..., -32.4999826 ,
        -18.99999823, -32.5       ]])
  * y        (y) float64 -1.275e+05 -1.262e+05 -1.249e+05 ... -1.287e+03 0.7419
  * x        (x) float64 -1.275e+05 -1.262e+05 -1.249e+05 ... -1.287e+03 0.5076
fig = plt.figure(figsize=(10, 8))
ax = plt.subplot(111, aspect="equal")
pm = gridded.plot()
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
Text(0, 0.5, 'Northing (m)')


Get more info about the function wradlib.comp.togrid().

Adjustment by rain gage observations#

Adjustment normally refers to using rain gage observations on the ground to correct for errors in the radar-based rainfall estimating. Goudenhoofdt et al., 2009 provide an excellent overview of adjustment procedures. A typical approach is to quantify the error of the radar-based rainfall estimate at the rain gage locations, assuming the rain gage observation to be accurate. The error can be assumed to be additive, multiplicative, or a mixture of both. Most approaches assume the error to be heterogeneous in space. Hence, the error at the rain gage locations will be interpolated to the radar bin (or grid) locations and then used to adjust (correct) the raw radar rainfall estimates.

In the following example, we will use an illustrative one-dimensional example with synthetic data (just imagine radar rainfall estimates and rain gage observations along one radar beam).

First, we create the synthetic “true” rainfall (truth):

import numpy as np

radar_coords = np.arange(0, 101)
truth = np.abs(1.5 + np.sin(0.075 * radar_coords)) + np.random.uniform(
    -0.1, 0.1, len(radar_coords)

The radar rainfall estimate radar is then computed by imprinting a multiplicative error on truth and adding some noise:

error = 0.75 + 0.015 * radar_coords
radar = error * truth + np.random.uniform(-0.1, 0.1, len(radar_coords))

Synthetic gage observations obs are then created by selecting arbitrary “true” values:

obs_coords = np.array([5, 10, 15, 20, 30, 45, 65, 70, 77, 90])
obs = truth[obs_coords]

Now we adjust the radar rainfall estimate by using the gage observations. First, you create an “adjustment object” from the approach you want to use for adjustment. After that, you can call the object with the actual data that is to be adjusted. Here, we use a multiplicative error model with spatially heterogenous error (see wradlib.adjust.AdjustMultiply():

adjuster = wrl.adjust.AdjustMultiply(obs_coords, radar_coords, nnear_raws=3)
adjusted = adjuster(obs, radar)

Let’s compare the truth, the radar rainfall estimate and the adjusted product:

plt.plot(radar_coords, truth, "k-", label="True rainfall", linewidth=2.0)
plt.xlabel("Distance (km)")
plt.ylabel("Rainfall intensity (mm/h)")
    label="Raw radar rainfall",
    label="Gage observation",
    label="Multiplicative adjustment",
plt.legend(prop={"size": 12})
<matplotlib.legend.Legend at 0x7f3421e3ef50>


Get more info in the library reference section Gage Adjustment. There, you will also learn how to use the built-in cross-validation in order to evaluate the performance of the adjustment approach.

Verification and quality control#

Typically, radar-based precipitation estimation and the effectiveness of the underlying correction and adjustment methods are verified by comparing the results against rain gage observations on the ground. wradlib.verify module provides procedures not only to extract the radar values at specific gauge locations, but also a set of error metrics which are computed from gage observations and the corresponding radar-based precipitation estimates (including standard metrics such as RMSE, mean error, Nash-Sutcliffe Efficiency). In the following, we will illustrate the usage of error metrics by comparing the “true” rainfall against the raw and adjusted radar rainfall estimates from the above example:

raw_error = wrl.verify.ErrorMetrics(truth, radar)
adj_error = wrl.verify.ErrorMetrics(truth, adjusted)

Error metrics can be reported e.g. as follows:

print("Error metrics for unadjusted radar rainfall estimates:")
print("\nError metrics for adjusted radar rainfall estimates:")
Error metrics for unadjusted radar rainfall estimates:
{'corr': 0.69,
 'mas': 0.77,
 'meanerr': 0.69,
 'mse': 1.19,
 'nash': -1.54,
 'pbias': 43.4,
 'r2': 0.48,
 'ratio': 1.5,
 'rmse': 1.09,
 'spearman': 0.77,
 'sse': 120.29}

Error metrics for adjusted radar rainfall estimates:
{'corr': 1.0,
 'mas': 0.06,
 'meanerr': 0.03,
 'mse': 0.01,
 'nash': 0.98,
 'pbias': 1.9,
 'r2': 0.99,
 'ratio': 1.01,
 'rmse': 0.1,
 'spearman': 0.99,
 'sse': 0.79}


Get more info in the library reference section Verification.

Visualisation and mapping#

In the above sections Reading the data, Clutter removal, and Gridding you already saw examples of the \(\omega radlib's\) plotting capabilities.


Get more info in the library reference section Visualisation.

Data export to other applications#

Once you created a dataset which meets your requirements, you might want to export it to other applications or archives. \(\omega radlib\) does not favour or support a specific output format. Basically, you have all the freedom of choice offered by Python and its packages in order to export your data. Arrays can be stored as text or binary files by using numpy functions. You can use the package NetCDF4 to write NetCDF files, and the packages h5py or PyTables to write hdf5 files.

Using xradar a standardized data export is available. Of course, you can also export data as images. See Visualisation for some options.

Export your data array as a text file:

np.savetxt("mydata.txt", data_attcorr)

Or as a gzip-compressed text file:

np.savetxt("mydata.gz", data_attcorr)

Or as a NetCDF file:

data_attcorr.name = "reflectivity corrected"
data_attcorr.to_netcdf("test.nc", group="sweep_0")
with xr.open_dataset("test.nc", group="sweep_0") as ds:
Dimensions:                 (range: 128, azimuth: 360)
Coordinates: (12/14)
  * range                   (range) float64 500.0 1.5e+03 ... 1.275e+05
  * azimuth                 (azimuth) float64 0.0 1.0 2.0 ... 357.0 358.0 359.0
    elevation               (azimuth) float64 ...
    longitude               float64 ...
    latitude                float64 ...
    altitude                float64 ...
    ...                      ...
    y                       (azimuth, range) float64 ...
    z                       (azimuth, range) float64 ...
    gr                      (azimuth, range) float64 ...
    rays                    (azimuth, range) float64 ...
    bins                    (azimuth, range) float64 ...
    crs_wkt                 int64 ...
Data variables:
    reflectivity corrected  (azimuth, range) float64 ...


An example for hdf5 export will follow.