.ipynb .pdf

repository open issue suggest edit

On this page

wradlib data quality

---

Overview

Within this notebook, we will cover:

1. Reading radar volume data into xarray based RadarVolume

2. Wrapping numpy-based functions to work with Xarray

3. Clutter detection

4. Beam Blockage calculation

Prerequisites

Concepts	Importance	Notes
Xarray Basics	Helpful	Basic Dataset/DataArray
Matplotlib Basics	Helpful	Basic Plotting
Intro to Cartopy	Helpful	Projections

• Time to learn: 10 minutes

---

Imports

import glob
import os

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import wradlib as wrl

/srv/conda/envs/notebook/lib/python3.9/site-packages/requests/__init__.py:102: RequestsDependencyWarning: urllib3 (1.26.8) or chardet (5.2.0)/charset_normalizer (2.0.10) doesn't match a supported version!
  warnings.warn("urllib3 ({}) or chardet ({})/charset_normalizer ({}) doesn't match a supported "

Import Swiss Radar Data from CfRadial1 Volumes

We use some of the pyrad example data here. Sweeps are provided as single files, so we open each file separately and create the RadarVolume from the open Datasets.

fglob = "../pyart/data/example_pyrad/22179/MLL22179/MLL2217907250U*.nc"
flist = glob.glob(fglob)
flist.sort()
print("Files available: {}".format(len(flist)))

Files available: 20

ds = [xr.open_dataset(f, group="sweep_1", engine="cfradial1", chunks={}) for f in flist]

vol = wrl.io.RadarVolume(engine="cfradial1")
vol.extend(ds)
vol.sort(key=lambda x: x.time.min())

display(vol)

<wradlib.RadarVolume>
Dimension(s): (sweep: 20)
Elevation(s): (2.5, 4.5, 6.5, -0.2, 1.0, 7.5, 8.5, 11.0, 16.0, 25.0, 35.0, 0.4, 1.6, 3.5, 5.5, 9.5, 13.0, 20.0, 30.0, 40.0)

vol.root

<xarray.Dataset>
Dimensions:              (sweep: 20)
Coordinates:
    sweep_mode           <U20 'azimuth_surveillance'
    longitude            float32 8.833
    altitude             float32 1.626e+03
    time                 datetime64[ns] 2022-06-28T07:19:28
    latitude             float32 46.04
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 0
    platform_type        <U5 'fixed'
    instrument_type      <U5 'radar'
    primary_axis         <U6 'axis_z'
    time_coverage_start  <U20 '2022-06-28T07:19:28Z'
    time_coverage_end    <U20 '2022-06-28T07:23:44Z'
    sweep_group_name     (sweep) <U8 'sweep_0' 'sweep_1' ... 'sweep_19'
    sweep_fixed_angle    (sweep) float32 2.5 4.5 6.5 -0.2 ... 20.0 30.0 40.0
Attributes:
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using wradlib
    instrument_name:  None
    fixed_angle:      2.5

xarray.Dataset

• Dimensions:
  • sweep: 20
• Coordinates: (5)
  • sweep_mode
    ()
    <U20
    'azimuth_surveillance'

    array('azimuth_surveillance', dtype='<U20')

  • longitude
    ()
    float32
    8.833

    long_name :

    Longitude

    units :

    degrees_east

    standard_name :

    Longitude

    array(8.833217, dtype=float32)

  • altitude
    ()
    float32
    1.626e+03

    long_name :

    Altitude

    units :

    meters

    standard_name :

    Altitude

    positive :

    up

    array(1626., dtype=float32)

  • time
    ()
    datetime64[ns]
    2022-06-28T07:19:28

    long_name :

    time_in_seconds_since_volume_start

    standard_name :

    time

    array('2022-06-28T07:19:28.000000000', dtype='datetime64[ns]')

  • latitude
    ()
    float32
    46.04

    long_name :

    Latitude

    units :

    degrees_north

    standard_name :

    Latitude

    array(46.04076, dtype=float32)

• Data variables: (8)
  • volume_number
    ()
    int64
    0

    array(0)

  • platform_type
    ()
    <U5
    'fixed'

    array('fixed', dtype='<U5')

  • instrument_type
    ()
    <U5
    'radar'

    array('radar', dtype='<U5')

  • primary_axis
    ()
    <U6
    'axis_z'

    array('axis_z', dtype='<U6')

  • time_coverage_start
    ()
    <U20
    '2022-06-28T07:19:28Z'

    array('2022-06-28T07:19:28Z', dtype='<U20')

  • time_coverage_end
    ()
    <U20
    '2022-06-28T07:23:44Z'

    array('2022-06-28T07:23:44Z', dtype='<U20')

  • sweep_group_name
    (sweep)
    <U8
    'sweep_0' 'sweep_1' ... 'sweep_19'

    array(['sweep_0', 'sweep_1', 'sweep_2', 'sweep_3', 'sweep_4', 'sweep_5',
           'sweep_6', 'sweep_7', 'sweep_8', 'sweep_9', 'sweep_10', 'sweep_11',
           'sweep_12', 'sweep_13', 'sweep_14', 'sweep_15', 'sweep_16',
           'sweep_17', 'sweep_18', 'sweep_19'], dtype='<U8')

  • sweep_fixed_angle
    (sweep)
    float32
    2.5 4.5 6.5 -0.2 ... 20.0 30.0 40.0

    array([ 2.5,  4.5,  6.5, -0.2,  1. ,  7.5,  8.5, 11. , 16. , 25. , 35. ,
            0.4,  1.6,  3.5,  5.5,  9.5, 13. , 20. , 30. , 40. ],
          dtype=float32)

• Attributes: (9)

  version :

  None

  title :

  None

  institution :

  None

  references :

  None

  source :

  None

  history :

  None

  comment :

  im/exported using wradlib

  instrument_name :

  None

  fixed_angle :

  2.5

sweep_number = 3
display(vol[sweep_number])

<xarray.Dataset>
Dimensions:                              (azimuth: 360, range: 492, frequency: 1)
Coordinates:
    rtime                                (azimuth) datetime64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
  * range                                (range) float32 250.0 ... 2.457e+05
  * azimuth                              (azimuth) float32 0.5246 ... 359.5
    elevation                            (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    sweep_mode                           <U20 ...
  * frequency                            (frequency) float32 5.451e+09
    latitude                             float32 ...
    longitude                            float32 ...
    altitude                             float32 ...
    time                                 datetime64[ns] ...
Data variables: (12/20)
    reflectivity                         (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    signal_to_noise_ratio                (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    reflectivity_vv                      (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    differential_reflectivity            (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    uncorrected_cross_correlation_ratio  (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    uncorrected_differential_phase       (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    ...                                   ...
    radar_beam_width_v                   float32 ...
    pulse_width                          (azimuth) timedelta64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
    nyquist_velocity                     (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    time_coverage_start                  |S32 ...
    time_coverage_end                    |S32 ...
    volume_number                        int32 ...
Attributes:
    fixed_angle:  -0.2

xarray.Dataset

• Dimensions:
  • azimuth: 360
  • range: 492
  • frequency: 1
• Coordinates: (10)
  • rtime
    (azimuth)
    datetime64[ns]
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    time_in_seconds_since_volume_start

    standard_name :

    time

    comment :

    Coordinate variable for time. Time at the center of each ray, in fractional seconds since the global variable time_coverage_start

    Array Chunk Bytes 2.81 kiB 2.81 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type datetime64[ns] numpy.ndarray

  • range
    (range)
    float32
    250.0 750.0 ... 2.452e+05 2.457e+05

    long_name :

    range_to_measurement_volume

    units :

    meters

    standard_name :

    projection_range_coordinate

    axis :

    radial_range_coordinate

    spacing_is_constant :

    true

    comment :

    Coordinate variable for range. Range to center of each bin.

    array([   249.99901,    749.997  ,   1249.995  , ..., 244749.02   ,
           245249.02   , 245749.02   ], dtype=float32)

  • azimuth
    (azimuth)
    float32
    0.5246 1.527 2.524 ... 358.5 359.5

    long_name :

    azimuth_angle_from_true_north

    units :

    degrees

    standard_name :

    beam_azimuth_angle

    axis :

    radial_azimuth_coordinate

    comment :

    Azimuth of antenna relative to true north

    array([  0.524605,   1.527123,   2.524148, ..., 357.52805 , 358.5278  ,
           359.5276  ], dtype=float32)

  • elevation
    (azimuth)
    float32
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    elevation_angle_from_horizontal_plane

    units :

    degrees

    standard_name :

    beam_elevation_angle

    axis :

    radial_elevation_coordinate

    comment :

    Elevation of antenna relative to the horizontal plane

    Array Chunk Bytes 1.41 kiB 1.41 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • sweep_mode
    ()
    <U20
    ...

    array('azimuth_surveillance', dtype='<U20')

  • frequency
    (frequency)
    float32
    5.451e+09

    long_name :

    Radiation frequency

    units :

    s-1

    meta_group :

    instrument_parameters

    array([5.450772e+09], dtype=float32)

  • latitude
    ()
    float32
    ...

    long_name :

    Latitude

    units :

    degrees_north

    standard_name :

    Latitude

    array(46.04076, dtype=float32)

  • longitude
    ()
    float32
    ...

    long_name :

    Longitude

    units :

    degrees_east

    standard_name :

    Longitude

    array(8.833217, dtype=float32)

  • altitude
    ()
    float32
    ...

    long_name :

    Altitude

    units :

    meters

    standard_name :

    Altitude

    positive :

    up

    array(1626., dtype=float32)

  • time
    ()
    datetime64[ns]
    ...

    long_name :

    time_in_seconds_since_volume_start

    standard_name :

    time

    array('2022-06-28T07:21:36.000000000', dtype='datetime64[ns]')

• Data variables: (20)
  • reflectivity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Reflectivity

    units :

    dBZ

    standard_name :

    equivalent_reflectivity_factor

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • signal_to_noise_ratio
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Signal to noise ratio

    units :

    dB

    standard_name :

    signal_to_noise_ratio

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • reflectivity_vv
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • differential_reflectivity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Differential reflectivity

    units :

    dB

    standard_name :

    log_differential_reflectivity_hv

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • uncorrected_cross_correlation_ratio
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • uncorrected_differential_phase
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • velocity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Mean dopper velocity

    units :

    meters_per_second

    standard_name :

    radial_velocity_of_scatterers_away_from_instrument

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • spectrum_width
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Doppler spectrum width

    units :

    meters_per_second

    standard_name :

    doppler_spectrum_width

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • reflectivity_hh_clut
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • sweep_number
    ()
    int64
    ...

    long_name :

    Sweep number

    units :

    count

    standard_name :

    sweep_number

    array(0)

  • fixed_angle
    ()
    float32
    ...

    long_name :

    Target angle for sweep

    units :

    degrees

    standard_name :

    target_fixed_angle

    array(-0.197757, dtype=float32)

  • rays_are_indexed
    ()
    |S32
    ...

    long_name :

    Flag for indexed rays

    units :

    unitless

    options :

    true: rays are indexed to a regular grid, false: rays are not indexed to a regular grid

    array(b'true', dtype='|S32')

  • ray_angle_res
    ()
    float32
    ...

    long_name :

    Angular resolution between rays

    units :

    degrees

    comment :

    Only applicable when rays_are_indexed variable is true

    array(1., dtype=float32)

  • radar_beam_width_h
    ()
    float32
    ...

    long_name :

    Antenna beam width H polarization

    units :

    degrees

    meta_group :

    radar_parameters

    array(1., dtype=float32)

  • radar_beam_width_v
    ()
    float32
    ...

    long_name :

    Antenna beam width V polarization

    units :

    degrees

    meta_group :

    radar_parameters

    array(1., dtype=float32)

  • pulse_width
    (azimuth)
    timedelta64[ns]
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    Pulse width

    comments :

    Pulse width

    meta_group :

    instrument_parameters

    Array Chunk Bytes 2.81 kiB 2.81 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type timedelta64[ns] numpy.ndarray

  • nyquist_velocity
    (azimuth)
    float32
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    Nyquist velocity

    units :

    meters_per_second

    comments :

    Unambiguous velocity

    meta_group :

    instrument_parameters

    Array Chunk Bytes 1.41 kiB 1.41 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • time_coverage_start
    ()
    |S32
    ...

    long_name :

    UTC time of first ray in the file

    units :

    unitless

    array(b'2022-06-28T07:21:36Z', dtype='|S32')

  • time_coverage_end
    ()
    |S32
    ...

    long_name :

    UTC time of last ray in the file

    units :

    unitless

    array(b'2022-06-28T07:21:36Z', dtype='|S32')

  • volume_number
    ()
    int32
    ...

    long_name :

    Volume number

    units :

    unitless

    array(0, dtype=int32)

• Attributes: (1)

  fixed_angle :

  -0.2

Clutter detection with Gabella

While in Switzerland, why not use the well-known clutter detection scheme by Marco Gabella et. al.

Wrap Gabella Clutter detection in Xarray apply_ufunc

The routine is implemented in wradlib in pure Numpy. Numpy based processing routines can be transformed to a first class Xarray citizen with the help of xr.apply_ufunc.

def extract_clutter(da, wsize=5, thrsnorain=0, tr1=6.0, n_p=6, tr2=1.3, rm_nans=False):
    return xr.apply_ufunc(
        wrl.clutter.filter_gabella,
        da,
        input_core_dims=[["azimuth", "range"]],
        output_core_dims=[["azimuth", "range"]],
        dask="parallelized",
        kwargs=dict(
            wsize=wsize,
            thrsnorain=thrsnorain,
            tr1=tr1,
            n_p=n_p,
            tr2=tr2,
            rm_nans=rm_nans,
        ),
    )

Calculate clutter map

Now we apply Gabella scheme and add the result to the Dataset.

swp = vol[sweep_number]
clmap = swp.reflectivity_hh_clut.pipe(
    extract_clutter, wsize=5, thrsnorain=0.0, tr1=21.0, n_p=23, tr2=1.3, rm_nans=False
)
swp = swp.assign({"CMAP": clmap})
display(swp)

<xarray.Dataset>
Dimensions:                              (azimuth: 360, range: 492, frequency: 1)
Coordinates:
    rtime                                (azimuth) datetime64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
  * range                                (range) float32 250.0 ... 2.457e+05
  * azimuth                              (azimuth) float32 0.5246 ... 359.5
    elevation                            (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    sweep_mode                           <U20 ...
  * frequency                            (frequency) float32 5.451e+09
    latitude                             float32 ...
    longitude                            float32 ...
    altitude                             float32 ...
    time                                 datetime64[ns] ...
Data variables: (12/21)
    reflectivity                         (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    signal_to_noise_ratio                (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    reflectivity_vv                      (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    differential_reflectivity            (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    uncorrected_cross_correlation_ratio  (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    uncorrected_differential_phase       (azimuth, range) float32 dask.array<chunksize=(360, 492), meta=np.ndarray>
    ...                                   ...
    pulse_width                          (azimuth) timedelta64[ns] dask.array<chunksize=(360,), meta=np.ndarray>
    nyquist_velocity                     (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    time_coverage_start                  |S32 ...
    time_coverage_end                    |S32 ...
    volume_number                        int32 ...
    CMAP                                 (azimuth, range) bool dask.array<chunksize=(360, 492), meta=np.ndarray>
Attributes:
    fixed_angle:  -0.2

xarray.Dataset

• Dimensions:
  • azimuth: 360
  • range: 492
  • frequency: 1
• Coordinates: (10)
  • rtime
    (azimuth)
    datetime64[ns]
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    time_in_seconds_since_volume_start

    standard_name :

    time

    comment :

    Coordinate variable for time. Time at the center of each ray, in fractional seconds since the global variable time_coverage_start

    Array Chunk Bytes 2.81 kiB 2.81 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type datetime64[ns] numpy.ndarray

  • range
    (range)
    float32
    250.0 750.0 ... 2.452e+05 2.457e+05

    long_name :

    range_to_measurement_volume

    units :

    meters

    standard_name :

    projection_range_coordinate

    axis :

    radial_range_coordinate

    spacing_is_constant :

    true

    comment :

    Coordinate variable for range. Range to center of each bin.

    array([   249.99901,    749.997  ,   1249.995  , ..., 244749.02   ,
           245249.02   , 245749.02   ], dtype=float32)

  • azimuth
    (azimuth)
    float32
    0.5246 1.527 2.524 ... 358.5 359.5

    long_name :

    azimuth_angle_from_true_north

    units :

    degrees

    standard_name :

    beam_azimuth_angle

    axis :

    radial_azimuth_coordinate

    comment :

    Azimuth of antenna relative to true north

    array([  0.524605,   1.527123,   2.524148, ..., 357.52805 , 358.5278  ,
           359.5276  ], dtype=float32)

  • elevation
    (azimuth)
    float32
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    elevation_angle_from_horizontal_plane

    units :

    degrees

    standard_name :

    beam_elevation_angle

    axis :

    radial_elevation_coordinate

    comment :

    Elevation of antenna relative to the horizontal plane

    Array Chunk Bytes 1.41 kiB 1.41 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • sweep_mode
    ()
    <U20
    ...

    array('azimuth_surveillance', dtype='<U20')

  • frequency
    (frequency)
    float32
    5.451e+09

    long_name :

    Radiation frequency

    units :

    s-1

    meta_group :

    instrument_parameters

    array([5.450772e+09], dtype=float32)

  • latitude
    ()
    float32
    ...

    long_name :

    Latitude

    units :

    degrees_north

    standard_name :

    Latitude

    array(46.04076, dtype=float32)

  • longitude
    ()
    float32
    ...

    long_name :

    Longitude

    units :

    degrees_east

    standard_name :

    Longitude

    array(8.833217, dtype=float32)

  • altitude
    ()
    float32
    ...

    long_name :

    Altitude

    units :

    meters

    standard_name :

    Altitude

    positive :

    up

    array(1626., dtype=float32)

  • time
    ()
    datetime64[ns]
    ...

    long_name :

    time_in_seconds_since_volume_start

    standard_name :

    time

    array('2022-06-28T07:21:36.000000000', dtype='datetime64[ns]')

• Data variables: (21)
  • reflectivity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Reflectivity

    units :

    dBZ

    standard_name :

    equivalent_reflectivity_factor

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • signal_to_noise_ratio
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Signal to noise ratio

    units :

    dB

    standard_name :

    signal_to_noise_ratio

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • reflectivity_vv
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • differential_reflectivity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Differential reflectivity

    units :

    dB

    standard_name :

    log_differential_reflectivity_hv

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • uncorrected_cross_correlation_ratio
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • uncorrected_differential_phase
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • velocity
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Mean dopper velocity

    units :

    meters_per_second

    standard_name :

    radial_velocity_of_scatterers_away_from_instrument

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • spectrum_width
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    long_name :

    Doppler spectrum width

    units :

    meters_per_second

    standard_name :

    doppler_spectrum_width

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • reflectivity_hh_clut
    (azimuth, range)
    float32
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 691.88 kiB 691.88 kiB Shape (360, 492) (360, 492) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • sweep_number
    ()
    int64
    ...

    long_name :

    Sweep number

    units :

    count

    standard_name :

    sweep_number

    array(0)

  • fixed_angle
    ()
    float32
    ...

    long_name :

    Target angle for sweep

    units :

    degrees

    standard_name :

    target_fixed_angle

    array(-0.197757, dtype=float32)

  • rays_are_indexed
    ()
    |S32
    ...

    long_name :

    Flag for indexed rays

    units :

    unitless

    options :

    true: rays are indexed to a regular grid, false: rays are not indexed to a regular grid

    array(b'true', dtype='|S32')

  • ray_angle_res
    ()
    float32
    ...

    long_name :

    Angular resolution between rays

    units :

    degrees

    comment :

    Only applicable when rays_are_indexed variable is true

    array(1., dtype=float32)

  • radar_beam_width_h
    ()
    float32
    ...

    long_name :

    Antenna beam width H polarization

    units :

    degrees

    meta_group :

    radar_parameters

    array(1., dtype=float32)

  • radar_beam_width_v
    ()
    float32
    ...

    long_name :

    Antenna beam width V polarization

    units :

    degrees

    meta_group :

    radar_parameters

    array(1., dtype=float32)

  • pulse_width
    (azimuth)
    timedelta64[ns]
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    Pulse width

    comments :

    Pulse width

    meta_group :

    instrument_parameters

    Array Chunk Bytes 2.81 kiB 2.81 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type timedelta64[ns] numpy.ndarray

  • nyquist_velocity
    (azimuth)
    float32
    dask.array<chunksize=(360,), meta=np.ndarray>

    long_name :

    Nyquist velocity

    units :

    meters_per_second

    comments :

    Unambiguous velocity

    meta_group :

    instrument_parameters

    Array Chunk Bytes 1.41 kiB 1.41 kiB Shape (360,) (360,) Count 2 Graph Layers 1 Chunks Type float32 numpy.ndarray

  • time_coverage_start
    ()
    |S32
    ...

    long_name :

    UTC time of first ray in the file

    units :

    unitless

    array(b'2022-06-28T07:21:36Z', dtype='|S32')

  • time_coverage_end
    ()
    |S32
    ...

    long_name :

    UTC time of last ray in the file

    units :

    unitless

    array(b'2022-06-28T07:21:36Z', dtype='|S32')

  • volume_number
    ()
    int32
    ...

    long_name :

    Volume number

    units :

    unitless

    array(0, dtype=int32)

  • CMAP
    (azimuth, range)
    bool
    dask.array<chunksize=(360, 492), meta=np.ndarray>

    Array Chunk Bytes 172.97 kiB 172.97 kiB Shape (360, 492) (360, 492) Count 5 Graph Layers 1 Chunks Type bool numpy.ndarray

• Attributes: (1)

  fixed_angle :

  -0.2

Plot Reflectivities, Clutter and Cluttermap

from osgeo import osr

wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)

0

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

swpx = swp.sel(range=slice(0, 100000)).pipe(wrl.georef.georeference_dataset, proj=wgs84)

ax1 = fig.add_subplot(221)
swpx.reflectivity_hh_clut.plot(x="x", y="y", ax=ax1, vmin=0, vmax=60)
ax1.set_title("Reflectivity raw")

ax2 = fig.add_subplot(222)
swpx.CMAP.plot(x="x", y="y", ax=ax2)
ax2.set_title("Cluttermap")

ax3 = fig.add_subplot(223)
swpx.reflectivity_hh_clut.where(swpx.CMAP == 1).plot(
    x="x", y="y", ax=ax3, vmin=0, vmax=60
)
ax3.set_title("Clutter")

ax4 = fig.add_subplot(224)
swpx.reflectivity_hh_clut.where(swpx.CMAP < 1).plot(
    x="x", y="y", ax=ax4, vmin=0, vmax=60
)
ax4.set_title("Reflectivity clutter removed")

Text(0.5, 1.0, 'Reflectivity clutter removed')

SRTM based clutter and beamblockage processing

Download needed SRTM data

For the course we already provide the needed SRTM tiles. For normal operation you would need a NASA EARTHDATA account and a connected bearer token.

The data will be loaded using GDAL machinery and transformed into an Xarray DataArray.

extent = wrl.zonalstats.get_bbox(swpx.x.values, swpx.y.values)
extent

{'left': 7.544804108214879,
 'right': 10.121619417254358,
 'bottom': 45.14348398510345,
 'top': 46.9378914870384}

# apply fake token, data is already available
os.environ["WRADLIB_EARTHDATA_BEARER_TOKEN"] = ""
# set location of wradlib-data, where wradlib will search for any available data
os.environ["WRADLIB_DATA"] = "data/wradlib-data"
# get the tiles
dem = wrl.io.get_srtm(extent.values())

elevation = wrl.georef.read_gdal_values(dem)
coords = wrl.georef.read_gdal_coordinates(dem)

elev = xr.DataArray(
    data=elevation,
    dims=["y", "x"],
    coords={"lat": (["y", "x"], coords[..., 1]), "lon": (["y", "x"], coords[..., 0])},
)
elev

<xarray.DataArray (y: 2401, x: 4801)>
array([[ 422,  422,  422, ..., 1846, 1798, 1745],
       [ 422,  422,  422, ..., 1856, 1807, 1748],
       [ 422,  422,  422, ..., 1815, 1781, 1740],
       ...,
       [2403, 2379, 2357, ...,   14,   15,   15],
       [2408, 2387, 2365, ...,   14,   15,   14],
       [2428, 2411, 2397, ...,   13,   13,   13]], dtype=int16)
Coordinates:
    lat      (y, x) float64 47.0 47.0 47.0 47.0 47.0 ... 45.0 45.0 45.0 45.0
    lon      (y, x) float64 7.0 7.001 7.002 7.003 7.003 ... 11.0 11.0 11.0 11.0
Dimensions without coordinates: y, x

xarray.DataArray

• y: 2401
• x: 4801

• 422 422 422 422 422 422 422 422 422 422 ... 14 14 14 13 13 13 13 13 13

  array([[ 422,  422,  422, ..., 1846, 1798, 1745],
         [ 422,  422,  422, ..., 1856, 1807, 1748],
         [ 422,  422,  422, ..., 1815, 1781, 1740],
         ...,
         [2403, 2379, 2357, ...,   14,   15,   15],
         [2408, 2387, 2365, ...,   14,   15,   14],
         [2428, 2411, 2397, ...,   13,   13,   13]], dtype=int16)

• Coordinates: (2)
  • lat
    (y, x)
    float64
    47.0 47.0 47.0 ... 45.0 45.0 45.0

    array([[47.        , 47.        , 47.        , ..., 47.        ,
            47.        , 47.        ],
           [46.99916667, 46.99916667, 46.99916667, ..., 46.99916667,
            46.99916667, 46.99916667],
           [46.99833333, 46.99833333, 46.99833333, ..., 46.99833333,
            46.99833333, 46.99833333],
           ...,
           [45.00166667, 45.00166667, 45.00166667, ..., 45.00166667,
            45.00166667, 45.00166667],
           [45.00083333, 45.00083333, 45.00083333, ..., 45.00083333,
            45.00083333, 45.00083333],
           [45.        , 45.        , 45.        , ..., 45.        ,
            45.        , 45.        ]])

  • lon
    (y, x)
    float64
    7.0 7.001 7.002 ... 11.0 11.0 11.0

    array([[ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ],
           [ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ],
           [ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ],
           ...,
           [ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ],
           [ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ],
           [ 7.        ,  7.00083333,  7.00166667, ..., 10.99833333,
            10.99916667, 11.        ]])

• Attributes: (0)

Plot Clutter on DEM

fig = plt.figure(figsize=(13, 10))
ax1 = fig.add_subplot(111)

swpx.CMAP.where(swpx.CMAP == 1).plot(
    x="x", y="y", ax=ax1, vmin=0, vmax=1, cmap="turbo", add_colorbar=False
)
ax1.set_title("Reflectivity corr")

ax1.plot(swpx.longitude.values, swpx.latitude.values, marker="*", c="r")

elev.plot(x="lon", y="lat", ax=ax1, zorder=-2, cmap="terrain")

<matplotlib.collections.QuadMesh at 0x7f3f4c8f44c0>

Use hvplot for interactive zooming and panning

Often it is desirable to quickly zoom and pan in the plots. Although matplotlib has that ability, it still is quite slow. Here hvplot, a holoviews based plotting framework, can be utilized. As frontend bokeh is used.

import hvplot
import hvplot.xarray

We need to rechunk the coordinates as hvplot needs chunked variables and coords.

todo # vergleichbar machen mit beam blockage angle/entfernung

cl = (
    swpx.CMAP.where(swpx.CMAP == 1)
    .chunk(chunks={})
    .hvplot.quadmesh(
        x="x", y="y", cmap="Reds", width=800, height=700, clim=(0, 1), alpha=0.6
    )
)
dm = elev.hvplot.quadmesh(
    x="lon",
    y="lat",
    cmap="terrain",
    width=800,
    height=700,
    clim=(0, 4000),
    rasterize=True,
)
dm * cl

Convert DEM to spherical coordinates

sitecoords = (swpx.longitude.values, swpx.latitude.values, swpx.altitude.values)
r = swpx.range.values
az = swpx.azimuth.values
bw = 0.8
beamradius = wrl.util.half_power_radius(r, bw)

rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
    dem, nodata=-32768.0
)

rlimits = (extent["left"], extent["bottom"], extent["right"], extent["top"])
# Clip the region inside our bounding box
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[ind[1] : ind[3], ind[0] : ind[2], ...]
rastervalues = rastervalues[ind[1] : ind[3], ind[0] : ind[2]]

polcoords = np.dstack([swpx.x.values, swpx.y.values])
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart_to_irregular_spline(
    rastercoords, rastervalues, polcoords, order=3, prefilter=False
)

Partial and Cumulative Beamblockage

PBB = wrl.qual.beam_block_frac(polarvalues, swpx.z.values, beamradius)
PBB = np.ma.masked_invalid(PBB)
print(PBB.shape)

(360, 200)

/srv/conda/envs/notebook/lib/python3.9/site-packages/wradlib/qual.py:127: RuntimeWarning: invalid value encountered in sqrt
  numer = (ya * np.sqrt(a**2 - y**2)) + (a * np.arcsin(ya)) + (np.pi * a / 2.0)
/srv/conda/envs/notebook/lib/python3.9/site-packages/wradlib/qual.py:127: RuntimeWarning: invalid value encountered in arcsin
  numer = (ya * np.sqrt(a**2 - y**2)) + (a * np.arcsin(ya)) + (np.pi * a / 2.0)

CBB = wrl.qual.cum_beam_block_frac(PBB)
print(CBB.shape)

(360, 200)

Plotting Beamblockage

# just a little helper function to style x and y axes of our maps
def annotate_map(ax, cm=None, title=""):
    xticks = ax.get_xticks()
    ticks = (xticks / 1000).astype(int)
    ax.set_xticks(xticks)
    ax.set_xticklabels(ticks)
    yticks = ax.get_yticks()
    ticks = (yticks / 1000).astype(int)
    ax.set_yticks(yticks)
    ax.set_yticklabels(ticks)
    ax.set_xlabel("Kilometers")
    ax.set_ylabel("Kilometers")
    if not cm is None:
        plt.colorbar(cm, ax=ax)
    if not title == "":
        ax.set_title(title)
    ax.grid()

alt = swpx.z.values
fig = plt.figure(figsize=(15, 12))

# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)

# azimuth angle
angle = 270

# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(
    polarvalues, ax=ax1, r=r, az=az, cmap=mpl.cm.terrain, vmin=0.0
)
ax1.plot(
    [0, np.sin(np.radians(angle)) * 1e5], [0, np.cos(np.radians(angle)) * 1e5], "r-"
)
ax1.plot(sitecoords[0], sitecoords[1], "ro")
annotate_map(ax1, dem, "Terrain within {0} km range".format(np.max(r / 1000.0) + 0.1))
ax1.set_xlim(-100000, 100000)
ax1.set_ylim(-100000, 100000)

# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(CBB, ax=ax2, r=r, az=az, cmap=mpl.cm.PuRd, vmin=0, vmax=1)
annotate_map(ax2, cbb, "Beam-Blockage Fraction")
ax2.set_xlim(-100000, 100000)
ax2.set_ylim(-100000, 100000)

# Plot single ray terrain profile on ax3
(bc,) = ax3.plot(r / 1000.0, alt[angle, :], "-b", linewidth=3, label="Beam Center")
(b3db,) = ax3.plot(
    r / 1000.0,
    (alt[angle, :] + beamradius),
    ":b",
    linewidth=1.5,
    label="3 dB Beam width",
)
ax3.plot(r / 1000.0, (alt[angle, :] - beamradius), ":b")
ax3.fill_between(r / 1000.0, 0.0, polarvalues[angle, :], color="0.75")
ax3.set_xlim(0.0, np.max(r / 1000.0) + 0.1)
ax3.set_ylim(0.0, 3000)
ax3.set_xlabel("Range (km)")
ax3.set_ylabel("Altitude (m)")
ax3.grid()

axb = ax3.twinx()
(bbf,) = axb.plot(r / 1000.0, CBB[angle, :], "-g", label="BBF")
axb.spines["right"].set_color("g")
axb.tick_params(axis="y", colors="g")
axb.set_ylabel("Beam-blockage fraction", c="g")
axb.set_ylim(0.0, 1.0)
axb.set_xlim(0.0, np.max(r / 1000.0) + 0.1)


legend = ax3.legend(
    (bc, b3db, bbf),
    ("Beam Center", "3 dB Beam width", "BBF"),
    loc="upper left",
    fontsize=10,
)

Plotting Beamblockage on Curvelinear Grid

Here you get an better impression of the actual beam progression .

def height_formatter(x, pos):
    x = (x - 6370000) / 1000
    fmt_str = "{:g}".format(x)
    return fmt_str


def range_formatter(x, pos):
    x = x / 1000.0
    fmt_str = "{:g}".format(x)
    return fmt_str

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

cgax, caax, paax = wrl.vis.create_cg(fig=fig, rot=0, scale=1)

# azimuth angle
angle = 270

# fix grid_helper
er = 6370000
gh = cgax.get_grid_helper()
gh.grid_finder.grid_locator2._nbins = 80
gh.grid_finder.grid_locator2._steps = [1, 2, 4, 5, 10]

# calculate beam_height and arc_distance for ke=1
# means line of sight
bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.0)
ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.0)
nn0 = np.zeros_like(r)
# for nice plotting we assume earth_radius = 6370000 m
ecp = nn0 + er
# theta (arc_distance sector angle)
thetap = -np.degrees(ade / er) + 90.0

# zero degree elevation with standard refraction
bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er)

# plot (ecp is earth surface normal null)
(bes,) = paax.plot(thetap, ecp, "-k", linewidth=3, label="Earth Surface NN")
(bc,) = paax.plot(thetap, ecp + alt[angle, :], "-b", linewidth=3, label="Beam Center")
(bc0r,) = paax.plot(thetap, ecp + bh0, "-g", label="0 deg Refraction")
(bc0n,) = paax.plot(thetap, ecp + bhe, "-r", label="0 deg line of sight")
(b3db,) = paax.plot(
    thetap, ecp + alt[angle, :] + beamradius, ":b", label="+3 dB Beam width"
)
paax.plot(thetap, ecp + alt[angle, :] - beamradius, ":b", label="-3 dB Beam width")

# orography
paax.fill_between(thetap, ecp, ecp + polarvalues[angle, :], color="0.75")

# shape axes
cgax.set_xlim(0, np.max(ade))
cgax.set_ylim([ecp.min() - 1000, ecp.max() + 2500])
caax.grid(True, axis="x")
cgax.grid(True, axis="y")
cgax.axis["top"].toggle(all=False)
caax.yaxis.set_major_locator(
    mpl.ticker.MaxNLocator(steps=[1, 2, 4, 5, 10], nbins=20, prune="both")
)
caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())
caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))
caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))

caax.set_xlabel("Range (km)")
caax.set_ylabel("Altitude (km)")

legend = paax.legend(
    (bes, bc0n, bc0r, bc, b3db),
    (
        "Earth Surface NN",
        "0 deg line of sight",
        "0 deg std refraction",
        "Beam Center",
        "3 dB Beam width",
    ),
    loc="lower left",
    fontsize=10,
)

Use Clutter and Beamblockage as Quality Index

Simple masking with cumulative beam blockage and Gabella.

swpx = swpx.assign({"CBB": (["azimuth", "range"], CBB)})
# recalculate georeferencing for AEQD
swpx = swpx.pipe(wrl.georef.georeference_dataset)

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

ax1 = fig.add_subplot(121)
swpx.reflectivity.plot(x="x", y="y", ax=ax1, cmap="turbo", vmin=0, vmax=60)
ax1.set_title(f"Signal Processor - {swpx.time.values.astype('M8[s]')}")
ax1.set_aspect("equal")

ax2 = fig.add_subplot(122)
# CBB > 0.5, CMAP == 1, RHOHV < 0.8 is masked
swpx.where(
    (swpx.CBB <= 0.5)
    & (swpx.CMAP < 1.0)
    & (swpx.uncorrected_cross_correlation_ratio >= 0.8)
).reflectivity_hh_clut.plot(x="x", y="y", ax=ax2, cmap="turbo", vmin=0, vmax=60)
ax2.set_title(f"Gabella+CBB+RHOHV - {swpx.time.values.astype('M8[s]')}")
ax2.set_aspect("equal")

---

Summary

We’ve just learned how to use \(\omega radlib\)’s Gabella clutter detection for single sweeps. Wrapping numpy based functions for use with xarray.apply_ufunc has been shown. We’ve looked into digital elevation maps and beam blockage calculations.

What’s next?

In the next notebook we dive into processing of differential phase.

Resources and references

• SRTM - NASA - EARTHDATA

• xarray

• apply_ufunc

• dask

• gdal

• hvplot

• wradlib xarray backends

• CfRadial1

• OPERA ODIM_H5

• WMO JET-OWR

previous

wradlib radar data io, visualisation, gridding and gis export

next

wradlib Phase Processing - System PhiDP - ZPHI-Method