Hydrometeorclassification - RHI#

The hydrometeorclassification code is based on the paper by Zrnic et.al 2001 utilizing 2D trapezoidal membership functions based on the paper by Straka et. al 2000 adapted by Evaristo et. al 2013 for X-Band.

RHI Stratiform (2014-10-07 02:45)
RHI Convective (2013-07-01 15:40)

import datetime as dt

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

Utility Functions#

import radarpolarimetry as rp

Precipitation Types#

pr_types = dict(wrl.classify.pr_types.values())
for k, v in enumerate(pr_types.items()):
    print(f"{k:02} {' - '.join(v)}")

LR - Light Rain
MR - Moderate Rain
HR - Heavy Rain
LD - Large Drops
HL - Hail
RH - Rain/Hail
GH - Graupel/Hail
DS - Dry Snow
WS - Wet Snow
HC - H Crystals
VC - V Crystals
NP - No Precip

Convective#

filename = "/radar/rp-data/workflow/tmp_20130701T1530_20130701T1600_workflow.nc"
with rp.rpol_data.open(filename) as f:
    filename = f.name
sweep_c = xr.open_dataset(
    filename,
    chunks=dict(time2=1),
).isel(time2=2)
display(sweep_c)

<xarray.Dataset> Size: 94MB
Dimensions:            (elevation: 450, range: 667)
Coordinates: (12/19)
  * elevation          (elevation) float64 4kB 0.1 0.3 0.5 ... 89.5 89.7 89.9
    azimuth            (elevation) float64 4kB dask.array<chunksize=(450,), meta=np.ndarray>
    time               (elevation) datetime64[ns] 4kB dask.array<chunksize=(450,), meta=np.ndarray>
  * range              (range) float32 3kB 37.5 112.5 ... 4.991e+04 4.999e+04
    x                  (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    y                  (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    latitude           float64 8B ...
    altitude           float64 8B ...
    crs_wkt            int64 8B ...
    stdev              float64 8B ...
    NOISE_LEVEL        float64 8B ...
    alpha              float64 8B ...
Data variables: (12/51)
    KDP                (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP              (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZV               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    RHOHV              (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    AH                 (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    KDP_AH             (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP_CORR5        (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH_CORR5         (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH_CORR5         (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ZDR_CORR5          (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>

Load Sounding Data#

# http://weather.uwyo.edu/upperair/sounding.html
rs_time_s = dt.datetime(2014, 10, 7, 2, 0)
rs_time_c = dt.datetime(2013, 7, 1, 12, 0)
print(rs_time_s)
date = rs_time_s
year = date.strftime("%Y")
month = date.strftime("%m")
day = date.strftime("%d")
hour = date.strftime("%H")
print(year, month, day, hour)
import urllib

try:
    rs_data_s, rs_meta_s = wrl.io.get_radiosonde(10410, rs_time_s, cols=(0, 1, 2))
    rs_data_c, rs_meta_c = wrl.io.get_radiosonde(10410, rs_time_c, cols=(0, 1, 2))
    print("downloaded from UWYO")
except urllib.error.HTTPError:
    raise
print(rs_meta_s, rs_meta_c)

2014-10-07 02:00:00
2014 10 07 02

downloaded from UWYO
{'Station identifier': 'EDZE', 'Station number': 10410, 'Observation time': datetime.datetime(2014, 10, 7, 0, 0), 'Station latitude': 51.4, 'Station longitude': 6.97, 'Station elevation': 147.0, 'Showalter index': 7.79, 'Lifted index': 4.55, 'LIFT computed using virtual temperature': 4.63, 'SWEAT index': 204.28, 'K index': 22.9, 'Cross totals index': 14.7, 'Vertical totals index': 24.7, 'Totals totals index': 39.4, 'Convective Available Potential Energy': 28.92, 'CAPE using virtual temperature': 36.38, 'Convective Inhibition': -70.37, 'CINS using virtual temperature': -11.48, 'Equilibrum Level': 689.06, 'Equilibrum Level using virtual temperature': 688.7, 'Level of Free Convection': 766.58, 'LFCT using virtual temperature': 881.01, 'Bulk Richardson Number': 1.33, 'Bulk Richardson Number using CAPV': 1.68, 'Temp [K] of the Lifted Condensation Level': 282.6, 'Pres [hPa] of the Lifted Condensation Level': 913.15, 'Equivalent potential temp [K] of the LCL': 313.42, 'Mean mixed layer potential temperature': 290.05, 'Mean mixed layer mixing ratio': 8.2, '1000 hPa to 500 hPa thickness': 5547.0, 'Precipitable water [mm] for entire sounding': 24.77, 'quantity': {'PRES': 'hPa', 'HGHT': 'm', 'TEMP': 'C'}} {'Station identifier': 'EDZE', 'Station number': 10410, 'Observation time': datetime.datetime(2013, 7, 1, 12, 0), 'Station latitude': 51.4, 'Station longitude': 6.97, 'Station elevation': 147.0, 'Showalter index': 6.1, 'Lifted index': 0.57, 'LIFT computed using virtual temperature': 0.51, 'SWEAT index': 77.7, 'K index': 11.7, 'Cross totals index': 13.7, 'Vertical totals index': 28.7, 'Totals totals index': 42.4, 'Convective Available Potential Energy': 7.03, 'CAPE using virtual temperature': 18.0, 'Convective Inhibition': 0.0, 'CINS using virtual temperature': 0.0, 'Equilibrum Level': 597.47, 'Equilibrum Level using virtual temperature': 589.23, 'Level of Free Convection': 931.56, 'LFCT using virtual temperature': 934.17, 'Bulk Richardson Number': 0.24, 'Bulk Richardson Number using CAPV': 0.62, 'Temp [K] of the Lifted Condensation Level': 284.17, 'Pres [hPa] of the Lifted Condensation Level': 934.17, 'Equivalent potential temp [K] of the LCL': 315.07, 'Mean mixed layer potential temperature': 289.76, 'Mean mixed layer mixing ratio': 8.92, '1000 hPa to 500 hPa thickness': 5543.0, 'Precipitable water [mm] for entire sounding': 19.02, 'quantity': {'PRES': 'hPa', 'HGHT': 'm', 'TEMP': 'C'}}

Extract Temperature Profile#

itemp_s = rp.get_temperature_profile(rs_data_s)
itemp_c = rp.get_temperature_profile(rs_data_c)

Plot Temperature Profile#

fig = plt.figure(figsize=(2.5, 5))
ax = fig.add_subplot(111)
itemp_s.plot(y="height", ax=ax, zorder=0, label="stratiform")
itemp_c.plot(y="height", ax=ax, zorder=0, label="convective")
plt.legend()
ax.grid(True)

Merge Temperature Profile#

hmc_s = sweep_s.pipe(rp.merge_radar_profile, itemp_s)
hmc_c = sweep_c.pipe(rp.merge_radar_profile, itemp_c)
display(hmc_s)
display(hmc_c)

<xarray.Dataset> Size: 106MB
Dimensions:            (elevation: 450, range: 667)
Coordinates: (12/19)
  * elevation          (elevation) float64 4kB 0.1 0.3 0.5 ... 89.5 89.7 89.9
    azimuth            (elevation) float64 4kB ...
    time               (elevation) datetime64[ns] 4kB ...
  * range              (range) float32 3kB 37.5 112.5 ... 4.991e+04 4.999e+04
    x                  (elevation, range) float64 2MB ...
    y                  (elevation, range) float64 2MB ...
    ...                 ...
    latitude           float64 8B ...
    altitude           float64 8B ...
    crs_wkt            int64 8B ...
    stdev              float64 8B ...
    NOISE_LEVEL        float64 8B ...
    alpha              float64 8B ...
Data variables: (12/52)
    KDP                (elevation, range) float32 1MB ...
    PHIDP              (elevation, range) float32 1MB ...
    DBZH               (elevation, range) float32 1MB ...
    DBZV               (elevation, range) float32 1MB ...
    RHOHV              (elevation, range) float32 1MB ...
    DBTH               (elevation, range) float32 1MB ...
    ...                 ...
    KDP_AH             (elevation, range) float64 2MB ...
    PHIDP_CORR5        (elevation, range) float64 2MB ...
    DBTH_CORR5         (elevation, range) float64 2MB ...
    DBZH_CORR5         (elevation, range) float64 2MB ...
    ZDR_CORR5          (elevation, range) float64 2MB ...
    TEMP               (elevation, range) float64 2MB 15.0 15.0 15.0 ... nan nan

<xarray.Dataset> Size: 98MB
Dimensions:            (elevation: 450, range: 667)
Coordinates: (12/20)
  * elevation          (elevation) float64 4kB 0.1 0.3 0.5 ... 89.5 89.7 89.9
    azimuth            (elevation) float64 4kB dask.array<chunksize=(450,), meta=np.ndarray>
    time               (elevation) datetime64[ns] 4kB dask.array<chunksize=(450,), meta=np.ndarray>
  * range              (range) float32 3kB 37.5 112.5 ... 4.991e+04 4.999e+04
    x                  (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    y                  (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    latitude           float64 8B ...
    altitude           float64 8B ...
    crs_wkt            int64 8B ...
    stdev              float64 8B ...
    NOISE_LEVEL        float64 8B ...
    alpha              float64 8B ...
Data variables: (12/52)
    KDP                (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP              (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZV               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    RHOHV              (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH               (elevation, range) float32 1MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    KDP_AH             (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP_CORR5        (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH_CORR5         (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH_CORR5         (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    ZDR_CORR5          (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>
    TEMP               (elevation, range) float64 2MB dask.array<chunksize=(450, 667), meta=np.ndarray>

HMC Workflow#

switch between stratiform and convective processing here

hmc_ds = hmc_s.copy()
# hmc_ds = hmc_c.copy()
# unify masks
hmc_ds = hmc_ds.where(
    ~np.isnan(hmc_ds.DBZH_CORR1)
    & ~np.isnan(hmc_ds.ZDR_CORR1)
    & ~np.isnan(hmc_ds.RHOHV_NC)
    & ~np.isnan(hmc_ds.KDP_VULP)
    & ~np.isnan(hmc_ds.TEMP)
)

Fuzzyfication#

fu = msf_val.wrl.classify.fuzzyfi(
    hmc_ds,
    dict(ZH="DBZH_CORR1", ZDR="ZDR_CORR1", RHO="RHOHV_NC", KDP="KDP_VULP", TEMP="TEMP"),
)
fu = fu.chunk({"hmc": 1})

display(fu)

<xarray.DataArray 'fuzzyfi' (hmc: 11, obs: 5, elevation: 450, range: 667)> Size: 132MB
dask.array<xarray-<this-array>, shape=(11, 5, 450, 667), dtype=float64, chunksize=(1, 5, 450, 667), chunktype=numpy.ndarray>
Coordinates: (12/21)
  * hmc          (hmc) <U2 88B 'LR' 'MR' 'HR' 'LD' 'HL' ... 'DS' 'WS' 'HC' 'VC'
  * obs          (obs) object 40B 'ZH' 'ZDR' 'RHO' 'KDP' 'TEMP'
  * elevation    (elevation) float64 4kB 0.1 0.3 0.5 0.7 ... 89.3 89.5 89.7 89.9
    azimuth      (elevation) float64 4kB dask.array<chunksize=(450,), meta=np.ndarray>
    time         (elevation) datetime64[ns] 4kB dask.array<chunksize=(450,), meta=np.ndarray>
  * range        (range) float32 3kB 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
    ...           ...
    latitude     float64 8B 50.73
    altitude     float64 8B 99.5
    crs_wkt      int64 8B 0
    stdev        float64 8B 0.003062
    NOISE_LEVEL  float64 8B -19.6
    alpha        float64 8B 0.28
Attributes:
    version:      1
    title:        2D Membershipfunctions Hydrometeorclassification for XBand
    institution:  Institute of Geosciences, Meteorology Section, University o...
    history:      created by Raquel Evaristo, adapted to netCDF by Kai Mühlbauer
    comment:      created with wradlib

Hydrometeorclassification - RHI

Contents

Hydrometeorclassification - RHI#

Utility Functions#

Precipitation Types#

Load 2D Membership Functions#

Stratiform#

Convective#

Load Sounding Data#

Extract Temperature Profile#

Plot Temperature Profile#

Merge Temperature Profile#

HMC Workflow#

Setup Independent Observable \(Z_H\)#

Fuzzyfication#

Score#

Classification#

Plot HMC Results#

Plot Score of HMC Types#

Plot maximum score#

Plot classification result#

Footer#