Track blocking events

.travis.yml

@@ -19,7 +19,7 @@ env:
     - BUILD_TYPE=Debug
     - TECA_DIR=/travis_teca_dir
     - TECA_PYTHON_VERSION=3
-    - TECA_DATA_REVISION=168
+    - TECA_DATA_REVISION=169
   jobs:
     - DOCKER_IMAGE=ubuntu IMAGE_VERSION=22.04 IMAGE_NAME=ubuntu_22_04 REQUIRE_NETCDF_MPI=TRUE
     - DOCKER_IMAGE=ubuntu IMAGE_VERSION=22.04 IMAGE_NAME=ubuntu_22_04 REQUIRE_NETCDF_MPI=FALSE

apps/CMakeLists.txt

@@ -47,6 +47,7 @@ teca_add_app(teca_cartesian_mesh_diff LIBS ${teca_app_link}
     FEATURES ${TECA_HAS_NETCDF})
 
 teca_add_python_app(teca_numpy_array_diff)
+teca_add_python_app(teca_pickle_diff)
 
 teca_add_app(teca_cpp_temporal_reduction LIBS ${teca_app_link}
     FEATURES ${TECA_HAS_NETCDF})
@@ -87,3 +88,5 @@ if(TECA_HAS_PYTHON)
 endif()
 teca_add_python_app(teca_blocking_event_detect_blobs FEATURES ${TECA_HAS_NETCDF})
 teca_add_python_app(teca_blocking_event_create_overlap_table FEATURES ${TECA_HAS_NETCDF})
+teca_add_python_app(teca_blocking_event_tracker)
+teca_add_python_app(teca_blocking_event_relabel FEATURES ${TECA_HAS_NETCDF})

apps/teca_blocking_event_create_overlap_table.in

@@ -825,6 +825,7 @@ def create_argument_parser():
         help=f'background ID value used in component labeling (default:'
              f' {DEFAULT_BACKGROUND_ID})')
 
+    # Verbosity and threading
     parser.add_argument('--verbose', action='store_true',
         help='enable verbose output for debugging')
 

apps/teca_blocking_event_detect_blobs.in

@@ -776,13 +776,15 @@ def create_argument_parser():
         required=False, default='threshold_H',
         help='name of variable with geopotential height threshold (default:'
              'threshold_H)')
+
     parser.add_argument('--x_axis_variable', type=str,
         default=DEFAULT_X_AXIS_VARIABLE, required=False,
         help=f'name of x coordinate variable (default: {DEFAULT_X_AXIS_VARIABLE})')
 
     parser.add_argument('--y_axis_variable', type=str,
         default=DEFAULT_Y_AXIS_VARIABLE, required=False,
         help=f'name of y coordinate variable (default: {DEFAULT_Y_AXIS_VARIABLE})')
+
     parser.add_argument('--z_axis_variable', type=str,
         default='lev', required=False,
         help='name of z coordinate variable (default: lev)')
@@ -817,6 +819,7 @@ def create_argument_parser():
 
     parser.add_argument('--last_step', type=int, default=-1, required=False,
         help='last time step to process. -1 means entire dataset (default: -1)')
+
     parser.add_argument('--calendar', type=str, required=False,
         help='time calendar')
 
@@ -852,6 +855,7 @@ def create_argument_parser():
         help=f'minimum component area in square kilometers for blocking events'
              f' (default: {DEFAULT_MIN_AREA})')
 
+    # Threading options
     parser.add_argument('--n_threads', type=int, default=-1,
         required=False,
         help='number of threads for the writer. -1 means one thread per physical'

apps/teca_blocking_event_operations.py

@@ -256,4 +256,106 @@ def __call__(self, port, data_in, req):
         # Create binary mask: 1 if component is not background, 0 otherwise
         binary_mask = (component_data != self.background_value).astype(np.byte)
         arrays[self.binary_varname] = binary_mask
-        return out_mesh
+        return out_mesh
+
+
+class ReLabelOp:
+    """
+    Operator for relabeling region IDs based on a provided mapping.
+
+    Region IDs are relabeled according to a mapping dictionary that specifies
+    how local IDs should be mapped to global IDs per timestep. If no mapping is
+    provided for a given timestep or local ID, a configurable background value
+    is used.
+
+    Attributes:
+        relabel_varname (str): Name of input variable to relabel.
+        background_value (int): Value used for the background in
+                                both input and output.
+        relabel_mapping (dict): A mapping dictionary with the format
+                                {<timestep>: {<local_id>: <global_id>}}.
+        relabeled_varname (str): Name of the output relabeled variable.
+    """
+
+    def __init__(self, relabel_varname, background_value, relabel_mapping, relabeled_varname):
+        """
+        Initialize the operator with the specified relabel mapping and background value.
+
+        Parameters
+        ----------
+        relabel_varname : str
+            Name of the variable to which relabeling is applied.
+        background_value : int
+            Value used for regions that do not have a corresponding mapping.
+        relabel_mapping : dict
+            Mapping for local ids to global ids per timestep.
+        relabeled_varname : str
+            Name of the output relabeled variable.
+        """
+        self.relabel_varname = relabel_varname
+        self.background_value = background_value
+        self.relabel_mapping = relabel_mapping
+        self.relabeled_varname = relabeled_varname
+
+    def __call__(self, port, data_in, req):
+        """
+        Relabel the region IDs in the mesh according to the mapping.
+
+        Parameters
+        ----------
+        port : int
+            The input port number (not used in this implementation).
+        data_in : list
+            List containing the input mesh data.
+        req : dict
+            Request dictionary that may contain 'device_id' for CUDA execution.
+
+        Returns
+        -------
+        teca_cartesian_mesh
+            A new mesh with relabeled region IDs based on the mapping.
+        """
+        dev = -1
+        np = numpy
+        if teca.get_teca_has_cuda() and teca.get_teca_has_cupy():
+            dev = req.get('device_id', -1)
+            if dev >= 0:
+                cupy.cuda.Device(dev).use()
+                np = cupy
+
+        # Get the current time index
+        index_request_key = req['index_request_key']
+        index_request = req[index_request_key]
+        t_index = index_request[0]
+        # ... double check that only one index was requested
+        if t_index != index_request[1]:
+            raise NotImplementedError('generate_overlap_table::execute: only one'
+                                      ' time index can be requested at a time')
+
+        in_mesh = teca.as_teca_cartesian_mesh(data_in[0])
+        out_mesh = teca.teca_cartesian_mesh.New()
+        out_mesh.shallow_copy(in_mesh)
+        arrays = out_mesh.get_point_arrays()
+
+        # Obtain the region labels array from the input mesh.
+        region_labels_array = arrays[self.relabel_varname]
+
+        if dev < 0:
+            region_labels_array = region_labels_array.get_host_accessible()
+        else:
+            region_labels_array = region_labels_array.get_cuda_accessible()
+
+        # Get mapping for the current timestep or use background_value
+        relabel_mapping_for_timestep = self.relabel_mapping.get(t_index, {})
+
+        # Map the region labels according to the mapping dictionary
+        unique_labels, inv_indices = np.unique(region_labels_array, return_inverse=True)
+        relabeled_ids = np.array(
+            [relabel_mapping_for_timestep.get(label, self.background_value)
+             for label in unique_labels],
+            dtype=np.int32
+            )
+        relabeled_ds = relabeled_ids[inv_indices]
+
+        arrays[self.relabeled_varname] = relabeled_ds
+        return out_mesh