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'.
---------------------------------------------------------------------------
RemoteDisconnected Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
788 conn,
789 method,
790 url,
791 timeout=timeout_obj,
792 body=body,
793 headers=headers,
794 chunked=chunked,
795 retries=retries,
796 response_conn=response_conn,
797 preload_content=preload_content,
798 decode_content=decode_content,
799 **response_kw,
800 )
802 # Everything went great!
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
533 try:
--> 534 response = conn.getresponse()
535 except (BaseSSLError, OSError) as e:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connection.py:571, in HTTPConnection.getresponse(self)
570 # Get the response from http.client.HTTPConnection
--> 571 httplib_response = super().getresponse()
573 try:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:1450, in HTTPConnection.getresponse(self)
1449 try:
-> 1450 response.begin()
1451 except ConnectionError:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:336, in HTTPResponse.begin(self)
335 while True:
--> 336 version, status, reason = self._read_status()
337 if status != CONTINUE:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:305, in HTTPResponse._read_status(self)
302 if not line:
303 # Presumably, the server closed the connection before
304 # sending a valid response.
--> 305 raise RemoteDisconnected("Remote end closed connection without"
306 " response")
307 try:
RemoteDisconnected: Remote end closed connection without response
During handling of the above exception, another exception occurred:
ProtocolError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/adapters.py:644, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
643 try:
--> 644 resp = conn.urlopen(
645 method=request.method,
646 url=url,
647 body=request.body,
648 headers=request.headers,
649 redirect=False,
650 assert_same_host=False,
651 preload_content=False,
652 decode_content=False,
653 retries=self.max_retries,
654 timeout=timeout,
655 chunked=chunked,
656 )
658 except (ProtocolError, OSError) as err:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connectionpool.py:841, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
839 new_e = ProtocolError("Connection aborted.", new_e)
--> 841 retries = retries.increment(
842 method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
843 )
844 retries.sleep()
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/util/retry.py:490, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
489 if read is False or method is None or not self._is_method_retryable(method):
--> 490 raise reraise(type(error), error, _stacktrace)
491 elif read is not None:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/util/util.py:38, in reraise(tp, value, tb)
37 if value.__traceback__ is not tb:
---> 38 raise value.with_traceback(tb)
39 raise value
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
788 conn,
789 method,
790 url,
791 timeout=timeout_obj,
792 body=body,
793 headers=headers,
794 chunked=chunked,
795 retries=retries,
796 response_conn=response_conn,
797 preload_content=preload_content,
798 decode_content=decode_content,
799 **response_kw,
800 )
802 # Everything went great!
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
533 try:
--> 534 response = conn.getresponse()
535 except (BaseSSLError, OSError) as e:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/urllib3/connection.py:571, in HTTPConnection.getresponse(self)
570 # Get the response from http.client.HTTPConnection
--> 571 httplib_response = super().getresponse()
573 try:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:1450, in HTTPConnection.getresponse(self)
1449 try:
-> 1450 response.begin()
1451 except ConnectionError:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:336, in HTTPResponse.begin(self)
335 while True:
--> 336 version, status, reason = self._read_status()
337 if status != CONTINUE:
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/http/client.py:305, in HTTPResponse._read_status(self)
302 if not line:
303 # Presumably, the server closed the connection before
304 # sending a valid response.
--> 305 raise RemoteDisconnected("Remote end closed connection without"
306 " response")
307 try:
ProtocolError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))
During handling of the above exception, another exception occurred:
ConnectionError Traceback (most recent call last)
Cell In[2], line 6
4 dopfile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-V.nc")
5 zdrfile = get_wradlib_data_file("netcdf/TAG-20120801-140046-02-D.nc")
----> 6 mapfile = get_wradlib_data_file("hdf5/TAG_cmap_sweeps_0204050607.hdf5")
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/wradlib/util.py:703, in get_wradlib_data_file(relfile)
696 def get_wradlib_data_file(relfile):
697 warn(
698 "Function get_wradlib_data_file is not part of public API,\n"
699 "it's use is deprecated and it will be removed in a future version.\n"
700 "Please see wradlib-data package for more information.",
701 DeprecationWarning,
702 )
--> 703 return _get_wradlib_data_file(relfile)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/wradlib/util.py:693, in _get_wradlib_data_file(relfile)
691 def _get_wradlib_data_file(relfile):
692 wradlib_data = import_optional("wradlib_data", dep="development")
--> 693 return wradlib_data.DATASETS.fetch(relfile)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/wradlib_data/dataset.py:48, in fetch(*args, **kwargs)
45 @wraps(DATASETS._fetch)
46 def fetch(*args, **kwargs):
47 kwargs.setdefault('downloader', wradlib_downloader)
---> 48 return DATASETS._fetch(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/pooch/core.py:598, in Pooch.fetch(self, fname, processor, downloader, progressbar)
595 if downloader is None:
596 downloader = choose_downloader(url, progressbar=progressbar)
--> 598 stream_download(
599 url,
600 full_path,
601 known_hash,
602 downloader,
603 pooch=self,
604 retry_if_failed=self.retry_if_failed,
605 )
607 if processor is not None:
608 return processor(str(full_path), action, self)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/pooch/core.py:823, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
819 try:
820 # Stream the file to a temporary so that we can safely check its
821 # hash before overwriting the original.
822 with temporary_file(path=str(fname.parent)) as tmp:
--> 823 downloader(url, tmp, pooch)
824 hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
825 shutil.move(tmp, str(fname))
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/wradlib_data/dataset.py:37, in wradlib_downloader(url, output_file, mypooch)
35 headers = {'User-Agent': f'wradlib-data {wradlib_data.__version__}'}
36 https = pooch.HTTPDownloader(headers=headers)
---> 37 https(url, output_file, mypooch)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/pooch/downloaders.py:230, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
228 # pylint: enable=consider-using-with
229 try:
--> 230 response = requests.get(url, timeout=timeout, **kwargs)
231 response.raise_for_status()
232 content = response.iter_content(chunk_size=self.chunk_size)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/api.py:73, in get(url, params, **kwargs)
62 def get(url, params=None, **kwargs):
63 r"""Sends a GET request.
64
65 :param url: URL for the new :class:`Request` object.
(...) 70 :rtype: requests.Response
71 """
---> 73 return request("get", url, params=params, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/api.py:59, in request(method, url, **kwargs)
55 # By using the 'with' statement we are sure the session is closed, thus we
56 # avoid leaving sockets open which can trigger a ResourceWarning in some
57 # cases, and look like a memory leak in others.
58 with sessions.Session() as session:
---> 59 return session.request(method=method, url=url, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
584 send_kwargs = {
585 "timeout": timeout,
586 "allow_redirects": allow_redirects,
587 }
588 send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
591 return resp
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/sessions.py:703, in Session.send(self, request, **kwargs)
700 start = preferred_clock()
702 # Send the request
--> 703 r = adapter.send(request, **kwargs)
705 # Total elapsed time of the request (approximately)
706 elapsed = preferred_clock() - start
File ~/checkouts/readthedocs.org/user_builds/wradlib/conda/2.7.0/lib/python3.13/site-packages/requests/adapters.py:659, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
644 resp = conn.urlopen(
645 method=request.method,
646 url=url,
(...) 655 chunked=chunked,
656 )
658 except (ProtocolError, OSError) as err:
--> 659 raise ConnectionError(err, request=request)
661 except MaxRetryError as e:
662 if isinstance(e.reason, ConnectTimeoutError):
663 # TODO: Remove this in 3.0.0: see #2811
ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))
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)
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()