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/2.7.0/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
  * range              (range) float32 4kB 50.0 150.0 ... 9.985e+04 9.995e+04
    sweep_mode         <U20 80B 'azimuth_surveillance'
    elevation          (azimuth) float64 3kB 1.505 1.505 1.505 ... 1.505 1.505
    time               (azimuth) datetime64[ns] 3kB ...
    longitude          float64 8B 7.072
    ...                 ...
    y                  (azimuth, range) float64 3MB 49.98 149.9 ... 9.988e+04
    z                  (azimuth, range) float64 3MB 100.8 103.4 ... 3.313e+03
    gr                 (azimuth, range) float64 3MB 49.98 149.9 ... 9.988e+04
    rays               (azimuth, range) float64 3MB 0.5054 0.5054 ... 359.5
    bins               (azimuth, range) float32 1MB 50.0 150.0 ... 9.995e+04
    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
  * range       (range) float32 4kB 50.0 150.0 250.0 ... 9.985e+04 9.995e+04
    sweep_mode  <U20 80B 'azimuth_surveillance'
    elevation   (azimuth) float64 3kB 1.505 1.505 1.505 ... 1.505 1.505 1.505
    time        (azimuth) datetime64[ns] 3kB ...
    longitude   float64 8B 7.072
    ...          ...
    y           (azimuth, range) float64 3MB 49.98 149.9 ... 9.978e+04 9.988e+04
    z           (azimuth, range) float64 3MB 100.8 103.4 ... 3.309e+03 3.313e+03
    gr          (azimuth, range) float64 3MB 49.98 149.9 ... 9.978e+04 9.988e+04
    rays        (azimuth, range) float64 3MB 0.5054 0.5054 ... 359.5 359.5
    bins        (azimuth, range) float32 1MB 50.0 150.0 ... 9.985e+04 9.995e+04
    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/1215a878914a76ce85387eec3c54417c41db2fea933081fd57207c0898db8462.png