Information gathering from plotted field

[MohitDahliya]

Hi, I am new to Tobac, and have used it for the first time. First I implemented only this much of code

from datetime import datetime, timedelta
from pathlib import Path

import numpy as np
import pandas as pd
import xarray as xr
import s3fs
import pyart
from pyproj import Proj, Geod

import matplotlib.pyplot as plt
from matplotlib.pyplot import cm as cmaps
import cartopy.crs as ccrs

%matplotlib inline

Disable a few warnings:

import warnings

warnings.filterwarnings("ignore", category=UserWarning, append=True)
warnings.filterwarnings("ignore", category=RuntimeWarning, append=True)
warnings.filterwarnings("ignore", category=FutureWarning, append=True)
warnings.filterwarnings("ignore", category=pd.io.pytables.PerformanceWarning)

read in radar data

radar_full = pyart.io.read("DLI241228000230-IMD-B_corrected.nc")
radar = radar_full.extract_sweeps([0])

First, we need to grid the input radar data.

radar_grid_pyart = pyart.map.grid_from_radars(
radar,
grid_shape=(41, 401, 401),
grid_limits=((0.0, 1.3e4), (-2e5, 2e5), (-2e5, 2e5),),
fields=['Z']
)

In Py-ART's graphing suite, there is a display class similar to RadarMapDisplay,

but for grids. To plot the grid:

fig = plt.figure(figsize=[12, 8])
pyart.graph.GridMapDisplay(radar_grid_pyart).plot_grid("Z", level=3, vmin=0, vmax=60, fig=fig)
plt.show()
and got the following plot

Afterwards, I applied feature detection and segmentatin to my radar data in the folowing manner :
Z = radar_grid_pyart.to_xarray().Z.drop_vars(['origin_altitude', 'origin_longitude', 'origin_latitude'], errors='ignore')
savedir = Path('C:\Users\mohit.DESKTOP-KCO663F\Documents\Paper2\Radar_data\28Dec2024')
import tobac
import tobac.utils
print("using tobac version", str(tobac.version))

Feature detection

We use multiple thresholds to detect features in the radar data.

feature_detection_params = dict()
feature_detection_params["threshold"] = [30, 40, 50]
feature_detection_params["target"] = "maximum"
feature_detection_params["position_threshold"] = "weighted_diff"
feature_detection_params["n_erosion_threshold"] = 2
feature_detection_params["sigma_threshold"] = 1
feature_detection_params["n_min_threshold"] = 4

Perform feature detection:

print('starting feature detection')
radar_features = tobac.feature_detection.feature_detection_multithreshold(
Z, 0, **feature_detection_params
)
radar_features.to_hdf(savedir / 'Features.h5', 'table')
print('feature detection performed and saved')

radar_features.head()
Z[0, 4,].where(Z[0, 4,] > 0).loc[-2e4:8e4, 0:1e5].plot(vmin=0, vmax=65, cmap="Spectral_r", figsize=(12,8))

plt.scatter(
radar_features["x"],
radar_features["y"],
70,
color="k",
)
plt.title("Reflectivity, 2 km AGL, 2021-05-26 15:56", size=15)

segmentation

parameters_segmentation = dict()
parameters_segmentation["method"] = "watershed"
parameters_segmentation["threshold"] = 240
parameters_segmentation["target"] = "minimum"
parameters_segmentation["seed_3D_flag"] = "box"
parameters_segmentation["seed_3D_size"] = 5

Perform segmentation and save results:

print('Starting segmentation.')
segmentation_mask, segmentation_features = tobac.segmentation.segmentation(
radar_features, Z, dxy=2000, **parameters_segmentation
)
print('segmentation performed, start saving results to files')
segmentation_mask.to_netcdf(savedir / 'Mask_Segmentation_sat.nc', encoding={"segmentation_mask":{"zlib":True, "complevel":4}})
segmentation_features.to_hdf(savedir / 'Features_sat.h5', 'table')
print('segmentation performed and saved')

segmentation_mask.coords

First, let's get the proper z-level (2 km = 2000 m)

z_level = 2000
Z_slice = Z.sel(z=z_level, method='nearest').squeeze()

Get the segmentation mask at the same level

seg_slice = segmentation_mask.sel(z=z_level, method='nearest').squeeze()

Create figure

fig = plt.figure(figsize=[10, 8])
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

Get lon/lat from the grid

lons = radar_grid_pyart.point_longitude['data'][0, :, :]
lats = radar_grid_pyart.point_latitude['data'][0, :, :]

Set extent based on actual data

lon_min, lon_max = lons.min(), lons.max()
lat_min, lat_max = lats.min(), lats.max()
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

Plot reflectivity

contoured = ax.contourf(
lons,
lats,
Z_slice,
transform=ccrs.PlateCarree(),
cmap='jet', vmin=0, vmax=60, levels=31
)

Get unique segments

unique_seg = np.unique(seg_slice)
unique_seg = unique_seg[unique_seg > 0] # Remove 0 and -1
color_map = cmaps.plasma(np.linspace(0, 1, len(unique_seg)))

Plot features without segmented area

curr_feat = radar_features[radar_features["feature"] == 1]
if len(curr_feat) > 0:
ax.scatter(
curr_feat["lon"],
curr_feat["lat"],
70,
transform=ccrs.PlateCarree(),
color="grey",
)

Plot segmentation contours and features

for seg_num, color in zip(unique_seg, color_map):
curr_seg = (seg_slice == seg_num).astype(int)

ax.contour(
    lons,
    lats,
    curr_seg,
    colors=[color],
    levels=[0.5],
    linewidths=3,
    transform=ccrs.PlateCarree(),
)

curr_feat = radar_features[radar_features["feature"] == seg_num]
if len(curr_feat) > 0:
    ax.scatter(
        curr_feat["lon"],
        curr_feat["lat"],
        70,
        transform=ccrs.PlateCarree(),
        color=color,
        edgecolors='black',
        zorder=10
    )

cb = plt.colorbar(contoured, ax=ax)
cb.set_label("Reflectivity (dBZ)", size=14)
cb.ax.tick_params(labelsize=14)
plt.title(f"Reflectivity at {z_level}m AGL with Segmentation", size=15)
plt.show()
and got the following plot

What I wish to know is that what information can I gather from this final plot I have created? Can anyone tell me what type of information can be gathere after applying radar_features and segmentatino to the radar data and plotting it. Thanks in advance for any help!