Fuzzy echo classification from dual-pol moments

Fuzzy echo classification from dual-pol moments#

import warnings

import matplotlib.pyplot as plt
import numpy as np
import wradlib
import xarray as xr
from IPython.display import display
from wradlib.util import get_wradlib_data_file

warnings.filterwarnings("ignore")

Setting the file paths#

rhofile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-R.nc")
phifile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-P.nc")
reffile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-Z.nc")
dopfile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-V.nc")
zdrfile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-D.nc")
mapfile = get_wradlib_data_file("hdf5/TAG_cmap_sweeps_0204050607.hdf5")

Downloading file 'netcdf/TAG-20120801-140046-02-R.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/netcdf/TAG-20120801-140046-02-R.nc' to '/home/docs/.cache/wradlib-data'.

Downloading file 'netcdf/TAG-20120801-140046-02-P.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/netcdf/TAG-20120801-140046-02-P.nc' to '/home/docs/.cache/wradlib-data'.

Downloading file 'netcdf/TAG-20120801-140046-02-Z.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/netcdf/TAG-20120801-140046-02-Z.nc' to '/home/docs/.cache/wradlib-data'.

Downloading file 'netcdf/TAG-20120801-140046-02-V.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/netcdf/TAG-20120801-140046-02-V.nc' to '/home/docs/.cache/wradlib-data'.

Downloading file 'netcdf/TAG-20120801-140046-02-D.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/netcdf/TAG-20120801-140046-02-D.nc' to '/home/docs/.cache/wradlib-data'.

Downloading file 'hdf5/TAG_cmap_sweeps_0204050607.hdf5' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/TAG_cmap_sweeps_0204050607.hdf5' to '/home/docs/.cache/wradlib-data'.

Read the data#

This reads teh moments and a precomputed static clutter map. The data is organized as a dictionary.

dat = {}
dat["rho"], attrs_rho = wradlib.io.read_edge_netcdf(rhofile)
dat["phi"], attrs_phi = wradlib.io.read_edge_netcdf(phifile)
dat["ref"], attrs_ref = wradlib.io.read_edge_netcdf(reffile)
dat["dop"], attrs_dop = wradlib.io.read_edge_netcdf(dopfile)
dat["zdr"], attrs_zdr = wradlib.io.read_edge_netcdf(zdrfile)
dat["map"] = wradlib.io.from_hdf5(mapfile)[0][0]

dat = {k: (["azimuth", "range"], v) for k, v in dat.items()}

Create an xarray.Dataset holding the data.

az, rng = dat["rho"][1].shape
swp = xr.Dataset(dat, coords={"azimuth": np.arange(az), "range": np.arange(rng)})
swp = swp.assign_coords(
    dict(
        longitude=7,
        latitude=53,
        altitude=0,
        elevation=1,
        sweep_mode="azimuth_surveillance",
    )
)
swp = swp.wrl.georef.georeference()
display(swp)

<xarray.Dataset> Size: 6MB
Dimensions:     (azimuth: 360, range: 240)
Coordinates: (12/14)
  * azimuth     (azimuth) int64 3kB 0 1 2 3 4 5 6 ... 354 355 356 357 358 359
  * range       (range) int64 2kB 0 1 2 3 4 5 6 ... 233 234 235 236 237 238 239
    x           (azimuth, range) float64 691kB 0.0 6.122e-17 ... -4.153 -4.17
    y           (azimuth, range) float64 691kB 0.0 0.9998 2.0 ... 237.9 238.9
    z           (azimuth, range) float64 691kB 0.0 0.01745 ... 4.157 4.174
    gr          (azimuth, range) float64 691kB 0.0 0.9998 2.0 ... 238.0 239.0
    ...          ...
    longitude   int64 8B 7
    latitude    int64 8B 53
    altitude    int64 8B 0
    elevation   int64 8B 1
    sweep_mode  <U20 80B 'azimuth_surveillance'
    crs_wkt     int64 8B 0
Data variables:
    rho         (azimuth, range) float32 346kB 0.8743 0.751 ... 0.2174 0.2788
    phi         (azimuth, range) float32 346kB 50.12 50.12 50.12 ... nan nan nan
    ref         (azimuth, range) float32 346kB -15.5 -18.5 -12.0 ... nan nan nan
    dop         (azimuth, range) float32 346kB -0.4297 -0.4297 ... nan nan
    zdr         (azimuth, range) float32 346kB 8.0 -7.059 -0.1882 ... nan nan
    map         (azimuth, range) bool 86kB True True True ... False False False

Identify non-meteorological echoes#

By defining weights for the used moments we can apply the final classifier. The algorithm returns the probability of meteorological echo.

See [Crisologo et al., 2014] and [Vulpiani et al., 2012] for details.

moments = dict(rho="rho", phi="phi", dop="dop", zdr="zdr", map="map")
weights = {"zdr": 0.4, "rho": 0.4, "rho2": 0.4, "phi": 0.1, "dop": 0.1, "map": 0.5}
prob, nanmask = swp.wrl.classify.classify_echo_fuzzy(moments, weights=weights)
thresh = 0.5
cmap = xr.where(prob < thresh, True, False)

Plot classification results#

fig = plt.figure(figsize=(12, 5))

# Horizontal reflectivity
ax = plt.subplot(121, aspect="equal")
pm = swp.ref.plot(x="x", y="y", ax=ax, cbar_kwargs=dict(label="dBZ"))
ax = wradlib.vis.plot_ppi_crosshair(site=(0, 0, 0), ranges=[80, 160, 240])
ax.set_xlim(-240, 240)
ax.set_ylim(-240, 240)
ax.set_xlabel("# bins from radar")
ax.set_ylabel("# bins from radar")
ax.set_title("Reflectivity")

# Echo classification
ax = plt.subplot(122, aspect="equal")
pm = cmap.where(~np.isnan(swp.ref)).plot(
    x="x",
    y="y",
    ax=ax,
    cmap="bwr",
    cbar_kwargs=dict(label="meterol. echo=0 - non-meteorol. echo=1"),
)
ax = wradlib.vis.plot_ppi_crosshair(site=(0, 0, 0), ranges=[80, 160, 240])
ax.set_xlim(-240, 240)
ax.set_ylim(-240, 240)
ax.set_xlabel("# bins from radar")
ax.set_ylabel("# bins from radar")
ax.set_title("Classification")
fig.tight_layout()
../../_images/8357198bd9d34d2b848943c35c5c358f447057ad6d332fae2a3e8c53c9196227.png