Ctrl+K

Hydrometeorclassification

Hydrometeorclassification#

[1]:

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

warnings.filterwarnings("ignore")
try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()
from scipy import interpolate
import datetime as dt
import glob

/home/runner/micromamba/envs/wradlib-tests/lib/python3.11/site-packages/h5py/__init__.py:36: UserWarning: h5py is running against HDF5 1.14.3 when it was built against 1.14.2, this may cause problems
  _warn(("h5py is running against HDF5 {0} when it was built against {1}, "

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.

Precipitation Types#

[2]:

pr_types = wrl.classify.pr_types
for k, v in pr_types.items():
    print(str(k) + " - ".join(v))

0LR - Light Rain
1MR - Moderate Rain
2HR - Heavy Rain
3LD - Large Drops
4HL - Hail
5RH - Rain/Hail
6GH - Graupel/Hail
7DS - Dry Snow
8WS - Wet Snow
9HC - H Crystals
10VC - V Crystals
11NP - No Precip

Membership Functions#

Load 2D Membership Functions#

[3]:

filename = wrl.util.get_wradlib_data_file("misc/msf_xband_v1.nc")
msf = xr.open_dataset(filename)
display(msf)

Downloading file 'misc/msf_xband_v1.nc' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/misc/msf_xband_v1.nc' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

<xarray.Dataset>
Dimensions:  (idp: 55, hmc: 11, trapezoid: 4)
Coordinates:
  * idp      (idp) float64 -10.0 -8.0 -6.0 -4.0 -2.0 ... 92.0 94.0 96.0 98.0
  * hmc      (hmc) <U2 'LR' 'MR' 'HR' 'LD' 'HL' 'RH' 'GH' 'DS' 'WS' 'HC' 'VC'
Dimensions without coordinates: trapezoid
Data variables:
    ZH       (hmc, idp, trapezoid) float64 ...
    ZDR      (hmc, idp, trapezoid) float64 ...
    RHO      (hmc, idp, trapezoid) float64 ...
    KDP      (hmc, idp, trapezoid) float64 ...
    TEMP     (hmc, idp, trapezoid) float64 ...
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

Plot 2D Membership Functions#

[4]:

minmax = [(-10, 100), (-1, 6), (0.0, 1.0), (-5, 35), (-65, 45)]

for i, pr in enumerate(pr_types.values()):
    if pr[0] == "NP":
        continue
    fig = plt.figure(figsize=(10, 8))
    t = fig.suptitle(" - ".join(pr))
    t.set_y(1.02)
    hmc = msf.sel(hmc=pr[0])
    for k, p in enumerate(hmc.data_vars.values()):
        p = p.where(p != 0)
        ax = fig.add_subplot(3, 2, k + 1)
        p.sel(trapezoid=0).plot(x="idp", c="k", lw=1.0, ax=ax)
        p.sel(trapezoid=1).plot(x="idp", c="k", lw=2.0, ax=ax)
        p.sel(trapezoid=2).plot(x="idp", c="k", lw=2.0, ax=ax)
        p.sel(trapezoid=3).plot(x="idp", c="k", lw=1.0, ax=ax)
        ax.set_xlim((hmc.idp.min(), hmc.idp.max()))
        ax.margins(x=0.05, y=0.05)
        t = ax.set_title(f"{p.long_name}")
        ax.set_ylim(minmax[k])
    fig.tight_layout()
plt.show()
../../_images/notebooks_classify_2d_hmc_10_0.png
../../_images/notebooks_classify_2d_hmc_10_1.png
../../_images/notebooks_classify_2d_hmc_10_2.png
../../_images/notebooks_classify_2d_hmc_10_3.png
../../_images/notebooks_classify_2d_hmc_10_4.png
../../_images/notebooks_classify_2d_hmc_10_5.png
../../_images/notebooks_classify_2d_hmc_10_6.png
../../_images/notebooks_classify_2d_hmc_10_7.png
../../_images/notebooks_classify_2d_hmc_10_8.png
../../_images/notebooks_classify_2d_hmc_10_9.png
../../_images/notebooks_classify_2d_hmc_10_10.png

Use Sounding Data#

Retrieve Sounding Data#

To get the temperature as additional discriminator we use radiosonde data from the University of Wyoming.

The function get_radiosonde tries to find the next next available radiosonde measurement on the given date.

[5]:

rs_time = dt.datetime(2014, 6, 10, 12, 0)
import urllib

try:
    rs_data, rs_meta = wrl.io.get_radiosonde(10410, rs_time)
except (urllib.error.HTTPError, urllib.error.URLError):
    dataf = wrl.util.get_wradlib_data_file("misc/radiosonde_10410_20140610_1200.h5")
    rs_data, _ = wrl.io.from_hdf5(dataf)
    metaf = wrl.util.get_wradlib_data_file("misc/radiosonde_10410_20140610_1200.json")
    with open(metaf, "r") as infile:
        import json

        rs_meta = json.load(infile)
rs_meta

[5]:

{'Station identifier': 'EDZE',
 'Station number': 10410,
 'Observation time': datetime.datetime(2014, 6, 10, 12, 0),
 'Station latitude': 51.4,
 'Station longitude': 6.97,
 'Station elevation': 147.0,
 'Showalter index': 1.65,
 'Lifted index': -5.85,
 'LIFT computed using virtual temperature': -6.24,
 'SWEAT index': 89.02,
 'K index': 23.7,
 'Cross totals index': 17.3,
 'Vertical totals index': 31.3,
 'Totals totals index': 48.6,
 'Convective Available Potential Energy': 1542.9,
 'CAPE using virtual temperature': 1644.44,
 'Convective Inhibition': -139.39,
 'CINS using virtual temperature': -79.9,
 'Equilibrum Level': 202.46,
 'Equilibrum Level using virtual temperature': 202.39,
 'Level of Free Convection': 736.36,
 'LFCT using virtual temperature': 773.09,
 'Bulk Richardson Number': 57.4,
 'Bulk Richardson Number using CAPV': 61.18,
 'Temp [K] of the Lifted Condensation Level': 288.22,
 'Pres [hPa] of the Lifted Condensation Level': 882.28,
 'Equivalent potential temp [K] of the LCL': 334.9,
 'Mean mixed layer potential temperature': 298.74,
 'Mean mixed layer mixing ratio': 12.42,
 '1000 hPa to 500 hPa thickness': 5663.0,
 'Precipitable water [mm] for entire sounding': 28.11,
 'quantity': {'PRES': 'hPa',
  'HGHT': 'm',
  'TEMP': 'C',
  'DWPT': 'C',
  'FRPT': 'C',
  'RELH': '%',
  'RELI': '%',
  'MIXR': 'g/kg',
  'DRCT': 'deg',
  'SKNT': 'knot',
  'THTA': 'K',
  'THTE': 'K',
  'THTV': 'K'}}

Extract Temperature and Height#

[6]:

stemp = rs_data["TEMP"]
sheight = rs_data["HGHT"]
# remove nans
idx = np.isfinite(stemp)
stemp = stemp[idx]
sheight = sheight[idx]

Create DataArray#

[7]:

stemp_da = xr.DataArray(
    data=stemp,
    dims=["height"],
    coords=dict(
        height=(["height"], sheight),
    ),
    attrs=dict(
        description="Temperature.",
        units="degC",
    ),
)
display(stemp_da)

<xarray.DataArray (height: 97)>
array([ 25.6,  19.8,  21.6,  21.6,  21.6,  19.7,  16.7,  16.4,  15.6,
        13.4,   5. ,   3.2,   1.8,  -5.3,  -6.7,  -9.3, -12.5, -14.9,
       -18.5, -18.9, -23.5, -25.5, -31.9, -32.7, -37. , -40.7, -41.9,
       -46.5, -50.5, -51.1, -60.7, -61.5, -65.7, -66.1, -64.1, -60.5,
       -59.1, -59.1, -59.4, -59.9, -60.4, -60.7, -58.5, -57.4, -55.5,
       -54.3, -54.9, -55.1, -55.3, -55.5, -55.9, -56.1, -55.4, -54.4,
       -53.4, -52.5, -52.7, -53.2, -55.6, -55.9, -55.5, -55.1, -54.9,
       -54.8, -54.6, -54.5, -54.3, -54.2, -54.2, -54.1, -53.9, -53.8,
       -52.5, -52.2, -51.8, -51.5, -50.7, -49.5, -48.9, -48.5, -48.6,
       -49.5, -49.7, -49.9, -49.1, -48.5, -47.6, -46.7, -45.8, -44.8,
       -43.9, -39.5, -37.3, -35.9, -36.1, -36.7, -35.5])
Coordinates:
  * height   (height) float64 147.0 744.0 828.0 ... 3.155e+04 3.228e+04
Attributes:
    description:  Temperature.
    units:        degC

Interpolate to higher resolution#

[8]:

hmax = 30000.0
ht = np.arange(0.0, hmax)
itemp_da = stemp_da.interp({"height": ht})
display(itemp_da)

<xarray.DataArray (height: 30000)>
array([         nan,          nan,          nan, ..., -38.47842779,
       -38.47440585, -38.47038391])
Coordinates:
  * height   (height) float64 0.0 1.0 2.0 3.0 4.0 ... 3e+04 3e+04 3e+04 3e+04
Attributes:
    description:  Temperature.
    units:        degC

Fix Temperature below first measurement#

[9]:

itemp_da = itemp_da.bfill(dim="height")

Plot Temperature Profile#

[10]:

fig = plt.figure(figsize=(5, 10))
ax = fig.add_subplot(111)
itemp_da.plot(y="height", ax=ax, marker="o", zorder=0, c="r")
stemp_da.to_dataset(name="stemp").plot.scatter(
    x="stemp", y="height", ax=ax, marker="o", c="b", zorder=1
)
ax.grid(True)
../../_images/notebooks_classify_2d_hmc_23_0.png

Prepare Radar Data#

Load Radar Data#

[11]:

# read the radar volume scan
filename = "hdf5/2014-06-09--185000.rhi.mvol"
filename = wrl.util.get_wradlib_data_file(filename)

Downloading file 'hdf5/2014-06-09--185000.rhi.mvol' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/2014-06-09--185000.rhi.mvol' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

Extract data for georeferencing#

[12]:

import xradar as xd

swp = xr.open_dataset(
    filename, engine=xd.io.backends.GamicBackendEntrypoint, group="sweep_0", chunks={}
)
swp = xd.util.remove_duplicate_rays(swp)
swp = xd.util.reindex_angle(
    swp, start_angle=0, stop_angle=90, angle_res=0.2, direction=1
)
swp

[12]:

<xarray.Dataset>
Dimensions:            (range: 667, elevation: 450)
Coordinates:
  * range              (range) float32 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
  * elevation          (elevation) float64 0.1 0.3 0.5 0.7 ... 89.5 89.7 89.9
    azimuth            (elevation) float64 dask.array<chunksize=(450,), meta=np.ndarray>
    time               (elevation) datetime64[ns] dask.array<chunksize=(450,), meta=np.ndarray>
    longitude          float64 ...
    latitude           float64 ...
    altitude           float64 ...
Data variables: (12/17)
    KDP                (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP              (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZV               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    RHOHV              (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    ZDR                (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    sweep_mode         <U3 ...
    sweep_number       int64 ...
    prt_mode           <U7 ...
    follow_mode        <U7 ...
    sweep_fixed_angle  float64 ...
[13]:

swp.azimuth.load(), swp.time.load()

[13]:

(<xarray.DataArray 'azimuth' (elevation: 450)>
 array([225.        , 225.        , 225.        , 225.00549316,
        405.00549316, 225.00549316, 405.00549316, 225.        ,
        225.00549316, 225.01098633, 405.00549316, 225.00549316,
        405.00549316, 225.00549316, 225.01098633, 225.01098633,
        405.00549316, 225.00549316, 405.00549316, 225.        ,
        225.00549316, 225.01098633, 405.00549316, 225.00549316,
        225.01098633, 405.00549316, 225.        , 225.01098633,
        225.        , 225.        , 225.        , 225.        ,
        225.        , 225.        , 225.00549316, 225.01098633,
        405.00549316, 225.00549316, 225.01098633, 405.00549316,
        225.00549316, 405.00549316, 225.00549316, 225.00549316,
        405.00549316, 225.00549316, 405.00549316, 225.        ,
        225.00549316, 405.00549316, 225.        , 225.        ,
        225.00549316, 405.00549316, 225.00549316, 405.00549316,
        225.        , 225.        , 225.        , 225.        ,
        225.00549316, 225.        , 225.        , 225.        ,
        225.        , 225.00549316, 405.00549316, 225.00549316,
        225.01098633, 405.00549316, 225.        , 225.        ,
        225.        , 225.        , 225.        , 225.        ,
        225.        , 225.        , 225.        , 225.        ,
 ...
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367, 224.98901367, 224.98901367,
        224.98901367, 224.98901367])
 Coordinates:
   * elevation  (elevation) float64 0.1 0.3 0.5 0.7 0.9 ... 89.3 89.5 89.7 89.9
     azimuth    (elevation) float64 225.0 225.0 225.0 225.0 ... 225.0 225.0 225.0
     time       (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 2014-06-09T...
     longitude  float64 7.072
     latitude   float64 50.73
     altitude   float64 99.5
 Attributes:
     standard_name:  ray_azimuth_angle
     long_name:      azimuth_angle_from_true_north
     units:          degrees
     axis:           radial_azimuth_coordinate,
 <xarray.DataArray 'time' (elevation: 450)>
 array(['2014-06-09T18:50:01.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:04.000000000', '2014-06-09T18:50:04.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:05.000000000', '2014-06-09T18:50:05.000000000',
        '2014-06-09T18:50:06.000000000', '2014-06-09T18:50:06.000000000',
 ...
        '2014-06-09T18:50:24.000000000', '2014-06-09T18:50:24.000000000',
        '2014-06-09T18:50:24.000000000', '2014-06-09T18:50:24.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:25.000000000', '2014-06-09T18:50:25.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:26.000000000',
        '2014-06-09T18:50:26.000000000', '2014-06-09T18:50:27.000000000'],
       dtype='datetime64[ns]')
 Coordinates:
   * elevation  (elevation) float64 0.1 0.3 0.5 0.7 0.9 ... 89.3 89.5 89.7 89.9
     azimuth    (elevation) float64 225.0 225.0 225.0 225.0 ... 225.0 225.0 225.0
     time       (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 2014-06-09T...
     longitude  float64 7.072
     latitude   float64 50.73
     altitude   float64 99.5
 Attributes:
     standard_name:  time)

[14]:

swp.azimuth.values

[14]:

array([225.        , 225.        , 225.        , 225.00549316,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 225.01098633, 405.00549316, 225.00549316,
       405.00549316, 225.00549316, 225.01098633, 225.01098633,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 225.01098633, 405.00549316, 225.00549316,
       225.01098633, 405.00549316, 225.        , 225.01098633,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.00549316, 225.01098633,
       405.00549316, 225.00549316, 225.01098633, 405.00549316,
       225.00549316, 405.00549316, 225.00549316, 225.00549316,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 405.00549316, 225.        , 225.        ,
       225.00549316, 405.00549316, 225.00549316, 405.00549316,
       225.        , 225.        , 225.        , 225.        ,
       225.00549316, 225.        , 225.        , 225.        ,
       225.        , 225.00549316, 405.00549316, 225.00549316,
       225.01098633, 405.00549316, 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.00549316, 405.00549316, 225.        ,
       225.        , 225.00549316, 225.01098633, 405.00549316,
       225.00549316, 225.01098633, 405.00549316, 225.        ,
       225.        , 225.        , 225.00549316, 405.00549316,
       225.00549316, 405.00549316, 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.00549316, 225.01098633,
       405.00549316, 225.        , 225.00549316, 225.01098633,
       405.00549316, 405.00549316, 225.        , 225.00549316,
       405.00549316, 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       404.99450684, 404.98352051, 224.98352051, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.99450684, 404.99450684, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367])

Get Heights of Radar Bins#

[15]:

swp = swp.wrl.georef.georeference()
swp.azimuth.values

[15]:

array([225.        , 225.        , 225.        , 225.00549316,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 225.01098633, 405.00549316, 225.00549316,
       405.00549316, 225.00549316, 225.01098633, 225.01098633,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 225.01098633, 405.00549316, 225.00549316,
       225.01098633, 405.00549316, 225.        , 225.01098633,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.00549316, 225.01098633,
       405.00549316, 225.00549316, 225.01098633, 405.00549316,
       225.00549316, 405.00549316, 225.00549316, 225.00549316,
       405.00549316, 225.00549316, 405.00549316, 225.        ,
       225.00549316, 405.00549316, 225.        , 225.        ,
       225.00549316, 405.00549316, 225.00549316, 405.00549316,
       225.        , 225.        , 225.        , 225.        ,
       225.00549316, 225.        , 225.        , 225.        ,
       225.        , 225.00549316, 405.00549316, 225.00549316,
       225.01098633, 405.00549316, 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.00549316, 405.00549316, 225.        ,
       225.        , 225.00549316, 225.01098633, 405.00549316,
       225.00549316, 225.01098633, 405.00549316, 225.        ,
       225.        , 225.        , 225.00549316, 405.00549316,
       225.00549316, 405.00549316, 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.00549316, 225.01098633,
       405.00549316, 225.        , 225.00549316, 225.01098633,
       405.00549316, 405.00549316, 225.        , 225.00549316,
       405.00549316, 225.        , 225.        , 225.        ,
       225.        , 225.        , 225.        , 225.        ,
       404.99450684, 404.98352051, 224.98352051, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.99450684, 404.99450684, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367, 224.98901367, 224.98901367,
       224.98901367, 224.98901367])

Plot RHI of Heights#

[16]:

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111)
cmap = mpl.cm.viridis
swp.z.plot(x="gr", y="z", ax=ax, cbar_kwargs=dict(label="Height [m]"))
ax.set_xlabel("Range [m]")
ax.set_ylabel("Height [m]")
ax.grid(True)
plt.show()
../../_images/notebooks_classify_2d_hmc_34_0.png

Get Index into High Res Height Array#

[17]:

def merge_radar_profile(rds, cds):
    cds = cds.interp({"height": rds.z}, method="linear")
    rds = rds.assign({"TEMP": cds})
    return rds

[18]:

hmc_ds = swp.pipe(merge_radar_profile, itemp_da)
display(hmc_ds)

<xarray.Dataset>
Dimensions:            (range: 667, elevation: 450)
Coordinates: (12/15)
  * range              (range) float32 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
  * elevation          (elevation) float64 0.1 0.3 0.5 0.7 ... 89.5 89.7 89.9
    azimuth            (elevation) float64 225.0 225.0 225.0 ... 225.0 225.0
    time               (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 201...
    longitude          float64 7.072
    latitude           float64 50.73
    ...                 ...
    z                  (elevation, range) float64 99.63 99.76 ... 5.009e+04
    gr                 (elevation, range) float64 18.39 93.39 ... 67.49 67.62
    rays               (elevation, range) float64 0.1 0.1 0.1 ... 89.9 89.9 89.9
    bins               (elevation, range) float32 37.5 112.5 ... 4.999e+04
    crs_wkt            int64 0
    height             (elevation, range) float64 99.63 99.76 ... 5.009e+04
Data variables: (12/18)
    KDP                (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    PHIDP              (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZH               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBZV               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    RHOHV              (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    DBTH               (elevation, range) float32 dask.array<chunksize=(450, 667), meta=np.ndarray>
    ...                 ...
    sweep_mode         <U3 ...
    sweep_number       int64 ...
    prt_mode           <U7 ...
    follow_mode        <U7 ...
    sweep_fixed_angle  float64 ...
    TEMP               (elevation, range) float64 25.6 25.6 25.6 ... nan nan nan
[19]:

fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
hmc_ds.TEMP.plot(
    x="gr",
    y="z",
    cmap=cmap,
    ax=ax,
    add_colorbar=True,
    cbar_kwargs=dict(label="Temperature [°C]"),
)
ax.set_xlabel("Range [m]")
ax.set_ylabel("Range [m]")
ax.set_aspect("equal")
ax.set_ylim(0, 30000)
plt.show()
../../_images/notebooks_classify_2d_hmc_38_0.png

HMC Workflow#

Setup Independent Observable \(Z_H\)#

Retrieve membership function values based on independent observable

[20]:

%%time
msf_val = msf.wrl.classify.msf_index_indep(swp.DBZH)
display(msf_val)

<xarray.DataArray 'msf_index_indep' (hmc: 11, obs: 5, elevation: 450,
                                     range: 667, trapezoid: 4)>
dask.array<transpose, shape=(11, 5, 450, 667, 4), dtype=int32, chunksize=(11, 5, 450, 667, 4), chunktype=numpy.ndarray>
Coordinates: (12/16)
  * hmc        (hmc) <U2 'LR' 'MR' 'HR' 'LD' 'HL' 'RH' 'GH' 'DS' 'WS' 'HC' 'VC'
  * obs        (obs) object 'ZH' 'ZDR' 'RHO' 'KDP' 'TEMP'
  * range      (range) float32 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
  * elevation  (elevation) float64 0.1 0.3 0.5 0.7 0.9 ... 89.3 89.5 89.7 89.9
    azimuth    (elevation) float64 225.0 225.0 225.0 225.0 ... 225.0 225.0 225.0
    time       (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 2014-06-09T...
    ...         ...
    y          (elevation, range) float64 -26.52 -79.55 -132.6 ... -61.25 -61.34
    z          (elevation, range) float64 99.63 99.76 ... 5.001e+04 5.009e+04
    gr         (elevation, range) float64 18.39 93.39 168.4 ... 67.49 67.62
    rays       (elevation, range) float64 0.1 0.1 0.1 0.1 ... 89.9 89.9 89.9
    bins       (elevation, range) float32 37.5 112.5 ... 4.991e+04 4.999e+04
    crs_wkt    int64 0
Dimensions without coordinates: trapezoid
CPU times: user 56.4 ms, sys: 105 µs, total: 56.5 ms
Wall time: 56.5 ms

Fuzzyfication#

[21]:

%%time
fu = msf_val.wrl.classify.fuzzyfi(
    hmc_ds, dict(ZH="DBZH", ZDR="ZDR", RHO="RHOHV", KDP="KDP", TEMP="TEMP")
)

CPU times: user 9.46 ms, sys: 0 ns, total: 9.46 ms
Wall time: 9.47 ms

Probability#

[22]:

# weights dataset
w = xr.Dataset(dict(ZH=2.0, ZDR=1.0, RHO=1.0, KDP=1.0, TEMP=1.0))
display(w)

<xarray.Dataset>
Dimensions:  ()
Data variables:
    ZH       float64 2.0
    ZDR      float64 1.0
    RHO      float64 1.0
    KDP      float64 1.0
    TEMP     float64 1.0
[23]:

%%time
prob = fu.wrl.classify.probability(w).compute()
display(prob)

<xarray.DataArray (hmc: 11, elevation: 450, range: 667)>
array([[[0.16666667, 0.5       , 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.5       , 0.66666667, ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.5       , 0.5       , ..., 0.        ,
         0.        , 0.        ],
        ...,
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ]],

       [[0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
...
        [0.29812113, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.31296539, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.32917684, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ]],

       [[0.29480082, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.25280711, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.25808055, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        ...,
        [0.29812113, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.31296539, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.32917684, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ]]])
Coordinates: (12/16)
  * hmc        (hmc) <U2 'LR' 'MR' 'HR' 'LD' 'HL' 'RH' 'GH' 'DS' 'WS' 'HC' 'VC'
  * range      (range) float32 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
  * elevation  (elevation) float64 0.1 0.3 0.5 0.7 0.9 ... 89.3 89.5 89.7 89.9
    azimuth    (elevation) float64 225.0 225.0 225.0 225.0 ... 225.0 225.0 225.0
    time       (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 2014-06-09T...
    longitude  float64 7.072
    ...         ...
    z          (elevation, range) float64 99.63 99.76 ... 5.001e+04 5.009e+04
    gr         (elevation, range) float64 18.39 93.39 168.4 ... 67.49 67.62
    rays       (elevation, range) float64 0.1 0.1 0.1 0.1 ... 89.9 89.9 89.9
    bins       (elevation, range) float32 37.5 112.5 ... 4.991e+04 4.999e+04
    crs_wkt    int64 0
    height     (elevation, range) float64 99.63 99.76 ... 5.001e+04 5.009e+04
CPU times: user 663 ms, sys: 214 ms, total: 877 ms
Wall time: 797 ms

[24]:

# prob = prob.compute()

Classification#

[25]:

cl_res = prob.wrl.classify.classify(threshold=0.0)
display(cl_res)

<xarray.DataArray (hmc: 12, elevation: 450, range: 667)>
array([[[0.16666667, 0.5       , 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.5       , 0.66666667, ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.5       , 0.5       , ..., 0.        ,
         0.        , 0.        ],
        ...,
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ],
        [0.16666667, 0.66666667, 0.5       , ..., 0.        ,
         0.        , 0.        ]],

       [[0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
...
        [0.29812113, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.31296539, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ],
        [0.32917684, 0.33333333, 0.33333333, ..., 0.        ,
         0.        , 0.        ]],

       [[0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ],
        [0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ],
        [0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ],
        ...,
        [0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ],
        [0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ],
        [0.        , 0.        , 0.        , ..., 1.        ,
         1.        , 1.        ]]])
Coordinates: (12/16)
  * hmc        (hmc) <U2 'LR' 'MR' 'HR' 'LD' 'HL' ... 'DS' 'WS' 'HC' 'VC' 'NP'
  * range      (range) float32 37.5 112.5 187.5 ... 4.991e+04 4.999e+04
  * elevation  (elevation) float64 0.1 0.3 0.5 0.7 0.9 ... 89.3 89.5 89.7 89.9
    azimuth    (elevation) float64 225.0 225.0 225.0 225.0 ... 225.0 225.0 225.0
    time       (elevation) datetime64[ns] 2014-06-09T18:50:01 ... 2014-06-09T...
    longitude  float64 7.072
    ...         ...
    z          (elevation, range) float64 99.63 99.76 ... 5.001e+04 5.009e+04
    gr         (elevation, range) float64 18.39 93.39 168.4 ... 67.49 67.62
    rays       (elevation, range) float64 0.1 0.1 0.1 0.1 ... 89.9 89.9 89.9
    bins       (elevation, range) float32 37.5 112.5 ... 4.991e+04 4.999e+04
    crs_wkt    int64 0
    height     (elevation, range) float64 99.63 99.76 ... 5.001e+04 5.009e+04

Compute#

[26]:

%%time
cl_res = cl_res.compute()
cl_res = cl_res.assign_coords(sweep_mode="rhi")

CPU times: user 2.65 ms, sys: 0 ns, total: 2.65 ms
Wall time: 2.65 ms

HMC Results#

Plot Probability of HMC Types#

[27]:

prob = prob.assign_coords(hmc=np.array(list(pr_types.values())).T[1][:11])
prob = prob.where(prob > 0)
prob.plot(x="x", y="y", col="hmc", col_wrap=4, cbar_kwargs=dict(label="Probability"))

[27]:

<xarray.plot.facetgrid.FacetGrid at 0x7fbc9cb82b90>
../../_images/notebooks_classify_2d_hmc_54_1.png

Plot maximum probability#

[28]:

fig = plt.figure(figsize=(10, 6))
cmap = "cubehelix"
im = cl_res.max("hmc").wrl.vis.plot(
    ax=111,
    crs={"angular_spacing": 20.0, "radial_spacing": 12.0, "latmin": 2.5},
    cmap=cmap,
    fig=fig,
)
cgax = plt.gca()
cbar = plt.colorbar(im, ax=cgax, fraction=0.046, pad=0.05)
cbar.set_label("Probability")
cgax.set_xlim(0, 40000)
cgax.set_ylim(0, 14000)
t = cgax.set_title("Hydrometeorclassification", y=1.05)

caax = cgax.parasites[0]
caax.set_xlabel("Range [m]")
caax.set_ylabel("Range [m]")
plt.show()
../../_images/notebooks_classify_2d_hmc_56_0.png

Plot classification result#

[29]:

bounds = np.arange(-0.5, prob.shape[0] + 0.6, 1)
ticks = np.arange(0, prob.shape[0] + 1)
cmap = mpl.cm.get_cmap("cubehelix", len(ticks))
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

[30]:

hydro = cl_res.argmax("hmc")
hydro.attrs = dict(long_name="Hydrometeorclassification")
hydro = hydro.assign_coords(sweep_mode="rhi")

[31]:

fig = plt.figure(figsize=(10, 8))
im = hydro.wrl.vis.plot(
    ax=111,
    crs={"angular_spacing": 20.0, "radial_spacing": 12.0, "latmin": 2.5},
    norm=norm,
    cmap=cmap,
    fig=fig,
)
cgax = plt.gca()
caax = cgax.parasites[0]
paax = cgax.parasites[1]

cbar = plt.colorbar(im, ticks=ticks, ax=cgax, fraction=0.046, norm=norm, pad=0.05)
cbar.set_label("Hydrometeorclass")
caax.set_xlabel("Range [km]")
caax.set_ylabel("Range [km]")
labels = [pr_types[i][1] for i, _ in enumerate(pr_types)]
labels = cbar.ax.set_yticklabels(labels)
t = cgax.set_title("Hydrometeorclassification", y=1.05)
cgax.set_xlim(0, 40000)
cgax.set_ylim(0, 14000)
plt.tight_layout()
../../_images/notebooks_classify_2d_hmc_60_0.png
Contents