Clutter detection by using space-born cloud images

Clutter detection by using space-born cloud images#

[1]:

import numpy as np

import xarray as xr
import xradar as xd
import wradlib as wrl
import matplotlib.pyplot as plt
from osgeo import osr

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()

Read the radar data into DataTree#

[2]:

# read the radar volume scan
filename = "hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf"
filename = wrl.util.get_wradlib_data_file(filename)
pvol = xd.io.open_odim_datatree(filename)
display(pvol)

/home/runner/micromamba/envs/wradlib-notebooks/lib/python3.12/site-packages/wradlib_data/__init__.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  from pkg_resources import DistributionNotFound, get_distribution
Downloading file 'hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf' to '/home/runner/work/wradlib-notebooks/wradlib-notebooks/wradlib-data'.

<xarray.DatasetView> Size: 312B
Dimensions:              (sweep: 5)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 8B 0
    platform_type        <U5 20B 'fixed'
    instrument_type      <U5 20B 'radar'
    time_coverage_start  <U20 80B '2013-04-29T04:30:00Z'
    time_coverage_end    <U20 80B '2013-04-29T04:31:39Z'
    longitude            float64 8B 5.506
    altitude             float64 8B 592.0
    latitude             float64 8B 49.91
    sweep_fixed_angle    (sweep) float64 40B 0.3 0.9 1.8 3.3 6.0
    sweep_group_name     (sweep) int64 40B 0 1 2 3 4
Attributes:
    Conventions:      ODIM_H5/V2_2
    instrument_name:  None
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar

Georeference sweeps#

[3]:

pvol1 = pvol.match("sweep*")
display(pvol1)
vol = []
for sweep in pvol1.values():
    vol.append(sweep.to_dataset().pipe(wrl.georef.georeference))
vol = xr.concat(vol, dim="tilt")
vol = vol.assign_coords(sweep_mode=vol.sweep_mode)
display(vol)

<xarray.DatasetView> Size: 0B
Dimensions:  ()
Data variables:
    *empty*
<xarray.Dataset> Size: 72MB
Dimensions:            (tilt: 5, azimuth: 360, range: 960)
Coordinates: (12/15)
    sweep_mode         (tilt) <U20 400B 'azimuth_surveillance' ... 'azimuth_s...
    elevation          (tilt, azimuth) float64 14kB 0.3 0.3 0.3 ... 6.0 6.0 6.0
    time               (tilt, azimuth) datetime64[ns] 14kB 2013-04-29T04:30:0...
  * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
    longitude          float64 8B 5.506
    latitude           float64 8B 49.91
    ...                 ...
    y                  (tilt, azimuth, range) float64 14MB 125.0 ... 2.378e+05
    z                  (tilt, azimuth, range) float64 14MB 592.7 ... 2.901e+04
    gr                 (tilt, azimuth, range) float64 14MB 125.0 ... 2.378e+05
    rays               (azimuth, range) float32 1MB 0.5 0.5 0.5 ... 359.5 359.5
    bins               (azimuth, range) float32 1MB 125.0 375.0 ... 2.399e+05
    crs_wkt            int64 8B 0
Dimensions without coordinates: tilt
Data variables:
    DBZH               (tilt, azimuth, range) float64 14MB -32.0 -32.0 ... -32.0
    sweep_number       (tilt) int64 40B 0 1 2 3 4
    prt_mode           (tilt) <U7 140B 'not_set' 'not_set' ... 'not_set'
    follow_mode        (tilt) <U7 140B 'not_set' 'not_set' ... 'not_set'
    sweep_fixed_angle  (tilt) float64 40B 0.3 0.9 1.8 3.3 6.0

Construct collocated satellite data#

[4]:

proj_radar = osr.SpatialReference()
proj_radar.ImportFromWkt(vol.crs_wkt.attrs["crs_wkt"])

[4]:

0

[5]:

filename = "hdf5/SAFNWC_MSG3_CT___201304290415_BEL_________.h5"
filename = wrl.util.get_wradlib_data_file(filename)

Downloading file 'hdf5/SAFNWC_MSG3_CT___201304290415_BEL_________.h5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/SAFNWC_MSG3_CT___201304290415_BEL_________.h5' to '/home/runner/work/wradlib-notebooks/wradlib-notebooks/wradlib-data'.

[6]:

sat_gdal = wrl.io.read_safnwc(filename)
val_sat = wrl.georef.read_gdal_values(sat_gdal)
coord_sat = wrl.georef.read_gdal_coordinates(sat_gdal)
proj_sat = wrl.georef.read_gdal_projection(sat_gdal)
coord_sat = wrl.georef.reproject(coord_sat, src_crs=proj_sat, trg_crs=proj_radar)

[7]:

coord_radar = np.stack((vol.x, vol.y), axis=-1)
coord_sat[..., 0:2].reshape(-1, 2).shape, coord_radar[..., 0:2].reshape(-1, 2).shape

[7]:

((180000, 2), (1728000, 2))

[8]:

interp = wrl.ipol.Nearest(
    coord_sat[..., 0:2].reshape(-1, 2), coord_radar[..., 0:2].reshape(-1, 2)
)

[9]:

val_sat = interp(val_sat.ravel()).reshape(coord_radar.shape[:-1])

Estimate localisation errors#

[10]:

timelag = 9 * 60
wind = 10
error = np.absolute(timelag) * wind

Identify clutter based on collocated cloudtype#

[11]:

rscale = vol.range.diff("range").median().values
clutter = wrl.classify.filter_cloudtype(
    vol.DBZH, val_sat, scale=rscale, smoothing=error
)

Assign to vol#

[12]:

vol = vol.assign(sat=(["tilt", "azimuth", "range"], val_sat))
vol = vol.assign(clutter=(["tilt", "azimuth", "range"], clutter.values))
display(vol)

<xarray.Dataset> Size: 75MB
Dimensions:            (tilt: 5, azimuth: 360, range: 960)
Coordinates: (12/15)
    sweep_mode         (tilt) <U20 400B 'azimuth_surveillance' ... 'azimuth_s...
    elevation          (tilt, azimuth) float64 14kB 0.3 0.3 0.3 ... 6.0 6.0 6.0
    time               (tilt, azimuth) datetime64[ns] 14kB 2013-04-29T04:30:0...
  * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
    longitude          float64 8B 5.506
    latitude           float64 8B 49.91
    ...                 ...
    y                  (tilt, azimuth, range) float64 14MB 125.0 ... 2.378e+05
    z                  (tilt, azimuth, range) float64 14MB 592.7 ... 2.901e+04
    gr                 (tilt, azimuth, range) float64 14MB 125.0 ... 2.378e+05
    rays               (azimuth, range) float32 1MB 0.5 0.5 0.5 ... 359.5 359.5
    bins               (azimuth, range) float32 1MB 125.0 375.0 ... 2.399e+05
    crs_wkt            int64 8B 0
Dimensions without coordinates: tilt
Data variables:
    DBZH               (tilt, azimuth, range) float64 14MB -32.0 -32.0 ... -32.0
    sweep_number       (tilt) int64 40B 0 1 2 3 4
    prt_mode           (tilt) <U7 140B 'not_set' 'not_set' ... 'not_set'
    follow_mode        (tilt) <U7 140B 'not_set' 'not_set' ... 'not_set'
    sweep_fixed_angle  (tilt) float64 40B 0.3 0.9 1.8 3.3 6.0
    sat                (tilt, azimuth, range) uint8 2MB 6 6 6 6 ... 10 10 10 10
    clutter            (tilt, azimuth, range) bool 2MB False False ... False

Plot the results#

[13]:

fig = plt.figure(figsize=(16, 8))

tilt = 0

ax = fig.add_subplot(131)
pm = vol.DBZH[tilt].wrl.vis.plot(ax=ax)
# plt.colorbar(pm, shrink=0.5)
plt.title("Radar reflectivity")

ax = fig.add_subplot(132)
pm = vol.sat[tilt].wrl.vis.plot(ax=ax)
# plt.colorbar(pm, shrink=0.5)
plt.title("Satellite cloud classification")

ax = fig.add_subplot(133)
pm = vol.clutter[tilt].wrl.vis.plot(ax=ax)
# plt.colorbar(pm, shrink=0.5)
plt.title("Detected clutter")

fig.tight_layout()
../../_images/notebooks_classify_clutter_cloud_21_0.png