Errors in cell to-track merger? Pulls in very distant features on 17 June.

[deeplycloudy]

I was looking at the Track dataset for 17 June and noticed that there are very distant features that are indicated as part of the same track. I think this a bug in the cell to track merger process, because the individual cells make sense.

I have definitely been tripped up by indexing in xarray before, which will sometimes use the labeled index variables when you think you're indexing by integer position in an array. It might be that, or something else…

Here's some code that looks at Track 726, which has the cell I looked at for our ARM ASR poster, @kelcyno.

import xarray as xr
from pyxlma.lmalib.traversal import OneToManyTraversal

trackpath = '/data/Houston/TRACER_TRACKING_graupel/tobac_Save_20220617/'
track = xr.open_dataset(trackpath+"Track_features_merges.nc")

# Find tracks that might be for a cell of interest
candidates = ((track.feature_projection_x_coordinate > -50e3) & (track.feature_projection_x_coordinate < 0e3)
             & (track.feature_projection_y_coordinate > 0e3) & (track.feature_projection_y_coordinate < 100e3)
             & (track.feature_time > np.datetime64('2022-06-17T20:00')) & (track.feature_time < np.datetime64('2022-06-17T21:30'))
             )
candidate_features = track[{'feature':candidates}]
candidate_cell_ids = set(np.unique(candidate_features.feature_parent_cell_id))
candidate_cell_ids.remove(-1)

traversal = OneToManyTraversal(track, ['track', 'cell', 'feature_id'], ['cell_parent_track_id', 'feature_parent_cell_id'])
candidate_ds = traversal.reduce_to_entities('cell', list(candidate_cell_ids))

feature_min_time_idx = candidate_ds.feature_time_index.min()
feature_max_time_idx = candidate_ds.feature_time_index.max()

from matplotlib.ticker import FuncFormatter

track_fig, track_ax = plt.subplots(1,1)

def max_feature_range(ds):
    r_sq = one_track.feature_projection_x_coordinate**2 + one_track.feature_projection_y_coordinate**2
    return np.sqrt(r_sq).max()

def time_for_index(time_index, pos):
    time = radar.time.data[int(time_index)]
    return str(time).split('T')[-1].split('.')[0]

track_of_interest = None

for track_id in candidate_ds.track:
    if track_id != -1:
        one_track = traversal.reduce_to_entities('track', [track_id])
        print(int(track_id), one_track.dims, "\n    ",
              list(one_track.cell.data), max_feature_range(one_track).data/1000.0
             )
        if track_id == 726:
            track_of_interest = one_track
            label = "T"+str(int(track_id)) + "-C" + str(list(np.unique(one_track.cell)))
            scart = track_ax.scatter(one_track.feature_projection_x_coordinate, one_track.feature_projection_y_coordinate,
                          c=one_track.feature_time_index, linewidth=1, vmin=feature_min_time_idx, vmax=feature_max_time_idx)
            track_ax.text(one_track.feature_projection_x_coordinate[0], one_track.feature_projection_y_coordinate[0],
                          label, fontdict={'fontsize':8})
track_ax.set_aspect(1.0)
cax = track_fig.colorbar(scart)
cax.ax.yaxis.set_major_formatter(FuncFormatter(time_for_index))
track_fig.tight_layout()

This is a plot of every feature in track 726:

# Cross-check we found all cells for this track.
print(track[{'cells':(track.cell_parent_track_id == 726)}].cell)
print(track_of_interest)

# Print all the feature IDs and locations for this track
for v in (list(
      zip(track_of_interest.feature_projection_x_coordinate.data/1000.0,
          track_of_interest.feature_projection_y_coordinate.data/1000.0,
          track_of_interest.feature_parent_cell_id.data,
          )
     )): print(v)

The last loop shows every feature in track 726, some of which are hundreds of km away from the "real" track near (-25, 30).

(-22.085710573660748, 26.20737207314187, 7138)
(-22.148923531226643, 26.235096117006776, 7138)
(-22.72264029709649, 26.500807371266433, 7138)
(-23.662392191245374, 26.736950800015165, 7138)
(-24.27681369284474, 26.844425436679046, 7138)
(-24.456049130850403, 27.676398083976437, 7138)
(-24.749926381686436, 27.384684187123128, 7138)
(-25.03196335685965, 28.02687769403201, 7138)
(-24.970209604911958, 27.9796324161494, 7138)
(-26.805567196329974, 29.0771293228201, 7138)
(-28.13334763836238, 29.673138059791597, 7138)
(-29.1511257962961, 29.908180540102478, 7138)
(-29.817523340240257, 30.61094324440154, 7138)
(-29.92329613221233, 31.112879617787428, 7138)
(-29.281225747178183, 31.026177912864625, 7138)
(-29.821445259893267, 31.430987361698214, 7138)
(-28.816767123907, 31.55097489334014, 7138)
(-28.5504252414855, 31.695691209669747, 7138)
(-28.925395552670068, 32.196950014023344, 7138)
(39.364427369996406, 35.29095252548075, 7139)
(37.53013882843942, 36.11846364341226, 7139)
(37.36943430790154, 36.85440709367066, 7139)
(36.78761819413245, 37.21294454456631, 7139)
(34.32245948340682, 38.5, 7139)
(35.5, 38.85872561777984, 7139)
(34.37427493719031, 39.48932195287347, 7139)
(35.29960859072696, 39.70277563035893, 7139)
(35.11398553662764, 40.09786959763488, 7139)
(158.75456886619173, 118.37988413756904, 7145)
(158.2381499809266, 119.4122559228818, 7145)
(156.56001306792047, 121.30793664215764, 7145)
(155.740109769247, 122.3846404054047, 7145)
(154.73023707477302, 123.09741404513147, 7145)
(-56.92851108672937, 158.57312020659447, 7148)
(-57.923898371618975, 158.70569027485425, 7148)
(-67.22173289535002, 158.68971314710956, 7148)
(-67.25763641261125, 158.05151469201456, 7148)
(-67.43292471987695, 154.68015739291167, 7148)
(-68.1868336776281, 153.78470971109272, 7148)
(-67.86720580700978, 153.91876177192182, 7148)
(-69.72153344668166, 153.33221852216292, 7148)
(-70.99393955964882, 154.2116604761444, 7148)
(-3.337778233902299, 175.56337027642468, 7151)
(-4.045376363125314, 175.72783805334063, 7151)
(-4.414234146972802, 176.62453704950968, 7151)
(-4.731451949891039, 176.89471634891987, 7151)
(-5.3684767693121955, 177.26531505614543, 7151)
(170.81677938645782, 206.2245011814723, 7152)
(171.0, 205.5, 7152)
(170.49203585211427, 204.44465162215397, 7152)
(170.0, 205.5, 7152)
(168.0216433502951, 204.6012950754265, 7152)

[kelcyno]

I found the error in the merge_split algorithm. The error is in how the track numbers are reassigned an ID counting up from 0.

`import xarray as xr trackpath = '/Users/kelcy/PYTHON/TRACER-PAWS-NEXRAD-LMA/scripts/Tracking/tobac_Save_20220617/' track = xr.open_dataset(trackpath+"Track_features_merges.nc") track = both_ds from matplotlib.ticker import FuncFormatter track_fig, track_ax = plt.subplots(1,1) t_ids = [] for i in [7138, 7139, 7145, 7148, 7151, 7152]: ids = np.where(track['feature_parent_cell_id'] == i) t_ids.append(track['feature_parent_track_id'][ids].values) t_ids = list(flatten(t_ids)) print(np.unique(t_ids)) for track_id in t_ids:#candidate_ds.track: ids = np.where(track['feature_parent_track_id'] == track_id) cells = np.unique(track['feature_parent_cell_id'][ids]) label = "T"+str(int(track_id)) + "-C" + str(list(np.unique(cells))) scart = track_ax.scatter(track['feature_projection_x_coordinate'][ids], track['feature_projection_y_coordinate'][ids], c=track['feature_time_index'][ids], linewidth=1)#, vmin=feature_min_time_idx, vmax=feature_max_time_idx) track_ax.text(track['feature_projection_x_coordinate'][ids][-1],track['feature_projection_y_coordinate'][ids][-1], f'{int(track_id)}', fontsize = 'small',rotation = 'vertical') track_ax.set_aspect(1.0) cax = track_fig.colorbar(scart) track_fig.tight_layout()'

The above gives me the following graph for each cell in @deeplycloudy 's track 726. The cells are labeled in the figure with their track number (now, all individual tracks and are not merged together).

[Liyrrr]

Hi， good job! But do I have to download radar data with a scan mode of RHI to use this code? Thanks！

[kelcyno]

@Liyrrr The merge_split.py file is a post processing step to the tobac tracking code. The NEXRAD tracking file knb_tobac_tracking.py ingests gridded NEXRAD data. Whatever version of radar data you use, it must be gridded to use tobac and the merge_split function.