Clutter detection using the Gabella approach

Clutter detection using the Gabella approach#

import warnings

import matplotlib.pyplot as plt
import numpy as np
import wradlib as wrl
import wradlib_data
import xarray as xr

warnings.filterwarnings("ignore")
/home/docs/checkouts/readthedocs.org/user_builds/wradlib/conda/latest/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Read the data#

filename = wradlib_data.DATASETS.fetch("hdf5/2014-08-10--182000.ppi.mvol")
swp = xr.open_dataset(filename, engine="gamic", group="sweep_0")
swp = swp.set_coords("sweep_mode")
swp = swp.wrl.georef.georeference()
display(swp)
Downloading file 'hdf5/2014-08-10--182000.ppi.mvol' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/2014-08-10--182000.ppi.mvol' to '/home/docs/.cache/wradlib-data'.
<xarray.Dataset> Size: 33MB
Dimensions:            (azimuth: 360, range: 1000)
Coordinates: (12/15)
  * azimuth            (azimuth) float64 3kB 0.5054 1.522 2.527 ... 358.5 359.5
    elevation          (azimuth) float64 3kB 1.505 1.505 1.505 ... 1.505 1.505
    time               (azimuth) datetime64[ns] 3kB ...
  * range              (range) float32 4kB 50.0 150.0 ... 9.985e+04 9.995e+04
    x                  (azimuth, range) float64 3MB 0.4409 1.323 ... -866.6
    y                  (azimuth, range) float64 3MB 49.98 149.9 ... 9.988e+04
    ...                 ...
    bins               (azimuth, range) float32 1MB 50.0 150.0 ... 9.995e+04
    sweep_mode         <U20 80B 'azimuth_surveillance'
    longitude          float64 8B 7.072
    latitude           float64 8B 50.73
    altitude           float64 8B 99.5
    crs_wkt            int64 8B 0
Data variables: (12/17)
    KDP                (azimuth, range) float32 1MB ...
    PHIDP              (azimuth, range) float32 1MB ...
    DBZH               (azimuth, range) float32 1MB ...
    DBZV               (azimuth, range) float32 1MB ...
    RHOHV              (azimuth, range) float32 1MB ...
    DBTH               (azimuth, range) float32 1MB ...
    ...                 ...
    ZDR                (azimuth, range) float32 1MB ...
    sweep_number       int64 8B ...
    prt_mode           <U7 28B ...
    follow_mode        <U7 28B ...
    sweep_fixed_angle  float64 8B ...
    nyquist_velocity   object 8B ...
Attributes:
    source:       gamic
    pulse_width:  2

Apply filter#

clmap = swp.DBTH.wrl.classify.filter_gabella(
    wsize=5, thrsnorain=0.0, tr1=6.0, n_p=8, tr2=1.3
)
clmap
<xarray.DataArray 'filter_gabella' (azimuth: 360, range: 1000)> Size: 360kB
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]], shape=(360, 1000))
Coordinates: (12/15)
  * azimuth     (azimuth) float64 3kB 0.5054 1.522 2.527 ... 357.5 358.5 359.5
    elevation   (azimuth) float64 3kB 1.505 1.505 1.505 ... 1.505 1.505 1.505
    time        (azimuth) datetime64[ns] 3kB ...
  * range       (range) float32 4kB 50.0 150.0 250.0 ... 9.985e+04 9.995e+04
    x           (azimuth, range) float64 3MB 0.4409 1.323 ... -865.7 -866.6
    y           (azimuth, range) float64 3MB 49.98 149.9 ... 9.978e+04 9.988e+04
    ...          ...
    bins        (azimuth, range) float32 1MB 50.0 150.0 ... 9.985e+04 9.995e+04
    sweep_mode  <U20 80B 'azimuth_surveillance'
    longitude   float64 8B 7.072
    latitude    float64 8B 50.73
    altitude    float64 8B 99.5
    crs_wkt     int64 8B 0
Attributes:
    standard_name:  clutter_map
    long_name:      Clutter Map
    short_name:     CMAP
    units:          unitless

Plot results#

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(121)
pm = swp.DBTH.wrl.vis.plot(ax=ax1)
ax1.set_title("Reflectivity")
ax2 = fig.add_subplot(122)
pm = clmap.wrl.vis.plot(ax=ax2)
ax2.set_title("Cluttermap")
fig.tight_layout()
../../../_images/284311b8dad6df670cec890bb19647405bd8d19086f81a1e0299a6797ab7cdfc.png