Repository: tobac-project/tobac

Issue:
3D feature detection with min_distance fails if dz is not provided 
If dz is not given, `filter_min_distance` tries to find the vertical coordinate from the dataframe but as this occurs before coordinate interpolate it fails, causing an error.

Easiest solution is to move this after coordinate interpolation I think

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[22], line 1
----> 1 features_w = tobac.feature_detection_multithreshold(
      2     w_iris, **parameters_features
      3 )

File ~/python/tobac/tobac/utils/decorators.py:168, in xarray_to_iris.<locals>.wrapper(*args, **kwargs)
    165             output = xarray.DataArray.from_iris(output)
    167 else:
--> 168     output = func(*args, **kwargs)
    169 # print(output)
    170 return output

File ~/python/tobac/tobac/feature_detection.py:1375, in feature_detection_multithreshold(field_in, dxy, threshold, min_num, target, position_threshold, sigma_threshold, n_erosion_threshold, n_min_threshold, min_distance, feature_number_start, PBC_flag, vertical_coord, vertical_axis, detect_subset, wavelength_filtering, dz, strict_thresholding, statistic)
   1372     # Loop over DataFrame to remove features that are closer than distance_min to each
   1373     # other:
   1374     if min_distance > 0:
-> 1375         features_thresholds = filter_min_distance(
   1376             features_thresholds,
   1377             dxy=dxy,
   1378             dz=dz,
   1379             min_distance=min_distance,
   1380             z_coordinate_name=vertical_coord,
   1381             target=target,
   1382             PBC_flag=PBC_flag,
   1383             min_h1=0,
   1384             max_h1=hdim1_max,
   1385             min_h2=0,
   1386             max_h2=hdim2_max,
   1387         )
   1388 list_features_timesteps.append(features_thresholds)
   1390 logging.debug(
   1391     "Finished feature detection for " + time_i.strftime("%Y-%m-%d_%H:%M:%S")
   1392 )

File ~/python/tobac/tobac/feature_detection.py:1492, in filter_min_distance(features, dxy, dz, min_distance, x_coordinate_name, y_coordinate_name, z_coordinate_name, target, PBC_flag, min_h1, max_h1, min_h2, max_h2)
   1489 is_3D = "vdim" in features
   1491 if is_3D and dz is None:
-> 1492     z_coordinate_name = internal_utils.find_dataframe_vertical_coord(
   1493         features, z_coordinate_name
   1494     )
   1496 # Check if both dxy and their coordinate names are specified.
   1497 # If they are, warn that we will use dxy.
   1498 if dxy is not None and (
   1499     x_coordinate_name in features and y_coordinate_name in features
   1500 ):

File ~/python/tobac/tobac/utils/internal/basic.py:187, in find_dataframe_vertical_coord(variable_dataframe, vertical_coord)
    185     return vertical_coord
    186 else:
--> 187     raise ValueError("Please specify vertical coordinate")

ValueError: Please specify vertical coordinate
```

* [x] Have you searched the issue tracker for the same problem?
* [x] Have you checked if you're using the latest version? If not, which version are you using?
* [x] Have you mentioned the steps to reproduce the issue?
* [x] Have you, if applicable, included error messages?

Hints from issue comments:
@w-k-jones: Resolved by #452 

Solution patch:
diff --git a/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb b/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
index 3ed2bb69..a8e17947 100644
--- a/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
+++ b/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
@@ -15,42 +15,37 @@
   {
    "cell_type": "code",
    "execution_count": 1,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:27.087494Z",
-     "iopub.status.busy": "2024-03-22T17:17:27.087185Z",
-     "iopub.status.idle": "2024-03-22T17:17:31.199088Z",
-     "shell.execute_reply": "2024-03-22T17:17:31.194250Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "# Import libraries\n",
     "import iris\n",
     "import numpy as np\n",
     "import pandas as pd\n",
-    "import xarray as xr\n",
+    "import xarray as xr \n",
     "import matplotlib.pyplot as plt\n",
     "import datetime\n",
     "import shutil\n",
     "from six.moves import urllib\n",
     "from pathlib import Path\n",
-    "\n",
     "%matplotlib inline"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:31.211608Z",
-     "iopub.status.busy": "2024-03-22T17:17:31.210964Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.419147Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.416362Z"
-    }
-   },
+   "metadata": {},
    "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/glade/u/home/kukulies/tobac/tobac/utils/decorators.py:498: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n",
+      "/glade/u/home/kukulies/tobac/tobac/utils/decorators.py:498: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n"
+     ]
+    },
     {
      "name": "stdout",
      "output_type": "stream",
@@ -62,30 +57,21 @@
    "source": [
     "# Import tobac itself\n",
     "import tobac\n",
-    "\n",
-    "print(\"using tobac version\", str(tobac.__version__))"
+    "print('using tobac version', str(tobac.__version__))"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.455941Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.455446Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.460672Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.460009Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "# Disable a few warnings:\n",
     "import warnings\n",
-    "\n",
-    "warnings.filterwarnings(\"ignore\", category=UserWarning, append=True)\n",
-    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning, append=True)\n",
-    "warnings.filterwarnings(\"ignore\", category=FutureWarning, append=True)\n",
-    "warnings.filterwarnings(\"ignore\", category=pd.io.pytables.PerformanceWarning)"
+    "warnings.filterwarnings('ignore', category=UserWarning, append=True)\n",
+    "warnings.filterwarnings('ignore', category=RuntimeWarning, append=True)\n",
+    "warnings.filterwarnings('ignore', category=FutureWarning, append=True)\n",
+    "warnings.filterwarnings('ignore',category=pd.io.pytables.PerformanceWarning)"
    ]
   },
   {
@@ -98,21 +84,14 @@
   {
    "cell_type": "code",
    "execution_count": 4,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.464189Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.463883Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.468048Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.467231Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "# Set up directory to save output:\n",
-    "savedir = Path(\"Save\")\n",
+    "#Set up directory to save output:\n",
+    "savedir=Path(\"Save\")\n",
     "if not savedir.is_dir():\n",
     "    savedir.mkdir()\n",
-    "plot_dir = Path(\"Plot\")\n",
+    "plot_dir=Path(\"Plot\")\n",
     "if not plot_dir.is_dir():\n",
     "    plot_dir.mkdir()"
    ]
@@ -120,115 +99,75 @@
   {
    "cell_type": "code",
    "execution_count": 5,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.471915Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.471614Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.479102Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.478486Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "data_out = Path(\"../\")\n",
+    "data_out=Path('../')\n",
     "# Download the data: This only has to be done once for all tobac examples and can take a while\n",
-    "data_file = list(data_out.rglob(\"data/Example_input_Precip.nc\"))\n",
+    "data_file = list(data_out.rglob('data/Example_input_Precip.nc'))\n",
     "if len(data_file) == 0:\n",
-    "    file_path = \"https://zenodo.org/records/3195910/files/climate-processes/tobac_example_data-v1.0.1.zip\"\n",
-    "    # file_path='http://zenodo..'\n",
-    "    tempfile = Path(\"temp.zip\")\n",
-    "    print(\"start downloading data\")\n",
-    "    request = urllib.request.urlretrieve(file_path, tempfile)\n",
-    "    print(\"start extracting data\")\n",
+    "    file_path='https://zenodo.org/records/3195910/files/climate-processes/tobac_example_data-v1.0.1.zip'\n",
+    "    #file_path='http://zenodo..'\n",
+    "    tempfile=Path('temp.zip')\n",
+    "    print('start downloading data')\n",
+    "    request=urllib.request.urlretrieve(file_path, tempfile)\n",
+    "    print('start extracting data')\n",
     "    shutil.unpack_archive(tempfile, data_out)\n",
     "    tempfile.unlink()\n",
-    "    print(\"data extracted\")\n",
-    "    data_file = list(data_out.rglob(\"data/Example_input_Precip.nc\"))"
+    "    print('data extracted')\n",
+    "    data_file = list(data_out.rglob('data/Example_input_Precip.nc'))"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 6,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.482460Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.482198Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.706441Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.705850Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "Precip = iris.load_cube(str(data_file[0]), \"surface_precipitation_average\")"
+    "Precip=iris.load_cube(str(data_file[0]),'surface_precipitation_average')"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 7,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.709868Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.709531Z",
-     "iopub.status.idle": "2024-03-22T17:17:36.716902Z",
-     "shell.execute_reply": "2024-03-22T17:17:36.715785Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "parameters_features = {}\n",
-    "parameters_features[\"position_threshold\"] = \"weighted_diff\"\n",
-    "parameters_features[\"sigma_threshold\"] = 0.5\n",
-    "parameters_features[\"min_distance\"] = 0\n",
-    "parameters_features[\"sigma_threshold\"] = 1\n",
-    "parameters_features[\"threshold\"] = [1, 2, 3, 4, 5, 10, 15]  # mm/h\n",
-    "parameters_features[\"n_erosion_threshold\"] = 0\n",
-    "parameters_features[\"n_min_threshold\"] = 3\n",
+    "parameters_features={}\n",
+    "parameters_features['position_threshold']='weighted_diff'\n",
+    "parameters_features['sigma_threshold']=0.5\n",
+    "parameters_features['min_distance']=0\n",
+    "parameters_features['sigma_threshold']=1\n",
+    "parameters_features['threshold']=[1,2,3,4,5,10,15] #mm/h\n",
+    "parameters_features['n_erosion_threshold']=0\n",
+    "parameters_features['n_min_threshold']=3\n",
     "\n",
     "# get temporal and spation resolution of the data\n",
-    "dxy, dt = tobac.get_spacings(Precip)"
+    "dxy,dt=tobac.get_spacings(Precip)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 8,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:36.720988Z",
-     "iopub.status.busy": "2024-03-22T17:17:36.720671Z",
-     "iopub.status.idle": "2024-03-22T17:17:42.032429Z",
-     "shell.execute_reply": "2024-03-22T17:17:42.031730Z"
-    }
-   },
+   "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "starting feature detection based on multiple thresholds\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "feature detection done\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
+      "starting feature detection based on multiple thresholds\n",
+      "feature detection done\n",
       "features saved\n"
      ]
     }
    ],
    "source": [
     "# Feature detection based on based on surface precipitation field and a range of thresholds\n",
-    "print(\"starting feature detection based on multiple thresholds\")\n",
-    "Features = tobac.feature_detection_multithreshold(Precip, dxy, **parameters_features)\n",
-    "print(\"feature detection done\")\n",
-    "Features.to_hdf(savedir / \"Features.h5\", \"table\")\n",
-    "print(\"features saved\")"
+    "print('starting feature detection based on multiple thresholds')\n",
+    "Features= tobac.feature_detection_multithreshold(Precip,dxy,**parameters_features) \n",
+    "print('feature detection done')\n",
+    "Features.to_hdf(savedir / 'Features.h5','table')\n",
+    "print('features saved')"
    ]
   },
   {
@@ -241,73 +180,41 @@
   {
    "cell_type": "code",
    "execution_count": 9,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:42.036502Z",
-     "iopub.status.busy": "2024-03-22T17:17:42.035726Z",
-     "iopub.status.idle": "2024-03-22T17:17:42.045228Z",
-     "shell.execute_reply": "2024-03-22T17:17:42.042712Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "# Dictionary containing keyword arguments for segmentation step:\n",
-    "parameters_segmentation = {}\n",
-    "parameters_segmentation[\"method\"] = \"watershed\"\n",
-    "parameters_segmentation[\"threshold\"] = 1  # mm/h mixing ratio\n",
+    "parameters_segmentation={}\n",
+    "parameters_segmentation['method']='watershed'\n",
+    "parameters_segmentation['threshold']=1  # mm/h mixing ratio\n",
     "\n",
     "# get temporal and spation resolution of the data\n",
-    "dxy, dt = tobac.get_spacings(Precip)"
+    "dxy,dt=tobac.get_spacings(Precip)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 10,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:42.049101Z",
-     "iopub.status.busy": "2024-03-22T17:17:42.048656Z",
-     "iopub.status.idle": "2024-03-22T17:17:52.522098Z",
-     "shell.execute_reply": "2024-03-22T17:17:52.514694Z"
-    }
-   },
+   "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Starting segmentation based on surface precipitation\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "segmentation based on surface precipitation performed, start saving results to files\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
+      "Starting segmentation based on surface precipitation\n",
+      "segmentation based on surface precipitation performed, start saving results to files\n",
       "segmentation surface precipitation performed and saved\n"
      ]
     }
    ],
    "source": [
     "# Perform Segmentation and save resulting mask to NetCDF file:\n",
-    "print(\"Starting segmentation based on surface precipitation\")\n",
-    "Mask_Precip, Features_Precip = tobac.segmentation_2D(\n",
-    "    Features, Precip, dxy, **parameters_segmentation\n",
-    ")\n",
-    "print(\n",
-    "    \"segmentation based on surface precipitation performed, start saving results to files\"\n",
-    ")\n",
-    "iris.save(\n",
-    "    [Mask_Precip], savedir / \"Mask_segmentation_precip.nc\", zlib=True, complevel=4\n",
-    ")\n",
-    "Features_Precip.to_hdf(savedir / \"Features_Precip.h5\", \"table\")\n",
-    "print(\"segmentation surface precipitation performed and saved\")"
+    "print('Starting segmentation based on surface precipitation')\n",
+    "Mask_Precip,Features_Precip=tobac.segmentation_2D(Features,Precip,dxy,**parameters_segmentation)\n",
+    "print('segmentation based on surface precipitation performed, start saving results to files')\n",
+    "iris.save([Mask_Precip], savedir / 'Mask_segmentation_precip.nc', zlib=True, complevel=4)                \n",
+    "Features_Precip.to_hdf(savedir / 'Features_Precip.h5', 'table')\n",
+    "print('segmentation surface precipitation performed and saved')"
    ]
   },
   {
@@ -322,14 +229,7 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:52.528489Z",
-     "iopub.status.busy": "2024-03-22T17:17:52.527447Z",
-     "iopub.status.idle": "2024-03-22T17:17:52.538706Z",
-     "shell.execute_reply": "2024-03-22T17:17:52.535020Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "from tobac.utils import get_statistics_from_mask"
@@ -352,20 +252,13 @@
   {
    "cell_type": "code",
    "execution_count": 12,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:52.548583Z",
-     "iopub.status.busy": "2024-03-22T17:17:52.548132Z",
-     "iopub.status.idle": "2024-03-22T17:17:52.554124Z",
-     "shell.execute_reply": "2024-03-22T17:17:52.553158Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "statistics = {}\n",
-    "statistics[\"mean_precip\"] = np.mean\n",
-    "statistics[\"total_precip\"] = np.sum\n",
-    "statistics[\"max_precip\"] = np.max"
+    "statistics['mean_precip'] = np.mean\n",
+    "statistics['total_precip'] = np.sum\n",
+    "statistics['max_precip'] = np.max"
    ]
   },
   {
@@ -378,35 +271,19 @@
   {
    "cell_type": "code",
    "execution_count": 13,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:52.562198Z",
-     "iopub.status.busy": "2024-03-22T17:17:52.561074Z",
-     "iopub.status.idle": "2024-03-22T17:17:52.571050Z",
-     "shell.execute_reply": "2024-03-22T17:17:52.569254Z"
-    }
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "statistics[\"percentiles\"] = (np.percentile, {\"q\": [95, 99]})"
+    "statistics['percentiles'] = (np.percentile, {'q': [95,99]})"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:52.575766Z",
-     "iopub.status.busy": "2024-03-22T17:17:52.575276Z",
-     "iopub.status.idle": "2024-03-22T17:17:57.003737Z",
-     "shell.execute_reply": "2024-03-22T17:17:57.001134Z"
-    }
-   },
+   "execution_count": 15,
+   "metadata": {},
    "outputs": [],
    "source": [
-    "features_with_stats = get_statistics_from_mask(\n",
-    "    Features_Precip, Mask_Precip, Precip, statistic=statistics\n",
-    ")"
+    "features_with_stats = get_statistics_from_mask(Features, Mask_Precip, Precip, statistic=statistics)"
    ]
   },
   {
@@ -418,15 +295,8 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:57.007064Z",
-     "iopub.status.busy": "2024-03-22T17:17:57.006774Z",
-     "iopub.status.idle": "2024-03-22T17:17:57.013033Z",
-     "shell.execute_reply": "2024-03-22T17:17:57.012425Z"
-    }
-   },
+   "execution_count": 16,
+   "metadata": {},
    "outputs": [
     {
      "data": {
@@ -439,7 +309,7 @@
        "Name: mean_precip, dtype: object"
       ]
      },
-     "execution_count": 15,
+     "execution_count": 16,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -450,28 +320,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:57.015793Z",
-     "iopub.status.busy": "2024-03-22T17:17:57.015558Z",
-     "iopub.status.idle": "2024-03-22T17:17:57.020851Z",
-     "shell.execute_reply": "2024-03-22T17:17:57.020126Z"
-    }
-   },
+   "execution_count": 17,
+   "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "0    16.296951\n",
+       "0    16.296949\n",
        "1    14.095468\n",
        "2    26.856783\n",
        "3    36.831512\n",
-       "4    49.737709\n",
+       "4    49.737713\n",
        "Name: total_precip, dtype: object"
       ]
      },
-     "execution_count": 16,
+     "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -482,15 +345,8 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-03-22T17:17:57.023633Z",
-     "iopub.status.busy": "2024-03-22T17:17:57.023395Z",
-     "iopub.status.idle": "2024-03-22T17:17:57.028693Z",
-     "shell.execute_reply": "2024-03-22T17:17:57.028166Z"
-    }
-   },
+   "execution_count": 18,
+   "metadata": {},
    "outputs": [
     {
      "data": {
@@ -503,7 +359,7 @@
        "Name: percentiles, dtype: object"
       ]
      },
-     "execution_count": 17,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -521,6 +377,11 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "tobac-interactive",
+   "language": "python",
+   "name": "tobac-interactive"
+  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
@@ -531,7 +392,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,
diff --git a/tobac/feature_detection.py b/tobac/feature_detection.py
index 075b8513..a866e897 100644
--- a/tobac/feature_detection.py
+++ b/tobac/feature_detection.py
@@ -1121,39 +1121,34 @@ def feature_detection_multithreshold_timestep(
         elif i_threshold == 0:
             regions_old = regions_i
 
-        if statistic:
-            # reconstruct the labeled regions based on the regions dict
-            labels = np.zeros(track_data.shape)
-            labels = labels.astype(int)
-            for key in regions_old.keys():
-                labels.ravel()[regions_old[key]] = key
-                # apply function to get statistics based on labeled regions and functions provided by the user
-                # the feature dataframe is updated by appending a column for each metric
-            if statistics_unsmoothed:
-                features_thresholds = get_statistics(
-                    features_thresholds,
-                    labels,
-                    data_i.core_data(),
-                    statistic=statistic,
-                    index=np.unique(labels[labels > 0]),
-                    id_column="idx",
-                )
-            else:
-                features_thresholds = get_statistics(
-                    features_thresholds,
-                    labels,
-                    track_data,
-                    statistic=statistic,
-                    index=np.unique(labels[labels > 0]),
-                    id_column="idx",
-                )
-
         logging.debug(
             "Finished feature detection for threshold "
             + str(i_threshold)
             + " : "
             + str(threshold_i)
         )
+
+    if statistic:
+        # reconstruct the labeled regions based on the regions dict
+        labels = np.zeros(track_data.shape)
+        labels = labels.astype(int)
+        for key in regions_old.keys():
+            labels.ravel()[regions_old[key]] = key
+            # apply function to get statistics based on labeled regions and functions provided by the user
+            # the feature dataframe is updated by appending a column for each metric
+        
+        # select which data to use according to statistics_unsmoothed option
+        stats_data = data_i.core_data() if statistics_unsmoothed else track_data
+        
+        features_thresholds = get_statistics(
+            features_thresholds,
+            labels,
+            stats_data, 
+            statistic=statistic,
+            index=np.unique(labels[labels > 0]),
+            id_column="idx",
+        )
+
     return features_thresholds
 
 
@@ -1296,6 +1291,9 @@ def feature_detection_multithreshold(
     if detect_subset is not None and ndim_time in detect_subset:
         raise NotImplementedError("Cannot subset on time")
 
+    # Remember if dz is set and not vertical coord for min distance filtering
+    use_dz_for_filtering = dz is not None
+
     if is_3D:
         # We need to determine the time axis so that we can determine the
         # vertical axis in each timestep if vertical_axis is not none.
@@ -1395,30 +1393,7 @@ def feature_detection_multithreshold(
             statistic=statistic,
             statistics_unsmoothed=statistics_unsmoothed,
         )
-        # check if list of features is not empty, then merge features from different threshold
-        # values into one DataFrame and append to list for individual timesteps:
-        if not features_thresholds.empty:
-            hdim1_ax, hdim2_ax = internal_utils.find_hdim_axes_3D(
-                field_in, vertical_coord=vertical_coord
-            )
-            hdim1_max = field_in.shape[hdim1_ax] - 1
-            hdim2_max = field_in.shape[hdim2_ax] - 1
-            # Loop over DataFrame to remove features that are closer than distance_min to each
-            # other:
-            if min_distance > 0:
-                features_thresholds = filter_min_distance(
-                    features_thresholds,
-                    dxy=dxy,
-                    dz=dz,
-                    min_distance=min_distance,
-                    z_coordinate_name=vertical_coord,
-                    target=target,
-                    PBC_flag=PBC_flag,
-                    min_h1=0,
-                    max_h1=hdim1_max,
-                    min_h2=0,
-                    max_h2=hdim2_max,
-                )
+
         list_features_timesteps.append(features_thresholds)
 
         logging.debug(
@@ -1440,9 +1415,41 @@ def feature_detection_multithreshold(
             )
         else:
             features = add_coordinates(features, field_in)
+
+        # Loop over DataFrame to remove features that are closer than distance_min to each
+        # other:
+        filtered_features = []
+        if min_distance > 0:
+            hdim1_ax, hdim2_ax = internal_utils.find_hdim_axes_3D(
+                field_in, vertical_coord=vertical_coord
+            )
+            hdim1_max = field_in.shape[hdim1_ax] - 1
+            hdim2_max = field_in.shape[hdim2_ax] - 1
+
+            for _, features_frame in features.groupby("frame"):
+                filtered_features.append(
+                    filter_min_distance(
+                        features_frame,
+                        dxy=dxy,
+                        dz=dz if use_dz_for_filtering else None,
+                        min_distance=min_distance,
+                        z_coordinate_name=(
+                            None if use_dz_for_filtering else vertical_coord
+                        ),
+                        target=target,
+                        PBC_flag=PBC_flag,
+                        min_h1=0,
+                        max_h1=hdim1_max,
+                        min_h2=0,
+                        max_h2=hdim2_max,
+                    )
+                )
+            features = pd.concat(filtered_features, ignore_index=True)
+
     else:
         features = None
         logging.debug("No features detected")
+
     logging.debug("feature detection completed")
     return features
 
@@ -1512,47 +1519,47 @@ def filter_min_distance(
     pandas DataFrame
         features after filtering
     """
+    # Optional coordinate names are not yet implemented, set to defaults here:
     if dxy is None:
         raise NotImplementedError("dxy currently must be set.")
 
-    # if PBC_flag != "none":
-    #    raise NotImplementedError("We haven't yet implemented PBCs into this.")
-
-    # if we are 3D, the vertical dimension is in features. if we are 2D, there
-    # is no vertical dimension in features.
-    is_3D = "vdim" in features
-
-    if is_3D and dz is None:
-        z_coordinate_name = internal_utils.find_dataframe_vertical_coord(
-            features, z_coordinate_name
-        )
-
     # Check if both dxy and their coordinate names are specified.
     # If they are, warn that we will use dxy.
-    if dxy is not None and (
-        x_coordinate_name in features and y_coordinate_name in features
-    ):
+    elif x_coordinate_name in features and y_coordinate_name in features:
         warnings.warn(
             "Both " + x_coordinate_name + "/" + y_coordinate_name + " and dxy "
             "set. Using constant dxy. Set dxy to None if you want to use the "
             "interpolated coordinates, or set `x_coordinate_name` and "
             "`y_coordinate_name` to None to use a constant dxy."
         )
+        y_coordinate_name = "hdim_1"
+        x_coordinate_name = "hdim_2"
+    # If dxy only, use hdim_1, hdim_1 as default horizontal dimensions
+    else:
+        y_coordinate_name = "hdim_1"
+        x_coordinate_name = "hdim_2"
 
-    # Check and if both dz is specified and altitude is available, warn that we will use dz.
-    if is_3D and (dz is not None and z_coordinate_name in features):
-        warnings.warn(
-            "Both "
-            + z_coordinate_name
-            + " and dz available to filter_min_distance; using constant dz. "
-            "Set dz to none if you want to use altitude or set `z_coordinate_name` to None to use "
-            "constant dz."
-        )
-
-    # As optional coordinate names are not yet implemented, set to defaults here:
-    z_coordinate_name = "vdim"
-    y_coordinate_name = "hdim_1"
-    x_coordinate_name = "hdim_2"
+    # if we are 3D, the vertical dimension is in features
+    is_3D = "vdim" in features
+    if is_3D:
+        if dz is None:
+            # Find vertical coord name and set dz to 1
+            z_coordinate_name = internal_utils.find_dataframe_vertical_coord(
+                variable_dataframe=features, vertical_coord=z_coordinate_name
+            )
+            dz = 1
+        else:
+            # Use dz, warn if both are set
+            if z_coordinate_name is not None:
+                warnings.warn(
+                    "Both "
+                    + z_coordinate_name
+                    + " and dz available to filter_min_distance; using constant dz. "
+                    "Set dz to none if you want to use altitude or set `z_coordinate_name` to None to use "
+                    "constant dz.",
+                    UserWarning,
+                )
+            z_coordinate_name = "vdim"
 
     if target not in ["minimum", "maximum"]:
         raise ValueError(
diff --git a/tobac/utils/bulk_statistics.py b/tobac/utils/bulk_statistics.py
index 2dbcf6f2..f9266dd8 100644
--- a/tobac/utils/bulk_statistics.py
+++ b/tobac/utils/bulk_statistics.py
@@ -18,14 +18,13 @@
     from numpy.core import multiarray as mu
 import pandas as pd
 import xarray as xr
-
 from tobac.utils import decorators
 
 
 def get_statistics(
     features: pd.DataFrame,
     labels: np.ndarray[int],
-    *fields: tuple[np.ndarray],
+    *fields: tuple[xr.DataArray],
     statistic: dict[str, Union[Callable, tuple[Callable, dict]]] = {
         "ncells": np.count_nonzero
     },
@@ -54,7 +53,7 @@ def get_statistics(
         Mask with labels of each regions to apply function to (e.g. output of
         segmentation for a specific timestep)
 
-    *fields : tuple[np.ndarray]
+    *fields : tuple[xr.DataArray]
         Fields to give as arguments to each function call. If the shape does not
         match that of labels, numpy-style broadcasting will be applied.
 
@@ -215,7 +214,7 @@ def get_statistics(
 def get_statistics_from_mask(
     features: pd.DataFrame,
     segmentation_mask: xr.DataArray,
-    *fields: xr.DataArray,
+    *fields: tuple[xr.DataArray],
     statistic: dict[str, tuple[Callable]] = {"Mean": np.mean},
     index: Union[None, list[int]] = None,
     default: Union[None, float] = None,
@@ -236,7 +235,7 @@ def get_statistics_from_mask(
     segmentation_mask : xr.DataArray
         Segmentation mask output
 
-    *fields : xr.DataArray[np.ndarray]
+    *fields : tuple[xr.DataArray]
         Field(s) with input data. If field does not have a time dimension it
         will be considered time invariant, and the entire field will be passed
         for each time step in segmentation_mask. If the shape does not match
@@ -291,26 +290,30 @@ def get_statistics_from_mask(
     else:
         collapse_axis = None
 
+    # check if any of the feature dataframe input values match with segmentaion mask IDs
+    if not np.any(np.isin(features[id_column], np.unique(segmentation_mask))):
+        raise ValueError(
+            "The labels of the segmentation mask and the feature dataframe do not seem to match. Please make sure you provide the correct input feature dataframe  to calculate the bulk statistics."
+        )
+
     # get bulk statistics for each timestep
     step_statistics = []
 
     for tt in pd.to_datetime(segmentation_mask.time):
         # select specific timestep
-        segmentation_mask_t = segmentation_mask.sel(time=tt).data
+        segmentation_mask_t = segmentation_mask.sel(time=tt, method = 'nearest').data
         fields_t = (
-            field.sel(time=tt).values if "time" in field.coords else field.values
+            field.sel(time=tt, method = 'nearest', tolerance = np.timedelta64(1000, 'us')).values if "time" in field.coords else field.values
             for field in fields
         )
-        features_t = features.loc[features.time == tt].copy()
 
+        features_t = features.loc[features.time == tt].copy()
         # make sure that the labels in the segmentation mask exist in feature dataframe
-        if (
-            np.intersect1d(np.unique(segmentation_mask_t), features_t[id_column]).size
-            > np.unique(segmentation_mask_t).size
-        ):
-            raise ValueError(
-                "The labels of the segmentation mask and the feature dataframe do not seem to match. Please make sure you provide the correct input feature dataframe to calculate the bulk statistics. "
-            )
+        # continue loop because not all timesteps might have matching IDs
+        if not np.any(np.isin(features_t[id_column], np.unique(segmentation_mask_t))):
+            warnings.warn("Not all timesteps have matching features", UserWarning)
+            step_statistics.append(features_t)
+            continue
         else:
             # make sure that features are not double-defined
             step_statistics.append(


Test patch:
diff --git a/tobac/tests/test_feature_detection.py b/tobac/tests/test_feature_detection.py
index 37b6b586..d872e2d8 100644
--- a/tobac/tests/test_feature_detection.py
+++ b/tobac/tests/test_feature_detection.py
@@ -319,6 +319,34 @@ def test_feature_detection_position(position_threshold):
             True,
             False,
         ),
+        (  # Test using z coord name
+            (0, 0, 0, 4, 1),
+            (1, 1, 1, 4, 1),
+            1000,
+            None,
+            1,
+            "maximum",
+            False,
+            False,
+            True,
+            "none",
+            True,
+            True,
+        ),
+        (  # Test using z coord name
+            (0, 0, 0, 5, 1),
+            (1, 1, 1, 4, 1),
+            1,
+            None,
+            101,
+            "maximum",
+            False,
+            False,
+            True,
+            "none",
+            True,
+            False,
+        ),
     ],
 )
 def test_filter_min_distance(
@@ -415,18 +443,6 @@ def test_filter_min_distance(
 
     feat_combined = pd.concat([feat_1_interp, feat_2_interp], ignore_index=True)
 
-    filter_dist_opts = dict()
-
-    if add_x_coords:
-        feat_combined[x_coord_name] = feat_combined["hdim_2"] * assumed_dxy
-        filter_dist_opts["x_coordinate_name"] = x_coord_name
-    if add_y_coords:
-        feat_combined[y_coord_name] = feat_combined["hdim_1"] * assumed_dxy
-        filter_dist_opts["y_coordinate_name"] = y_coord_name
-    if add_z_coords and is_3D:
-        feat_combined[z_coord_name] = feat_combined["vdim"] * assumed_dz
-        filter_dist_opts["z_coordinate_name"] = z_coord_name
-
     filter_dist_opts = {
         "features": feat_combined,
         "dxy": dxy,
@@ -439,6 +455,16 @@ def test_filter_min_distance(
         "min_h2": 0,
         "max_h2": 100,
     }
+    if add_x_coords:
+        feat_combined[x_coord_name] = feat_combined["hdim_2"] * assumed_dxy
+        filter_dist_opts["x_coordinate_name"] = x_coord_name
+    if add_y_coords:
+        feat_combined[y_coord_name] = feat_combined["hdim_1"] * assumed_dxy
+        filter_dist_opts["y_coordinate_name"] = y_coord_name
+    if add_z_coords and is_3D:
+        feat_combined[z_coord_name] = feat_combined["vdim"] * assumed_dz
+        filter_dist_opts["z_coordinate_name"] = z_coord_name
+
     if target not in ["maximum", "minimum"]:
         with pytest.raises(ValueError):
             out_feats = feat_detect.filter_min_distance(**filter_dist_opts)
diff --git a/tobac/tests/test_utils_bulk_statistics.py b/tobac/tests/test_utils_bulk_statistics.py
index c62ee821..2026f9c4 100644
--- a/tobac/tests/test_utils_bulk_statistics.py
+++ b/tobac/tests/test_utils_bulk_statistics.py
@@ -40,7 +40,12 @@ def test_bulk_statistics_fd(statistics_unsmoothed):
 
 
 @pytest.mark.parametrize(
-    "id_column, index", [("feature", [1]), ("feature_id", [1]), ("cell", [1])]
+    "id_column, index",
+    [
+        ("feature", [1]),
+        ("feature_id", [1]),
+        ("cell", [1]),
+    ],
 )
 def test_bulk_statistics(id_column, index):
     """
@@ -60,7 +65,6 @@ def test_bulk_statistics(id_column, index):
 
     # detect features
     threshold = 7
-    # test_data_iris = testing.make_dataset_from_arr(test_data, data_type="iris")
     fd_output = tobac.feature_detection.feature_detection_multithreshold(
         test_data_iris,
         dxy=1000,
@@ -81,6 +85,7 @@ def test_bulk_statistics(id_column, index):
 
     #### checks
     out_df = out_df.rename(columns={"feature": id_column})
+
     #  assure that bulk statistics in postprocessing give same result
     out_segmentation = tb_utils.get_statistics_from_mask(
         out_df, out_seg_mask, test_data_iris, statistic=stats, id_column=id_column
@@ -139,6 +144,60 @@ def test_bulk_statistics(id_column, index):
         )
 
 
+def test_bulk_statistics_missing_segments():
+    """
+    Test that output feature dataframe contains all the same time steps even though for some timesteps,
+    the statistics have not been calculated (in the case of unmatching labels or no segment labels for a given feature)
+    """
+
+    ### Test 2D data with time dimension
+    test_data = tb_test.make_simple_sample_data_2D().core_data()
+    common_dset_opts = {
+    "in_arr": test_data,
+    "data_type": "iris",}
+    
+    test_data_iris = tb_test.make_dataset_from_arr(
+    time_dim_num=0, y_dim_num=1, x_dim_num=2, **common_dset_opts)
+
+    # detect features
+    threshold = 7
+    # test_data_iris = testing.make_dataset_from_arr(test_data, data_type="iris")
+    fd_output = tobac.feature_detection.feature_detection_multithreshold(
+    test_data_iris,
+    dxy=1000,
+    threshold=[threshold],
+    n_min_threshold=100,
+    target="maximum",)
+
+    # perform segmentation with bulk statistics
+    stats = {
+    "segment_max": np.max,
+    "segment_min": min,
+    "percentiles": (np.percentile, {"q": 95}),}
+
+    out_seg_mask, out_df = tobac.segmentation.segmentation_2D(
+    fd_output, test_data_iris, dxy=1000, threshold=threshold)
+
+    # specify some timesteps we set to zero 
+    timesteps_to_zero = [1, 3, 10]  # 0-based indexing
+    modified_data = out_seg_mask.data.copy() 
+    # Set values to zero for the specified timesteps
+    for timestep in timesteps_to_zero:
+        modified_data[timestep, :, :] = 0  # Set all values for this timestep to zero
+
+    #  assure that bulk statistics in postprocessing give same result
+    out_segmentation = tb_utils.get_statistics_from_mask(
+    out_df, out_seg_mask, test_data_iris, statistic=stats)
+
+    assert out_df.time.unique().size == out_segmentation.time.unique().size
+
+
+
+
+    
+
+
+        
 def test_bulk_statistics_multiple_fields():
     """
     Test that multiple field input to bulk_statistics works as intended