Repository: tobac-project/tobac

Issue:
Floating point error in feature detection breaks interpn
It looks like there is a chance during the feature detection step for the identification of cells on boundary edges which result in hdim1 or hdim2 values slightly outside of the maximum boundary. I.e. the feature detection can create a cell with ix value of 218.00000003 when the maximum index is 218, and this throws an error in tobac/general.py:107 with `coordinate_points = interpn(points, values, xi)` raising `ValueError: One of the requested xi is out of bounds in dimension 0 `

Full Error:
```
Traceback (most recent call last):

  Cell In[282], line 1
    wrf_features = wrf_tobac_feature_id(wrf_cube, 'IC', CONFIG)

  File ~/CoMET/CoMET/wrf_tobac.py:40 in wrf_tobac_feature_id
    wrf_radar_features = tobac.feature_detection.feature_detection_multithreshold(cube, dxy=dxy, **CONFIG['wrf']['tobac']['feature_id'])

  File ~/.conda/envs/CoMET/lib/python3.11/site-packages/tobac/utils/decorators.py:271 in wrapper
    output = func(*args, **kwargs)

  File ~/.conda/envs/CoMET/lib/python3.11/site-packages/tobac/feature_detection.py:1409 in feature_detection_multithreshold
    features = add_coordinates(features, field_in)

  File ~/.conda/envs/CoMET/lib/python3.11/site-packages/tobac/utils/general.py:109 in add_coordinates
    coordinate_points = interpn(points, values, xi)

  File ~/.conda/envs/CoMET/lib/python3.11/site-packages/scipy/interpolate/_rgi.py:740 in interpn
    raise ValueError("One of the requested xi is out of bounds "

ValueError: One of the requested xi is out of bounds in dimension 0
```


This can be fixed by using np.clip to ensure there is no overflow:
```

elif variable_cube.coord(coord).ndim == 2:
            if variable_cube.coord_dims(coord) == (hdim_1, hdim_2):
                points = (dimvec_1, dimvec_2)
                values = variable_cube.coord(coord).points
                temp_hdim1 = np.clip(t['hdim_1'].values, 0, points[0].max())
                temp_hdim2 = np.clip(t['hdim_2'].values, 0, points[1].max())
                xi = np.column_stack((temp_hdim1, temp_hdim2))
                coordinate_points = interpn(points, values, xi)
```

Hints from issue comments:
@w-k-jones: Hi @travis-j-hahn , thanks for raising this issue. Do you have a minimal example that can recreate this error? (e.g. a subset of the data you were using). I'll have a look into this@travis-j-hahn: [wrf_tb_tobac_test.zip](https://github.com/user-attachments/files/15956912/wrf_tb_tobac_test.zip)
Hi @w-k-jones, I've attached a zip with a small iris cube here containing the problematic values; it is only 2mb. Here are the feature detection parameters I used to see if you can recreate the issue:
```

threshold: [250,225,200,175,150]
target: 'minimum'
position_threshold: 'weighted_diff'
sigma_threshold: 0.5
n_min_threshold: 4
```

Then I used the snippet:

`wrf_radar_features = tobac.feature_detection.feature_detection_multithreshold(test[0], dxy=1000, **params)` where test[0] is the loaded iris cube. Hope this helps!@freemansw1: Anecdotally, I have also seen this issue, and I agree with the proposed solution. Thanks much @travis-j-hahn @w-k-jones: Thanks for sending the test data, I've successfully recreated the error. One thing to watch out for with the fix is that with periodic boundary conditions there are expected situations where the feature point lies outside the range of the coordinate axis (e.g. a feature a 359.5 longitude when the coord ranges 0-359), which we should ensure is handled correctly@w-k-jones: @freemansw1 we will also have to update the new `get_coordinates` functions in #354 @freemansw1: > @freemansw1 we will also have to update the new `get_coordinates` functions in #354

Yes those need fixing anyway. Good note. @freemansw1: > Thanks for sending the test data, I've successfully recreated the error. One thing to watch out for with the fix is that with periodic boundary conditions there are expected situations where the feature point lies outside the range of the coordinate axis (e.g. a feature a 359.5 longitude when the coord ranges 0-359), which we should ensure is handled correctly

Yes, I believe our expected out would be 359.5 there, but I'm not sure how `add_coordinates` handles this with PBCs, not the least of which including the changes in #354. @w-k-jones: Add coordinates doesn't currently handle PBCs at all 😬 I have no idea how we've avoided these errors before now...@w-k-jones: Fixed with #434 

Solution patch:
diff --git a/.github/workflows/codecov-CI.yml b/.github/workflows/codecov-CI.yml
index ef2ef259..c20c96c5 100644
--- a/.github/workflows/codecov-CI.yml
+++ b/.github/workflows/codecov-CI.yml
@@ -27,7 +27,7 @@ jobs:
           mamba install --quiet --yes --file requirements.txt coverage pytest-cov &&
           python -m coverage run -m pytest --cov=./ --cov-report=xml
       - name: Upload Coverage to Codecov
-        uses: codecov/codecov-action@v3
+        uses: codecov/codecov-action@v4
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           fail_ci_if_error: true
diff --git a/.zenodo.json b/.zenodo.json
index 8d6ff6c8..65805358 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -23,7 +23,7 @@
 
         {
             "name": "Kukulies, Julia",
-            "affiliation": "University of Gothenburg (Sweden)",
+            "affiliation": "NSF National Center for Atmospheric Research",
             "orcid": "0000-0001-6084-0069"
         },
         {
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 6664863a..98aee8b8 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -2,16 +2,14 @@
 
 __Welcome! We are very happy that you are interested in our project and thanks for taking time to contribute! :)__
 
-We are currently reorganizing our project. So, please check our modifications later.
-
 
 ## Getting Started
 ### Installation & Environment details
-You will find them in the [README.md](https://github.com/climate-processes/tobac/blob/master/README.md).
+You will find them in the [README.md](https://github.com/tobac-project/tobac/blob/master/README.md).
 
 ### Tutorials
-Tutorials have been prepared to provide you further inside to `tobac`s functionality. Please visit the separate [tobac-tutorials](https://github.com/climate-processes/tobac-tutorials) repository here on Github.
-
+Tutorials have been prepared to provide you further inside to `tobac`s functionality. Please have a look in the 
+[examples folder](https://github.com/tobac-project/tobac/tree/main/examples). 
 
 ### Documentation
 You will find our documentation at [https://tobac.readthedocs.io](https://tobac.readthedocs.io).
@@ -19,56 +17,33 @@ You will find our documentation at [https://tobac.readthedocs.io](https://tobac.
 ### Testing
 The tests are located in the [tests folder](https://github.com/climate-processes/tobac/tree/master/tobac/tests).
 
-
 ## Reporting Bugs
-Please create a new issue on [GitHub](https://github.com/climate-processes/tobac/issues) if it is not listed there, yet.
+Please create a new issue on [GitHub](https://github.com/tobac-project/tobac/issues) if it is not listed there, yet.
 
 ### How to write a good Bug Report?
 * Give it a clear descriptive title.
 * Copy and paste the error message.
 * Describe the steps for reproducing the problem and give an specific example.  
-* Optional: Make a suggestion to fix it.
-
+* Optional: Make a suggestion to fix it. 
 
 ## How to Submit Changes
-* Please read the [README.md](https://github.com/climate-processes/tobac/blob/master/README.md) first, to learn about our project goals and check the [changelog.md]().
-* Before you start a pull request, please make sure that you added [numpydoc docstrings](#docstringExample) to your functions. This way the api documentation will be parsed properly.
-* If it is a larger change or an newly added feature or workflow, please place an example of use in the [tobac-tutorials](https://github.com/climate-processes/tobac-tutorials) repository or adapt the existing examples there.
-* If necessary add a folder or modify a file.
+* Have a look at [our roadmap](https://github.com/tobac-project/tobac-roadmap/blob/master/tobac-roadmap-main.md) first, 
+to learn about our project goals and check the 
+[changelog.md](https://github.com/tobac-project/tobac/blob/main/CHANGELOG.md).
+* More details on the code structure and further help for code contributions can be found in our [developer 
+guide](https://tobac.readthedocs.io/code_structure.html)
+* Before you start a pull request, please make sure that you added [numpydoc 
+docstrings](https://numpydoc.readthedocs.io/en/latest/format.html) to your 
+functions. See [docstring example in the developer guide](https://tobac.readthedocs.io/contributing.html). This way the 
+api documentation will be parsed properly.
+* If it is a larger change or an newly added feature or workflow, please add an example in the [example 
+folder](https://github.com/tobac-project/tobac/tree/main/examples) or adapt the existing examples there.
 * The code should be PEP 8 compliant, as this facilitates our collaboration. Please use the first stable version (22.6.0) of [black](https://black.readthedocs.io/en/stable/) to format your code. When you submit a pull request, all files are checked for formatting.
 * The tobac repository is set up with pre-commit hooks to automatically format your code when commiting changes. Please run the command "pre-commit install" in the root directory of tobac to set up pre-commit formatting.
 
-We aim to respond to all new issues/pull requests as soon as possible, however at times this is not possible due to work commitments.
-
-### Numpydoc Example <a name='docstringExample'>
-```python
-  
-   '''
-   calculate centre of gravity and mass forech individual tracked cell in the simulation
-
-
-    Parameters
-    ----------
-    tracks : pandas.DataFram
-        DataFrame containing trajectories of cell centres
-        
-    param mass : iris.cube.Cube
-        cube of quantity (need coordinates 'time', 'geopotential_height','projection_x_coordinate' and 
-        'projection_y_coordinate')
-        
-    param mask : iris.cube.Cube
-        cube containing mask (int > where belonging to cloud volume, 0 everywhere else )
-
+We aim to respond to all new issues/pull requests as soon as possible, however sometimes this is not possible due to work commitments.
 
-    Returns
-    -------
-    track_out : pandas.DataFrame
-        Dataframe containing t,x,y,z positions of centre of gravity and total cloud mass each tracked cells 
-        at each timestep
-    
-    '''
 
-```
 
 ## Slack
 In addition to the workflow here on Github, there's a tobac workspace on Slack [tobac-dev.slack.com](tobac-dev.slack.com) that we use for some additional communication around the project. Please join us there to stay updated about all things tobac that go beyond the detailed work on the code.
diff --git a/doc/big_datasets.rst b/doc/big_datasets.rst
index bac60e83..26851034 100644
--- a/doc/big_datasets.rst
+++ b/doc/big_datasets.rst
@@ -1,11 +1,55 @@
 Handling Large Datasets
 -------------------------------------
 
-Often, one desires to use *tobac* to identify and track features in large datasets ("big data"). This documentation strives to suggest various methods for doing so efficiently. Current versions of *tobac* do not allow for out-of-memory computation, meaning that these strategies may need to be employed for both computational and memory reasons. 
+Often, one desires to use *tobac* to identify and track features in large datasets ("big data"). This documentation strives to suggest various methods for doing so efficiently. Current versions of *tobac* do not support out-of-core (e.g., :code:`dask`) computation, meaning that these strategies may need to be employed for both computational and memory reasons.
 
 .. _Split Feature Detection:
 
 =======================
-Split Feature Detection
+Split Feature Detection and Run in Parallel
 =======================
-Current versions of threshold feature detection (see :doc:`feature_detection_overview`) are time independent, meaning that one can parallelize feature detection across all times (although not across space). *tobac* provides the :py:meth:`tobac.utils.combine_tobac_feats` function to combine a list of dataframes produced by a parallelization method (such as :code:`jug` or :code:`multiprocessing.pool`) into a single combined dataframe suitable to perform tracking with. 
+Current versions of threshold feature detection (see :doc:`feature_detection_overview`) are time independent, meaning that one can easily parallelize feature detection across all times (although not across space). *tobac* provides the :py:meth:`tobac.utils.combine_feature_dataframes` function to combine a list of dataframes produced by a parallelization method (such as :code:`jug`,  :code:`multiprocessing.pool`, or :code:`dask.bag`) into a single combined dataframe suitable to perform tracking with.
+
+Below is a snippet from a larger notebook demonstrating how to run feature detection in parallel ( :doc:`big_datasets_examples/notebooks/parallel_processing_tobac`):
+
+::
+
+    # build list of tracked variables using Dask.Bag
+
+    b = db.from_sequence(
+        [
+            combined_ds["data"][x : x + 1]
+            for x in range(len(combined_ds["time"]))
+        ],
+        npartitions=1,
+    )
+    out_feature_dfs = db.map(
+        lambda x: tobac.feature_detection_multithreshold(
+            x.to_iris(), 4000, **parameters_features
+        ),
+        b,
+    ).compute()
+
+    combined_dataframes = tobac.utils.general.combine_feature_dataframes(out_feature_dfs)
+
+
+.. _Split Segmentation:
+
+======================================
+Split Segmentation and Run in Parallel
+======================================
+Recall that the segmentation mask (see :doc:`segmentation_output`) is the same size as the input grid, which results in large files when handling large input datasets. The following strategies can help reduce the output size and make segmentation masks more useful for the analysis.
+
+The first strategy is to only segment on features *after tracking and quality control*. While this will not directly impact performance, waiting to run segmentation on the final set of features (after discarding, e.g., non-tracked cells) can make analysis of the output segmentation dataset easier.
+
+To enhance the speed at which segmentation runs, one can process multiple segmentation times in parallel independently, similar to feature detection. Unlike feature detection, however, there is currently no built-in *tobac* method to combine multiple segmentation times into a single file. While one can do this using typical NetCDF tools such as :code:`nccat` or with xarray utilities such as :code:`xr.concat`, you can also leave the segmentation mask output as separate files, opening them later with multiple file retrievals such as :code:`xr.open_mfdataset`.
+
+
+.. _Tracking Hanging:
+
+=====================================
+Tracking Hangs with too many Features
+=====================================
+
+When tracking on a large dataset, :code:`tobac.tracking.linking_trackpy` can hang using the default parameters. This is due to the tracking library :code:`trackpy` searching for the next timestep's feature in too large of an area. This can be solved *without impact to scientific output* by lowering the :code:`subnetwork_size` parameter in :code:`tobac.tracking.linking_trackpy`.
+
diff --git a/doc/big_datasets_examples/notebooks/parallel_processing_tobac.ipynb b/doc/big_datasets_examples/notebooks/parallel_processing_tobac.ipynb
new file mode 100644
index 00000000..c65f698a
--- /dev/null
+++ b/doc/big_datasets_examples/notebooks/parallel_processing_tobac.ipynb
@@ -0,0 +1,1160 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Parallel Processing of Feature Detection with `dask`"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook demonstrates how to run *tobac* feature detection in parallel using the `dask` library as the parallel processor."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Imports and Dask Cluster Setup"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/seanfreeman/mambaforge/envs/tobac_installed/lib/python3.11/site-packages/tobac/utils/decorators.py:366: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n",
+      "/Users/seanfreeman/mambaforge/envs/tobac_installed/lib/python3.11/site-packages/tobac/utils/decorators.py:366: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n"
+     ]
+    }
+   ],
+   "source": [
+    "import tobac\n",
+    "import dask.bag as db\n",
+    "import xarray as xr\n",
+    "import s3fs"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are many different ways to initialize a dask cluster. This is just one example, running two workers on a single local machine. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\"> </div>\n",
+       "    <div style=\"margin-left: 48px;\">\n",
+       "        <h3 style=\"margin-bottom: 0px;\">Client</h3>\n",
+       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Client-5a603f42-cd1a-11ee-a16f-06803d36ca4b</p>\n",
+       "        <table style=\"width: 100%; text-align: left;\">\n",
+       "\n",
+       "        <tr>\n",
+       "        \n",
+       "            <td style=\"text-align: left;\"><strong>Connection method:</strong> Cluster object</td>\n",
+       "            <td style=\"text-align: left;\"><strong>Cluster type:</strong> distributed.LocalCluster</td>\n",
+       "        \n",
+       "        </tr>\n",
+       "\n",
+       "        \n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\"></td>\n",
+       "            </tr>\n",
+       "        \n",
+       "\n",
+       "        </table>\n",
+       "\n",
+       "        \n",
+       "\n",
+       "        \n",
+       "            <details>\n",
+       "            <summary style=\"margin-bottom: 20px;\"><h3 style=\"display: inline;\">Cluster Info</h3></summary>\n",
+       "            <div class=\"jp-RenderedHTMLCommon jp-RenderedHTML jp-mod-trusted jp-OutputArea-output\">\n",
+       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\">\n",
+       "    </div>\n",
+       "    <div style=\"margin-left: 48px;\">\n",
+       "        <h3 style=\"margin-bottom: 0px; margin-top: 0px;\">LocalCluster</h3>\n",
+       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">e2f9b06d</p>\n",
+       "        <table style=\"width: 100%; text-align: left;\">\n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Workers:</strong> 2\n",
+       "                </td>\n",
+       "            </tr>\n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Total threads:</strong> 2\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Total memory:</strong> 64.00 GiB\n",
+       "                </td>\n",
+       "            </tr>\n",
+       "            \n",
+       "            <tr>\n",
+       "    <td style=\"text-align: left;\"><strong>Status:</strong> running</td>\n",
+       "    <td style=\"text-align: left;\"><strong>Using processes:</strong> True</td>\n",
+       "</tr>\n",
+       "\n",
+       "            \n",
+       "        </table>\n",
+       "\n",
+       "        <details>\n",
+       "            <summary style=\"margin-bottom: 20px;\">\n",
+       "                <h3 style=\"display: inline;\">Scheduler Info</h3>\n",
+       "            </summary>\n",
+       "\n",
+       "            <div style=\"\">\n",
+       "    <div>\n",
+       "        <div style=\"width: 24px; height: 24px; background-color: #FFF7E5; border: 3px solid #FF6132; border-radius: 5px; position: absolute;\"> </div>\n",
+       "        <div style=\"margin-left: 48px;\">\n",
+       "            <h3 style=\"margin-bottom: 0px;\">Scheduler</h3>\n",
+       "            <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Scheduler-b2b92b52-530a-4c41-b0b5-5da2c4cfcafb</p>\n",
+       "            <table style=\"width: 100%; text-align: left;\">\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Comm:</strong> tcp://127.0.0.1:51973\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Workers:</strong> 2\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Total threads:</strong> 2\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Started:</strong> Just now\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Total memory:</strong> 64.00 GiB\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "            </table>\n",
+       "        </div>\n",
+       "    </div>\n",
+       "\n",
+       "    <details style=\"margin-left: 48px;\">\n",
+       "        <summary style=\"margin-bottom: 20px;\">\n",
+       "            <h3 style=\"display: inline;\">Workers</h3>\n",
+       "        </summary>\n",
+       "\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 0</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:51981\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 1\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:51982/status\" target=\"_blank\">http://127.0.0.1:51982/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 32.00 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:51976\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /var/folders/bj/m6g82c6n41g83y3_dx02y7ch0000gp/T/dask-scratch-space/worker-pou0krcp\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 1</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:51980\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 1\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:51983/status\" target=\"_blank\">http://127.0.0.1:51983/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 32.00 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:51978\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /var/folders/bj/m6g82c6n41g83y3_dx02y7ch0000gp/T/dask-scratch-space/worker-cygts80d\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "\n",
+       "    </details>\n",
+       "</div>\n",
+       "\n",
+       "        </details>\n",
+       "    </div>\n",
+       "</div>\n",
+       "            </details>\n",
+       "        \n",
+       "\n",
+       "    </div>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "<Client: 'tcp://127.0.0.1:51973' processes=2 threads=2, memory=64.00 GiB>"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from dask.distributed import Client, progress\n",
+    "\n",
+    "client = Client(n_workers=2, threads_per_worker=1)\n",
+    "client"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Read in Data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here, we are using the NOAA Global Mosaic of Geostationary Satellite Imagery (GMGSI) as our input data source from AWS s3."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs = s3fs.S3FileSystem(anon=True)\n",
+    "aws_urls = [\n",
+    "    \"s3://noaa-gmgsi-pds/GMGSI_LW/2024/01/01/00/GLOBCOMPLIR_nc.2024010100\",\n",
+    "    \"s3://noaa-gmgsi-pds/GMGSI_LW/2024/01/01/01/GLOBCOMPLIR_nc.2024010101\",\n",
+    "]\n",
+    "\n",
+    "all_ds = list()\n",
+    "for aws_url in aws_urls:\n",
+    "    fileObj = fs.open(aws_url)\n",
+    "    all_ds.append(xr.open_dataset(fileObj, engine=\"h5netcdf\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We loaded in two files and we will use xarray to concatenate them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "combined_ds = xr.concat(all_ds, dim=\"time\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
+       "<defs>\n",
+       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
+       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
+       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
+       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
+       "</symbol>\n",
+       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
+       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
+       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
+       "</symbol>\n",
+       "</defs>\n",
+       "</svg>\n",
+       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
+       " *\n",
+       " */\n",
+       "\n",
+       ":root {\n",
+       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
+       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
+       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
+       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
+       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
+       "  --xr-background-color: var(--jp-layout-color0, white);\n",
+       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
+       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
+       "}\n",
+       "\n",
+       "html[theme=dark],\n",
+       "body[data-theme=dark],\n",
+       "body.vscode-dark {\n",
+       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
+       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
+       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
+       "  --xr-border-color: #1F1F1F;\n",
+       "  --xr-disabled-color: #515151;\n",
+       "  --xr-background-color: #111111;\n",
+       "  --xr-background-color-row-even: #111111;\n",
+       "  --xr-background-color-row-odd: #313131;\n",
+       "}\n",
+       "\n",
+       ".xr-wrap {\n",
+       "  display: block !important;\n",
+       "  min-width: 300px;\n",
+       "  max-width: 700px;\n",
+       "}\n",
+       "\n",
+       ".xr-text-repr-fallback {\n",
+       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-header {\n",
+       "  padding-top: 6px;\n",
+       "  padding-bottom: 6px;\n",
+       "  margin-bottom: 4px;\n",
+       "  border-bottom: solid 1px var(--xr-border-color);\n",
+       "}\n",
+       "\n",
+       ".xr-header > div,\n",
+       ".xr-header > ul {\n",
+       "  display: inline;\n",
+       "  margin-top: 0;\n",
+       "  margin-bottom: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-obj-type,\n",
+       ".xr-array-name {\n",
+       "  margin-left: 2px;\n",
+       "  margin-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-obj-type {\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-sections {\n",
+       "  padding-left: 0 !important;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input {\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input + label {\n",
+       "  color: var(--xr-disabled-color);\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input:enabled + label {\n",
+       "  cursor: pointer;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-section-item input:enabled + label:hover {\n",
+       "  color: var(--xr-font-color0);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary {\n",
+       "  grid-column: 1;\n",
+       "  color: var(--xr-font-color2);\n",
+       "  font-weight: 500;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary > span {\n",
+       "  display: inline-block;\n",
+       "  padding-left: 0.5em;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:disabled + label {\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in + label:before {\n",
+       "  display: inline-block;\n",
+       "  content: '►';\n",
+       "  font-size: 11px;\n",
+       "  width: 15px;\n",
+       "  text-align: center;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:disabled + label:before {\n",
+       "  color: var(--xr-disabled-color);\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked + label:before {\n",
+       "  content: '▼';\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked + label > span {\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary,\n",
+       ".xr-section-inline-details {\n",
+       "  padding-top: 4px;\n",
+       "  padding-bottom: 4px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-inline-details {\n",
+       "  grid-column: 2 / -1;\n",
+       "}\n",
+       "\n",
+       ".xr-section-details {\n",
+       "  display: none;\n",
+       "  grid-column: 1 / -1;\n",
+       "  margin-bottom: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-array-wrap {\n",
+       "  grid-column: 1 / -1;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 20px auto;\n",
+       "}\n",
+       "\n",
+       ".xr-array-wrap > label {\n",
+       "  grid-column: 1;\n",
+       "  vertical-align: top;\n",
+       "}\n",
+       "\n",
+       ".xr-preview {\n",
+       "  color: var(--xr-font-color3);\n",
+       "}\n",
+       "\n",
+       ".xr-array-preview,\n",
+       ".xr-array-data {\n",
+       "  padding: 0 5px !important;\n",
+       "  grid-column: 2;\n",
+       "}\n",
+       "\n",
+       ".xr-array-data,\n",
+       ".xr-array-in:checked ~ .xr-array-preview {\n",
+       "  display: none;\n",
+       "}\n",
+       "\n",
+       ".xr-array-in:checked ~ .xr-array-data,\n",
+       ".xr-array-preview {\n",
+       "  display: inline-block;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list {\n",
+       "  display: inline-block !important;\n",
+       "  list-style: none;\n",
+       "  padding: 0 !important;\n",
+       "  margin: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list li {\n",
+       "  display: inline-block;\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list:before {\n",
+       "  content: '(';\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list:after {\n",
+       "  content: ')';\n",
+       "}\n",
+       "\n",
+       ".xr-dim-list li:not(:last-child):after {\n",
+       "  content: ',';\n",
+       "  padding-right: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-has-index {\n",
+       "  font-weight: bold;\n",
+       "}\n",
+       "\n",
+       ".xr-var-list,\n",
+       ".xr-var-item {\n",
+       "  display: contents;\n",
+       "}\n",
+       "\n",
+       ".xr-var-item > div,\n",
+       ".xr-var-item label,\n",
+       ".xr-var-item > .xr-var-name span {\n",
+       "  background-color: var(--xr-background-color-row-even);\n",
+       "  margin-bottom: 0;\n",
+       "}\n",
+       "\n",
+       ".xr-var-item > .xr-var-name:hover span {\n",
+       "  padding-right: 5px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-list > li:nth-child(odd) > div,\n",
+       ".xr-var-list > li:nth-child(odd) > label,\n",
+       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
+       "  background-color: var(--xr-background-color-row-odd);\n",
+       "}\n",
+       "\n",
+       ".xr-var-name {\n",
+       "  grid-column: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-var-dims {\n",
+       "  grid-column: 2;\n",
+       "}\n",
+       "\n",
+       ".xr-var-dtype {\n",
+       "  grid-column: 3;\n",
+       "  text-align: right;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-var-preview {\n",
+       "  grid-column: 4;\n",
+       "}\n",
+       "\n",
+       ".xr-index-preview {\n",
+       "  grid-column: 2 / 5;\n",
+       "  color: var(--xr-font-color2);\n",
+       "}\n",
+       "\n",
+       ".xr-var-name,\n",
+       ".xr-var-dims,\n",
+       ".xr-var-dtype,\n",
+       ".xr-preview,\n",
+       ".xr-attrs dt {\n",
+       "  white-space: nowrap;\n",
+       "  overflow: hidden;\n",
+       "  text-overflow: ellipsis;\n",
+       "  padding-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-var-name:hover,\n",
+       ".xr-var-dims:hover,\n",
+       ".xr-var-dtype:hover,\n",
+       ".xr-attrs dt:hover {\n",
+       "  overflow: visible;\n",
+       "  width: auto;\n",
+       "  z-index: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs,\n",
+       ".xr-var-data,\n",
+       ".xr-index-data {\n",
+       "  display: none;\n",
+       "  background-color: var(--xr-background-color) !important;\n",
+       "  padding-bottom: 5px !important;\n",
+       "}\n",
+       "\n",
+       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
+       ".xr-var-data-in:checked ~ .xr-var-data,\n",
+       ".xr-index-data-in:checked ~ .xr-index-data {\n",
+       "  display: block;\n",
+       "}\n",
+       "\n",
+       ".xr-var-data > table {\n",
+       "  float: right;\n",
+       "}\n",
+       "\n",
+       ".xr-var-name span,\n",
+       ".xr-var-data,\n",
+       ".xr-index-name div,\n",
+       ".xr-index-data,\n",
+       ".xr-attrs {\n",
+       "  padding-left: 25px !important;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs,\n",
+       ".xr-var-attrs,\n",
+       ".xr-var-data,\n",
+       ".xr-index-data {\n",
+       "  grid-column: 1 / -1;\n",
+       "}\n",
+       "\n",
+       "dl.xr-attrs {\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "  display: grid;\n",
+       "  grid-template-columns: 125px auto;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt,\n",
+       ".xr-attrs dd {\n",
+       "  padding: 0;\n",
+       "  margin: 0;\n",
+       "  float: left;\n",
+       "  padding-right: 10px;\n",
+       "  width: auto;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt {\n",
+       "  font-weight: normal;\n",
+       "  grid-column: 1;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dt:hover span {\n",
+       "  display: inline-block;\n",
+       "  background: var(--xr-background-color);\n",
+       "  padding-right: 10px;\n",
+       "}\n",
+       "\n",
+       ".xr-attrs dd {\n",
+       "  grid-column: 2;\n",
+       "  white-space: pre-wrap;\n",
+       "  word-break: break-all;\n",
+       "}\n",
+       "\n",
+       ".xr-icon-database,\n",
+       ".xr-icon-file-text2,\n",
+       ".xr-no-icon {\n",
+       "  display: inline-block;\n",
+       "  vertical-align: middle;\n",
+       "  width: 1em;\n",
+       "  height: 1.5em !important;\n",
+       "  stroke-width: 0;\n",
+       "  stroke: currentColor;\n",
+       "  fill: currentColor;\n",
+       "}\n",
+       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
+       "Dimensions:  (time: 2, yc: 3000, xc: 4999)\n",
+       "Coordinates:\n",
+       "    lat      (yc, xc) float32 72.72 72.72 72.72 72.72 ... -72.74 -72.74 -72.74\n",
+       "    lon      (yc, xc) float32 180.0 -179.9 -179.9 -179.8 ... 179.8 179.8 179.9\n",
+       "  * time     (time) datetime64[ns] 2024-01-01 2024-01-01T01:00:00\n",
+       "Dimensions without coordinates: yc, xc\n",
+       "Data variables:\n",
+       "    data     (time, yc, xc) float32 206.0 204.0 204.0 ... 187.0 188.0 182.0\n",
+       "Attributes:\n",
+       "    Conventions:          CF-1.4\n",
+       "    Source:               McIDAS Area File\n",
+       "    Satellite Sensor:     DERIVED DATA\n",
+       "    time_coverage_start:  2024-01-01T00:00:00\n",
+       "    instrument_name:      GLOBCOMPLIR\n",
+       "    history:              Mon Jan  1 00:38:21 2024: ncks -d xc,0,4998 templir...\n",
+       "    NCO:                  netCDF Operators version 4.7.5 (Homepage = http://n...</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-ca79bc41-f441-44f3-b307-ccab45081cec' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-ca79bc41-f441-44f3-b307-ccab45081cec' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>time</span>: 2</li><li><span>yc</span>: 3000</li><li><span>xc</span>: 4999</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-d2ed5bb4-a962-40a1-a752-9bb5c4c4e5a8' class='xr-section-summary-in' type='checkbox'  checked><label for='section-d2ed5bb4-a962-40a1-a752-9bb5c4c4e5a8' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>lat</span></div><div class='xr-var-dims'>(yc, xc)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>72.72 72.72 72.72 ... -72.74 -72.74</div><input id='attrs-972988bb-a0b6-4621-9ccf-c947a4ae0ee6' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-972988bb-a0b6-4621-9ccf-c947a4ae0ee6' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-9f244669-15f1-437a-928b-e89b36e47ae1' class='xr-var-data-in' type='checkbox'><label for='data-9f244669-15f1-437a-928b-e89b36e47ae1' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>lat</dd><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>valid_range :</span></dt><dd>[-90.  90.]</dd><dt><span>standard_name :</span></dt><dd>latitude</dd></dl></div><div class='xr-var-data'><pre>array([[ 72.71541 ,  72.71541 ,  72.71541 , ...,  72.71541 ,  72.71541 ,\n",
+       "         72.71541 ],\n",
+       "       [ 72.694   ,  72.694   ,  72.694   , ...,  72.694   ,  72.694   ,\n",
+       "         72.694   ],\n",
+       "       [ 72.67256 ,  72.67256 ,  72.67256 , ...,  72.67256 ,  72.67256 ,\n",
+       "         72.67256 ],\n",
+       "       ...,\n",
+       "       [-72.69398 , -72.69398 , -72.69398 , ..., -72.69398 , -72.69398 ,\n",
+       "        -72.69398 ],\n",
+       "       [-72.715385, -72.715385, -72.715385, ..., -72.715385, -72.715385,\n",
+       "        -72.715385],\n",
+       "       [-72.73677 , -72.73677 , -72.73677 , ..., -72.73677 , -72.73677 ,\n",
+       "        -72.73677 ]], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>lon</span></div><div class='xr-var-dims'>(yc, xc)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>180.0 -179.9 -179.9 ... 179.8 179.9</div><input id='attrs-daa3e7a2-a40a-4cac-b03d-b4dbb545cfa9' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-daa3e7a2-a40a-4cac-b03d-b4dbb545cfa9' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-0af61e05-cdf1-4d4f-b277-fcf829b11a2d' class='xr-var-data-in' type='checkbox'><label for='data-0af61e05-cdf1-4d4f-b277-fcf829b11a2d' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>lon</dd><dt><span>units :</span></dt><dd>degrees_east</dd><dt><span>valid_range :</span></dt><dd>[-180.  180.]</dd><dt><span>standard_name :</span></dt><dd>longitude</dd></dl></div><div class='xr-var-data'><pre>array([[ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038],\n",
+       "       [ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038],\n",
+       "       [ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038],\n",
+       "       ...,\n",
+       "       [ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038],\n",
+       "       [ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038],\n",
+       "       [ 179.99962, -179.92838, -179.85637, ...,  179.75636,  179.82837,\n",
+       "         179.90038]], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>2024-01-01 2024-01-01T01:00:00</div><input id='attrs-58c34275-6565-4afe-ac44-2070de1cab69' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-58c34275-6565-4afe-ac44-2070de1cab69' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-58a7264f-f87c-4736-a7f7-32bb253df5ba' class='xr-var-data-in' type='checkbox'><label for='data-58a7264f-f87c-4736-a7f7-32bb253df5ba' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>seconds since 1970-1-1 0:0:0</dd></dl></div><div class='xr-var-data'><pre>array([&#x27;2024-01-01T00:00:00.000000000&#x27;, &#x27;2024-01-01T01:00:00.000000000&#x27;],\n",
+       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-a394f17e-c062-4169-b7f6-55d79349f9db' class='xr-section-summary-in' type='checkbox'  checked><label for='section-a394f17e-c062-4169-b7f6-55d79349f9db' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>data</span></div><div class='xr-var-dims'>(time, yc, xc)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>206.0 204.0 204.0 ... 188.0 182.0</div><input id='attrs-b7ae04f8-f563-48c2-964e-2ec7d905e58a' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-b7ae04f8-f563-48c2-964e-2ec7d905e58a' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-271f86ab-32ee-4f33-a34b-9b2ec1055dcb' class='xr-var-data-in' type='checkbox'><label for='data-271f86ab-32ee-4f33-a34b-9b2ec1055dcb' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>0-255 Brightness Temperature</dd><dt><span>type :</span></dt><dd>VISR</dd><dt><span>units :</span></dt><dd>none</dd></dl></div><div class='xr-var-data'><pre>array([[[206., 204., 204., ..., 207., 206., 206.],\n",
+       "        [206., 204., 204., ..., 207., 206., 206.],\n",
+       "        [204., 202., 204., ..., 206., 204., 204.],\n",
+       "        ...,\n",
+       "        [185., 186., 186., ..., 191., 190., 185.],\n",
+       "        [185., 186., 186., ..., 191., 190., 185.],\n",
+       "        [185., 186., 186., ..., 191., 190., 185.]],\n",
+       "\n",
+       "       [[205., 204., 204., ..., 205., 205., 205.],\n",
+       "        [205., 204., 204., ..., 205., 205., 205.],\n",
+       "        [204., 203., 203., ..., 205., 204., 204.],\n",
+       "        ...,\n",
+       "        [182., 183., 183., ..., 187., 188., 182.],\n",
+       "        [182., 183., 183., ..., 187., 188., 182.],\n",
+       "        [182., 183., 183., ..., 187., 188., 182.]]], dtype=float32)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-505e07bb-7e75-44d0-98ad-c1c9b85dc854' class='xr-section-summary-in' type='checkbox'  ><label for='section-505e07bb-7e75-44d0-98ad-c1c9b85dc854' class='xr-section-summary' >Indexes: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>time</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-4430f4a2-6202-4892-b45e-cf62d03469f2' class='xr-index-data-in' type='checkbox'/><label for='index-4430f4a2-6202-4892-b45e-cf62d03469f2' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(DatetimeIndex([&#x27;2024-01-01 00:00:00&#x27;, &#x27;2024-01-01 01:00:00&#x27;], dtype=&#x27;datetime64[ns]&#x27;, name=&#x27;time&#x27;, freq=None))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-0d052c3f-dce4-4c1f-ace7-668045bf7c55' class='xr-section-summary-in' type='checkbox'  checked><label for='section-0d052c3f-dce4-4c1f-ace7-668045bf7c55' class='xr-section-summary' >Attributes: <span>(7)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>Conventions :</span></dt><dd>CF-1.4</dd><dt><span>Source :</span></dt><dd>McIDAS Area File</dd><dt><span>Satellite Sensor :</span></dt><dd>DERIVED DATA</dd><dt><span>time_coverage_start :</span></dt><dd>2024-01-01T00:00:00</dd><dt><span>instrument_name :</span></dt><dd>GLOBCOMPLIR</dd><dt><span>history :</span></dt><dd>Mon Jan  1 00:38:21 2024: ncks -d xc,0,4998 templir.nc outlir.nc</dd><dt><span>NCO :</span></dt><dd>netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)</dd></dl></div></li></ul></div></div>"
+      ],
+      "text/plain": [
+       "<xarray.Dataset>\n",
+       "Dimensions:  (time: 2, yc: 3000, xc: 4999)\n",
+       "Coordinates:\n",
+       "    lat      (yc, xc) float32 72.72 72.72 72.72 72.72 ... -72.74 -72.74 -72.74\n",
+       "    lon      (yc, xc) float32 180.0 -179.9 -179.9 -179.8 ... 179.8 179.8 179.9\n",
+       "  * time     (time) datetime64[ns] 2024-01-01 2024-01-01T01:00:00\n",
+       "Dimensions without coordinates: yc, xc\n",
+       "Data variables:\n",
+       "    data     (time, yc, xc) float32 206.0 204.0 204.0 ... 187.0 188.0 182.0\n",
+       "Attributes:\n",
+       "    Conventions:          CF-1.4\n",
+       "    Source:               McIDAS Area File\n",
+       "    Satellite Sensor:     DERIVED DATA\n",
+       "    time_coverage_start:  2024-01-01T00:00:00\n",
+       "    instrument_name:      GLOBCOMPLIR\n",
+       "    history:              Mon Jan  1 00:38:21 2024: ncks -d xc,0,4998 templir...\n",
+       "    NCO:                  netCDF Operators version 4.7.5 (Homepage = http://n..."
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "combined_ds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "These feature detection parameters are just examples."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## *tobac* Feature Detection"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "parameters_features = {}\n",
+    "parameters_features[\"position_threshold\"] = \"weighted_diff\"\n",
+    "parameters_features[\"sigma_threshold\"] = 0.5\n",
+    "parameters_features[\"n_min_threshold\"] = 4\n",
+    "parameters_features[\"target\"] = \"minimum\"\n",
+    "parameters_features[\"threshold\"] = [180, 170]\n",
+    "parameters_features[\"PBC_flag\"] = \"hdim_2\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "While future versions (1.6 and greater) of *tobac* will support xarray natively in feature detection and segmentation, current versions of *tobac* rely on Iris for gridded data. Because of this, we have to make some conversions to have this data be compatible with iris. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# iris issues\n",
+    "combined_ds[\"data\"].attrs[\"units\"] = \"kelvin\"\n",
+    "combined_ds[\"data\"][\"time\"].attrs[\"long_name\"] = \"time\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, we will use a *dask bag* to parallelize our feature detection over time. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/seanfreeman/mambaforge/envs/tobac_installed/lib/python3.11/site-packages/distributed/client.py:3162: UserWarning: Sending large graph of size 57.21 MiB.\n",
+      "This may cause some slowdown.\n",
+      "Consider scattering data ahead of time and using futures.\n",
+      "  warnings.warn(\n",
+      "/Users/seanfreeman/mambaforge/envs/tobac_installed/lib/python3.11/site-packages/tobac/utils/decorators.py:366: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n",
+      "/Users/seanfreeman/mambaforge/envs/tobac_installed/lib/python3.11/site-packages/tobac/utils/decorators.py:366: UserWarning: Numba not able to be imported; periodic boundary calculations will be slower.Exception raised: ModuleNotFoundError(\"No module named 'numba'\")\n",
+      "  warnings.warn(\n"
+     ]
+    }
+   ],
+   "source": [
+    "b = db.from_sequence(\n",
+    "    [\n",
+    "        combined_ds[\"data\"][x : x + 1][0:500, 0:500]\n",
+    "        for x in range(len(combined_ds[\"time\"]))\n",
+    "    ],\n",
+    "    npartitions=1,\n",
+    ")\n",
+    "out_feature_dfs = db.map(\n",
+    "    lambda x: tobac.feature_detection_multithreshold(\n",
+    "        x.to_iris(), 4000, **parameters_features\n",
+    "    ),\n",
+    "    b,\n",
+    ").compute()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Combining parallel-detected features into one coherent DataFrame"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>frame</th>\n",
+       "      <th>idx</th>\n",
+       "      <th>hdim_1</th>\n",
+       "      <th>hdim_2</th>\n",
+       "      <th>num</th>\n",
+       "      <th>threshold_value</th>\n",
+       "      <th>feature</th>\n",
+       "      <th>time</th>\n",
+       "      <th>timestr</th>\n",
+       "      <th>latitude</th>\n",
+       "      <th>longitude</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>0</td>\n",
+       "      <td>1</td>\n",
+       "      <td>0.975325</td>\n",
+       "      <td>60.718132</td>\n",
+       "      <td>26</td>\n",
+       "      <td>180</td>\n",
+       "      <td>1</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>72.694528</td>\n",
+       "      <td>-175.628134</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>0</td>\n",
+       "      <td>2</td>\n",
+       "      <td>0.670018</td>\n",
+       "      <td>79.074818</td>\n",
+       "      <td>16</td>\n",
+       "      <td>180</td>\n",
+       "      <td>2</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>72.701065</td>\n",
+       "      <td>-174.306289</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>0</td>\n",
+       "      <td>5</td>\n",
+       "      <td>0.358115</td>\n",
+       "      <td>1492.979778</td>\n",
+       "      <td>15</td>\n",
+       "      <td>180</td>\n",
+       "      <td>3</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>72.707742</td>\n",
+       "      <td>-72.492469</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>0</td>\n",
+       "      <td>6</td>\n",
+       "      <td>0.482579</td>\n",
+       "      <td>1531.520215</td>\n",
+       "      <td>26</td>\n",
+       "      <td>180</td>\n",
+       "      <td>4</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>72.705077</td>\n",
+       "      <td>-69.717214</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>0</td>\n",
+       "      <td>7</td>\n",
+       "      <td>3.409896</td>\n",
+       "      <td>2113.185770</td>\n",
+       "      <td>285</td>\n",
+       "      <td>180</td>\n",
+       "      <td>5</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>2024-01-01 00:00:00</td>\n",
+       "      <td>72.642292</td>\n",
+       "      <td>-27.832080</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1565</th>\n",
+       "      <td>1</td>\n",
+       "      <td>1283</td>\n",
+       "      <td>498.180906</td>\n",
+       "      <td>30.005031</td>\n",
+       "      <td>16</td>\n",
+       "      <td>170</td>\n",
+       "      <td>1566</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>58.262049</td>\n",
+       "      <td>-177.839756</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1566</th>\n",
+       "      <td>1</td>\n",
+       "      <td>1284</td>\n",
+       "      <td>497.817803</td>\n",
+       "      <td>715.879216</td>\n",
+       "      <td>7</td>\n",
+       "      <td>170</td>\n",
+       "      <td>1567</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>58.275803</td>\n",
+       "      <td>-128.450668</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1567</th>\n",
+       "      <td>1</td>\n",
+       "      <td>1285</td>\n",
+       "      <td>498.158232</td>\n",
+       "      <td>3295.186171</td>\n",
+       "      <td>63</td>\n",
+       "      <td>170</td>\n",
+       "      <td>1568</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>58.262908</td>\n",
+       "      <td>57.282537</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1568</th>\n",
+       "      <td>1</td>\n",
+       "      <td>1288</td>\n",
+       "      <td>498.454404</td>\n",
+       "      <td>3793.095530</td>\n",
+       "      <td>5</td>\n",
+       "      <td>170</td>\n",
+       "      <td>1569</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>58.251682</td>\n",
+       "      <td>93.136465</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1569</th>\n",
+       "      <td>1</td>\n",
+       "      <td>1293</td>\n",
+       "      <td>498.648911</td>\n",
+       "      <td>3317.641704</td>\n",
+       "      <td>6</td>\n",
+       "      <td>170</td>\n",
+       "      <td>1570</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>2024-01-01 01:00:00</td>\n",
+       "      <td>58.244310</td>\n",
+       "      <td>58.899535</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>1570 rows × 11 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "      frame   idx      hdim_1       hdim_2  num  threshold_value  feature  \\\n",
+       "0         0     1    0.975325    60.718132   26              180        1   \n",
+       "1         0     2    0.670018    79.074818   16              180        2   \n",
+       "2         0     5    0.358115  1492.979778   15              180        3   \n",
+       "3         0     6    0.482579  1531.520215   26              180        4   \n",
+       "4         0     7    3.409896  2113.185770  285              180        5   \n",
+       "...     ...   ...         ...          ...  ...              ...      ...   \n",
+       "1565      1  1283  498.180906    30.005031   16              170     1566   \n",
+       "1566      1  1284  497.817803   715.879216    7              170     1567   \n",
+       "1567      1  1285  498.158232  3295.186171   63              170     1568   \n",
+       "1568      1  1288  498.454404  3793.095530    5              170     1569   \n",
+       "1569      1  1293  498.648911  3317.641704    6              170     1570   \n",
+       "\n",
+       "                     time              timestr   latitude   longitude  \n",
+       "0     2024-01-01 00:00:00  2024-01-01 00:00:00  72.694528 -175.628134  \n",
+       "1     2024-01-01 00:00:00  2024-01-01 00:00:00  72.701065 -174.306289  \n",
+       "2     2024-01-01 00:00:00  2024-01-01 00:00:00  72.707742  -72.492469  \n",
+       "3     2024-01-01 00:00:00  2024-01-01 00:00:00  72.705077  -69.717214  \n",
+       "4     2024-01-01 00:00:00  2024-01-01 00:00:00  72.642292  -27.832080  \n",
+       "...                   ...                  ...        ...         ...  \n",
+       "1565  2024-01-01 01:00:00  2024-01-01 01:00:00  58.262049 -177.839756  \n",
+       "1566  2024-01-01 01:00:00  2024-01-01 01:00:00  58.275803 -128.450668  \n",
+       "1567  2024-01-01 01:00:00  2024-01-01 01:00:00  58.262908   57.282537  \n",
+       "1568  2024-01-01 01:00:00  2024-01-01 01:00:00  58.251682   93.136465  \n",
+       "1569  2024-01-01 01:00:00  2024-01-01 01:00:00  58.244310   58.899535  \n",
+       "\n",
+       "[1570 rows x 11 columns]"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "tobac.utils.general.combine_feature_dataframes(out_feature_dfs)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "tobac_installed",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/bulk_statistics/notebooks/compute_statistics_during_feature_detection.ipynb b/doc/bulk_statistics/notebooks/compute_statistics_during_feature_detection.ipynb
index cb55918f..e9eb2a73 100644
--- a/doc/bulk_statistics/notebooks/compute_statistics_during_feature_detection.ipynb
+++ b/doc/bulk_statistics/notebooks/compute_statistics_during_feature_detection.ipynb
@@ -34,6 +34,7 @@
     "import shutil\n",
     "from six.moves import urllib\n",
     "from pathlib import Path\n",
+    "\n",
     "%matplotlib inline"
    ]
   },
@@ -60,7 +61,8 @@
    "source": [
     "# Import tobac itself\n",
     "import tobac\n",
-    "print('using tobac version', str(tobac.__version__))"
+    "\n",
+    "print(\"using tobac version\", str(tobac.__version__))"
    ]
   },
   {
@@ -78,10 +80,11 @@
    "source": [
     "# Disable a few warnings:\n",
     "import warnings\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)"
+    "\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)"
    ]
   },
   {
@@ -106,7 +109,7 @@
    },
    "outputs": [],
    "source": [
-    "data_out=Path('../')"
+    "data_out = Path(\"../\")"
    ]
   },
   {
@@ -123,18 +126,18 @@
    "outputs": [],
    "source": [
     "# 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\"))"
    ]
   },
   {
@@ -157,7 +160,7 @@
    },
    "outputs": [],
    "source": [
-    "Precip=iris.load_cube(str(data_file[0]),'surface_precipitation_average')"
+    "Precip = iris.load_cube(str(data_file[0]), \"surface_precipitation_average\")"
    ]
   },
   {
@@ -324,7 +327,7 @@
     }
    ],
    "source": [
-    "#display information about the iris cube containing the surface precipitation data:\n",
+    "# display information about the iris cube containing the surface precipitation data:\n",
     "display(Precip)"
    ]
   },
@@ -341,11 +344,11 @@
    },
    "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()"
    ]
@@ -380,17 +383,17 @@
    },
    "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)"
    ]
   },
   {
@@ -421,9 +424,9 @@
    "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"
    ]
   },
   {
@@ -446,7 +449,7 @@
    },
    "outputs": [],
    "source": [
-    "statistics['percentiles'] = (np.percentile, {'q': [95,99]})"
+    "statistics[\"percentiles\"] = (np.percentile, {\"q\": [95, 99]})"
    ]
   },
   {
@@ -485,11 +488,13 @@
    ],
    "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, statistic=statistics) \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(\n",
+    "    Precip, dxy, **parameters_features, statistic=statistics\n",
+    ")\n",
+    "print(\"feature detection done\")\n",
+    "Features.to_hdf(savedir / \"Features.h5\", \"table\")\n",
+    "print(\"features saved\")"
    ]
   },
   {
diff --git a/doc/bulk_statistics/notebooks/compute_statistics_during_segmentation.ipynb b/doc/bulk_statistics/notebooks/compute_statistics_during_segmentation.ipynb
index 21f51511..0ee1b71a 100644
--- a/doc/bulk_statistics/notebooks/compute_statistics_during_segmentation.ipynb
+++ b/doc/bulk_statistics/notebooks/compute_statistics_during_segmentation.ipynb
@@ -32,6 +32,7 @@
     "import shutil\n",
     "from six.moves import urllib\n",
     "from pathlib import Path\n",
+    "\n",
     "%matplotlib inline"
    ]
   },
@@ -50,10 +51,11 @@
    "source": [
     "# Disable a few warnings:\n",
     "import warnings\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)"
+    "\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)"
    ]
   },
   {
@@ -79,7 +81,8 @@
    "source": [
     "# Import tobac itself\n",
     "import tobac\n",
-    "print('using tobac version', str(tobac.__version__))"
+    "\n",
+    "print(\"using tobac version\", str(tobac.__version__))"
    ]
   },
   {
@@ -95,20 +98,20 @@
    },
    "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\"))"
    ]
   },
   {
@@ -124,11 +127,11 @@
    },
    "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()"
    ]
@@ -146,7 +149,7 @@
    },
    "outputs": [],
    "source": [
-    "Precip=iris.load_cube(str(data_file[0]),'surface_precipitation_average')"
+    "Precip = iris.load_cube(str(data_file[0]), \"surface_precipitation_average\")"
    ]
   },
   {
@@ -169,17 +172,17 @@
    },
    "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)"
    ]
   },
   {
@@ -218,11 +221,11 @@
    ],
    "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\")"
    ]
   },
   {
@@ -255,12 +258,12 @@
    "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)"
    ]
   },
   {
@@ -291,9 +294,9 @@
    "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"
    ]
   },
   {
@@ -316,7 +319,7 @@
    },
    "outputs": [],
    "source": [
-    "statistics['percentiles'] = (np.percentile, {'q': [95,99]})"
+    "statistics[\"percentiles\"] = (np.percentile, {\"q\": [95, 99]})"
    ]
   },
   {
@@ -355,12 +358,16 @@
    ],
    "source": [
     "# Perform Segmentation and save resulting mask to NetCDF file:\n",
-    "print('Starting segmentation based on surface precipitation')\n",
-    "Mask,Features_Precip=tobac.segmentation_2D(Features,Precip,dxy,**parameters_segmentation, statistic=statistics)\n",
-    "print('segmentation based on surface precipitation performed, start saving results to files')\n",
-    "iris.save([Mask], 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')"
+    "print(\"Starting segmentation based on surface precipitation\")\n",
+    "Mask, Features_Precip = tobac.segmentation_2D(\n",
+    "    Features, Precip, dxy, **parameters_segmentation, statistic=statistics\n",
+    ")\n",
+    "print(\n",
+    "    \"segmentation based on surface precipitation performed, start saving results to files\"\n",
+    ")\n",
+    "iris.save([Mask], 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\")"
    ]
   },
   {
diff --git a/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb b/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
index 112f07f2..3ed2bb69 100644
--- a/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
+++ b/doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb
@@ -29,12 +29,13 @@
     "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"
    ]
   },
@@ -61,7 +62,8 @@
    "source": [
     "# Import tobac itself\n",
     "import tobac\n",
-    "print('using tobac version', str(tobac.__version__))"
+    "\n",
+    "print(\"using tobac version\", str(tobac.__version__))"
    ]
   },
   {
@@ -79,10 +81,11 @@
    "source": [
     "# Disable a few warnings:\n",
     "import warnings\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)"
+    "\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)"
    ]
   },
   {
@@ -105,11 +108,11 @@
    },
    "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()"
    ]
@@ -127,20 +130,20 @@
    },
    "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\"))"
    ]
   },
   {
@@ -156,7 +159,7 @@
    },
    "outputs": [],
    "source": [
-    "Precip=iris.load_cube(str(data_file[0]),'surface_precipitation_average')"
+    "Precip = iris.load_cube(str(data_file[0]), \"surface_precipitation_average\")"
    ]
   },
   {
@@ -172,17 +175,17 @@
    },
    "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)"
    ]
   },
   {
@@ -221,11 +224,11 @@
    ],
    "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\")"
    ]
   },
   {
@@ -249,12 +252,12 @@
    "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)"
    ]
   },
   {
@@ -293,12 +296,18 @@
    ],
    "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(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')"
+    "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\")"
    ]
   },
   {
@@ -354,9 +363,9 @@
    "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"
    ]
   },
   {
@@ -379,7 +388,7 @@
    },
    "outputs": [],
    "source": [
-    "statistics['percentiles'] = (np.percentile, {'q': [95,99]})"
+    "statistics[\"percentiles\"] = (np.percentile, {\"q\": [95, 99]})"
    ]
   },
   {
@@ -395,7 +404,9 @@
    },
    "outputs": [],
    "source": [
-    "features_with_stats = get_statistics_from_mask(Features_Precip, Mask_Precip, Precip, statistic=statistics)"
+    "features_with_stats = get_statistics_from_mask(\n",
+    "    Features_Precip, Mask_Precip, Precip, statistic=statistics\n",
+    ")"
    ]
   },
   {
diff --git a/doc/code_reviews.rst b/doc/code_reviews.rst
new file mode 100644
index 00000000..da8c2c90
--- /dev/null
+++ b/doc/code_reviews.rst
@@ -0,0 +1,36 @@
+Code reviews
+------------------
+
+Before anything is merged into the release branch (:code:`RC_*`), we require that two reviewers accept the code changes of a pull request. 
+
+============================
+How to do a code review 
+============================
+
+* Checkout out pull request locally (`how to checkout a pull request locally <https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/reviewing-changes-in-pull-requests/checking-out-pull-requests-locally>`_)
+
+* Run tests locally
+  
+* Go through code and see if it is readable and easy to understand
+  
+* Not required, but often useful: test new features with your own data 
+
+
+============================
+Tips and expectations
+============================
+
+
+Doing a code review can be very challenging if you are unfamiliar with the process. Here is a set of documents which might provide a good resource on how to get started:
+
+https://github.com/google/eng-practices
+
+
+=========================
+Conventional comments
+=========================
+
+The comments in a code review should be clear and constructive.
+
+A useful way of highlighting the intention of specific comments is to label them according to `conventional comments <https://conventionalcomments.org/>`_.
+
diff --git a/doc/code_structure.rst b/doc/code_structure.rst
new file mode 100644
index 00000000..cbc1b23e
--- /dev/null
+++ b/doc/code_structure.rst
@@ -0,0 +1,72 @@
+Code structure and key design concepts 
+--------------------------------------
+
+==================================
+Modules
+==================================
+
+**tobac** aims to provide a flexible and modular framework which can be seen as a toolbox to create tracking algorithms according to the user's specific research needs. 
+
+The **tobac** package currently consists of three **main modules**:
+
+1. The :py:mod:`tobac.feature_detection` contains methods to identify objects (*features*) in 2D or 3D (3D or 4D when including the time dimensions) gridded data. This is done by identifying contiguous regions above or below one or multiple user-defined thresholds. The module makes use of :py:mod:`scipy.ndimage.label`, a generic image processing method that labels features in an array. The methods in :py:mod:`tobac.feature_detection` are high-level functions that enable a fast and effective feature detection and create easy-to-use output in form of a :py:mod:`pandas.DataFrame` that contains the coordinates and some basic information on each detected feature. The most high-level methods that is commonly used by users is :py:func:`tobac.feature_detection_multithreshold`. 
+
+2. The :py:mod:`tobac.segmentation` module contains methods to define the extent of the identified feature areas or volumes. This step is needed to create a mask of the identified features because the feature detection currently only saves the center points of the features. The segmentation procedure is performed by using the watershedding method, but more methods are to be implemented in the future. Just as the feature detection, this module can handle both 2D and 3D data. 
+
+3. The :py:mod:`tobac.tracking` module is responsible for linking identified features over time. This module makes primarily use of the python package :py:mod:`trackpy`. Note that the linking using :py:mod:`trackpy` is based on particle tracking principles which means that only the feature center positions (not the entire area or volume associated with each feature) are needed to link features over time. Other methods such as tracking based on overlapping areas from the segmented features are to be implemented.
+
+In addition to the main modules, there are three **postprocessing modules**: 
+
+4. The :py:mod:`tobac.merge_split` module provides functionality to identify mergers and splitters in the tracking output and to add labels such that one can reconstruct the parent and child tracks of each cell. 
+
+5. The :py:mod:`tobac.analysis` module contains methods to analyze the tracking output and derive statistics about individual tracks as well as summary statistics of the entire populations of tracks or subsets of the latter. 
+
+6. The :py:mod:`tobac.plotting` module provides methods to visualize the tracking output, for example for creating maps and animations of identified features, segmented areas and tracks.
+
+   
+Finally, there are two modules that are primarily **important for developers**:
+
+7. The :py:mod:`tobac.utils` module is a collection of smaller, not necessarily tracking-specific methods that facilitate and support the methods of the main modules. This module has multiple submodules. We separate methods that are rather generic and could also be practical for tobac users who build their own tracking algorithms (:py:mod:`tobac.utils.general`) and methods that mainly facilitate the development of **tobac** (:py:mod:`tobac.utils.internal`). Sometimes, new features come with the need of a whole set of new methods, so it could make sense to save these in their own submodule (see e.g. :py:mod:`tobac.periodic_boundaries`)
+
+8. The :py:mod:`tobac.testing` module provides support for writing of unit tests. This module contains several methods to create simplified test data sets on which the various methods and parameters for feature detection, segmentation, and tracking can be tested. 
+
+For more information on each submodule, refer to the respective source code documentation.
+
+One thing to note is that **tobac** as of now is purely functional. The plan is, however, to move towards a more object-oriented design with base classes for the main operations such as feature detection and tracking. 
+
+
+========
+Examples
+========
+
+To help users get started with **tobac** and to demonstrate the various functionalities, **tobac** hosts several detailed and **illustrated examples** in the form of Jupyter notebooks. They are hosted under the directory `examples/` and be executed by the user. Our readthedocs page also hosts a rendered version of our examples as `gallery <https://tobac.readthedocs.io/en/latest/examples.html>`_
+
+
+============================
+Migrating to xarray and dask
+============================
+
+Currently, **tobac** uses `iris cubes <https://scitools-iris.readthedocs.io/en/latest/userguide/iris_cubes.html>`_ as the 
+primary data container. However, we are currently working on migrating the source code to 
+`xarray <https://docs.xarray.dev/en/stable/>`_ such that all internal functions are based on `xr.DataArray 
+objects <https://docs.xarray.dev/en/stable/generated/xarray.DataArray.html>`_. 
+
+To ensure a robust transition from **iris** to **xarray**, we make use of various decorators that convert input and 
+output data for the main functions without changing their actual code. These decorators are located in the `decorator 
+submodule <https://github.com/tobac-project/tobac/blob/main/tobac/utils/decorators.py>`_. 
+
+In addition, one of our main goals for the future is to fully support `dask <https://www.dask.org/>`_, in order to scale 
+to large datasets and enable parallelization.  
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/conf.py b/doc/conf.py
index a0c1b913..02a7c5b6 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -20,6 +20,7 @@
     "sphinx.ext.napoleon",
     "nbsphinx",
     "sphinx_gallery.load_style",
+    "sphinx_toolbox.code",
 ]
 
 
@@ -28,6 +29,7 @@
 html_static_path = ["_static"]
 exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
 
+
 project = "tobac"
 
 master_doc = "index"
@@ -62,6 +64,8 @@ def setup(app):
 
 sys.path.insert(0, os.path.abspath("../"))
 
+pygments_style = "sphinx"
+
 # Napoleon settings for configuring the Napoleon extension
 # See documentation here:
 # https://www.sphinx-doc.org/en/master/usage/extensions/napoleon.html
diff --git a/doc/contributing.rst b/doc/contributing.rst
new file mode 100644
index 00000000..c458a69d
--- /dev/null
+++ b/doc/contributing.rst
@@ -0,0 +1,168 @@
+..
+    How to contribute to the tobac project
+
+How to contribute 
+-------------------------
+
+Step-by-step overview of most important points: https://github.com/tobac-project/tobac/blob/main/CONTRIBUTING.md
+
+=========================
+Code of conduct 
+=========================
+
+We are a multi-institutional and international community that aims to maintain and increase our diversity. We acknowledge that we all come with different experiences and capacities. Therefore, we strive to foster an inclusive and respectful environment where we help and support each other. We welcome any types of contributions and believe that we together can create accessible, reusable, and maintanable code that empowers researchers and enables groundbreaking science. 
+
+We would like to refer to the `Python code of conduct <https://www.python.org/psf/conduct/>`_ as we follow the same principles for communication and working with each other!
+
+=========================
+git basics
+=========================
+
+* **Create a Github account**: The first thing, you need to do is to `create a GitHub account <https://docs.github.com/en/get-started/signing-up-for-github/signing-up-for-a-new-github-account>`_ if you do not already have one. 
+
+* **Get familiar with the basics of GitHub and git**:
+   * Getting started with the `basics <https://docs.github.com/en/get-started/quickstart/hello-world>`_
+   * Learn about `branches <https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/about-branches>`_ 
+   * Learn about `forks <https://docs.github.com/en/get-started/quickstart/fork-a-repo>`_
+   * Learn about `pull requests <https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/about-pull-requests>`_
+   * Learn about `how to commit and push changes from your local repository <https://docs.github.com/en/migrations/importing-source-code/using-the-command-line-to-import-source-code/adding-locally-hosted-code-to-github>`_ 
+
+* **Create an issue**: If you have an idea for a new feature or a suggestion for any kind of code changes, please create an issue for this. We sort `our issues <https://github.com/tobac-project/tobac/issues>`_ into `milestones <https://github.com/tobac-project/tobac/milestones>`_ to priorize work and manage our workflow, i.e. the different versions of **tobac** to come.
+  
+  The issues act, therefore, not only as a place for reporting bugs, but also as a collection of *to do* points. 
+
+* **Work on an issue**: You can also work on any issue that was created by somebody else and is already out there. A tip is to look for the **good first issue** label, if you are a new developer. These issues are usually fairly easy to address and can be good to practice our GitHub workflow. 
+
+  
+* **Create a pull request from your fork:** We use our personal forks of the tobac repository to create pull requests. This means that you have to first commit and push your local changes to your personal fork and then create a pull request from that fork: https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/creating-a-pull-request-from-a-fork
+
+===================================
+Writing proper documentation
+===================================
+
+Please provide **Numpy Docstrings** for all new functions. 
+
+**Example**:
+
+.. code::
+
+   '''
+   Calculates centre of gravity and mass for each individually tracked cell in the simulation.
+
+
+    Parameters
+    ----------
+    tracks : pandas.DataFram
+        DataFrame containing trajectories of cell centres
+        
+    param mass : iris.cube.Cube
+        cube of quantity (need coordinates 'time', 'geopotential_height','projection_x_coordinate' and 
+        'projection_y_coordinate')
+        
+    param mask : iris.cube.Cube
+        cube containing mask (int > where belonging to cloud volume, 0 everywhere else )
+
+
+    Returns
+    -------
+    track_out : pandas.DataFrame
+        Dataframe containing t,x,y,z positions of centre of gravity and total cloud mass each tracked cells 
+        at each timestep
+    
+    '''
+
+
+
+===================================
+Tips on working on your local code
+===================================
+
+* Install tobac package with :code:`pip install -e` 
+    * This allows you to directly test your local code changes as you run tobac. Instead of using the **tobac** version of the latest release, your local version of tobac will be used when you import **tobac** in a python script or IDE. 
+    * *Note that* this way of installing a local package will use the code of the checked in branch, so this allows you also to test code while switching between branches.
+
+* You can locally **build the documentation page**: 
+    * see :doc:`testing_sphinx-based_rendering`
+
+* Writing `meaningful commit messages <https://www.conventionalcommits.org/en/v1.0.0/>`_ can be very helpful for you and people who review your code to better understand the code changes.
+
+
+=========================
+Our branching strategy
+=========================
+
+While you can use any type of branching strategy and naming as you work in your personal fork, we have three branches in the tobac repository: 
+
+* :code:`RC_*`
+* :code:`dev_*`
+* :code:`hotfix`
+
+:code:`RC_*` is the release candidate of the next tobac version. The asterisk stands here for the specific tobac version: RC_vx.x.x (e.g. RC_v1.5.0). Pull requests to this branch need two reviews to be accepted before it can be merged into main. 
+
+:code:`dev_*` is the development branch where we experiment with new features. This branch is perfectly suited to collaboratively work on a feature together with other **tobac** developers (see :doc:`mentoring`). In general, this branch is used for long-term, comprehensive code changes that might not be covered by a single pull request and where it might not be conceivable in which future **tobac** version to include the changes. There are no branch protection rules for this branch, which means that collaborators of our GitHub organization can directly push changes to this branch. Note that **dev_** can never directly merged into main, it has be merged into the release candidate branch :code:`RC_*` first! There can be more than one `dev_*` branch, therefore it we recommend to describe the feature to work on in the respective branch (e.g. :code:`dev_xarray_transition`). 
+
+:code:`hotfix` is the branch we use for hotfixes, i.e. bug fixes that need to be released as fast as possible because it influences people's code. This branch needs only one review before it can directly merged into :code:`main`.
+
+In brief: **Unless you are collaboratively working on a comprehensive feature or on a hotfix, the branch to submit your pull request to is the next release candidate RC_v.x.x.x**
+
+
+=========================
+GitHub workflow
+=========================
+
+We use several `GitHub actions <https://docs.github.com/en/actions/learn-github-actions/understanding-github-actions>`_ to 
+assure continuous integration and to enable an efficient code development and release process. Our workflow 
+configuration can be found in 
+`.github/workflows <https://github.com/tobac-project/tobac/tree/main/.github/workflows>`_ and encompass
+
+* check that code is formatted using the latest stable version of black 
+* linting of the latest code changes that checks the code quality and results in a score compared to the most recent released version 
+* check of the zenodo JSON file that ensures that the citation is correct 
+* check that all unit tests pass (including testing on multiple operating systems) and report test coverage 
+* check that the example jupyter notebooks run without problems 
+
+=========================
+Writing unit tests
+=========================
+
+We use unit tests that ensure that the functions of each module and submodule work properly. If you add a new 
+functionality, you should also add a unit test. All tests are located in the `test 
+folder <https://github.com/tobac-project/tobac/tree/main/tobac/tests>`_ The module :py:mod:`tobac.testing` may help to 
+create simple, idealized cases where objects can be tracked to test if the new features result in the expected outcome. 
+
+If you are unsure on how to construct tests and run tests locally, you can find additional documentation on 
+`pytest <https://docs.pytest.org/en/7.1.x/getting-started.html>`_ and `pytest 
+fixtures <https://docs.pytest.org/en/6.2.x/fixture.html>`_. 
+
+You will also notice that we report the test coverage, i.e. how much of our current code is triggered and thus tested by 
+the unit tests. When you submit a pull request, you will see if your code changes have increased or decreased the test 
+coverage. Ideally, test coverage should not decrease, so please make sure to add appropriate unit tests that cover 
+all newly added functions. 
+
+=========================
+Add examples 
+=========================
+
+In addition to the unit tests, we aim to provide examples on how to use all functionalities and how to choose different 
+tracking parameters. These `examples <https://github.com/tobac-project/tobac/tree/main/examples>`_ are in form of jupyter 
+notebooks and can be based on simple, idealized test cases or real data. We strongly encourage the use of real data that 
+is publicly accessible, but another option for new examples with real data is to either upload the data to our `zenodo 
+repository <https://zenodo.org/records/10863405>`_ or create your own data upload on zenodo. Please include the name "tobac" in the data title for the latter. 
+
+=========================
+Releasing a new version 
+=========================
+
+This is the checklist of steps for a release of a new **tobac** version:
+
+* Bump version in :code:`__init__.py` in :code:`RC_vXXX`
+* Add changelog in :code:`RC_vXXX` 
+* Regenerate example notebooks with the new version
+* Merge :code:`RC_vXXX` into :code:`main` 
+* Merge updated :code:`main` branch back into release and dev branches 
+* Delete :code:`RC_vXXX` branch
+* Create release
+* Push release to conda-forge: https://github.com/tobac-project/tobac-notes/blob/main/uploading_to_conda-forge.md
+* Create new tag
+* E-mail tobac mailing list
+
diff --git a/doc/examples.rst b/doc/examples.rst
index 85160f30..34aac497 100644
--- a/doc/examples.rst
+++ b/doc/examples.rst
@@ -32,4 +32,19 @@ tobac is provided with a set of Jupyter notebooks that show examples of the appl
    WRF Updrafts <examples/Example_Updraft_Tracking/Example_Updraft_Tracking>
    WRF Mesoscale Vorticity <examples/Example_vorticity_tracking_model/Example_vorticity_tracking_model>
 
+
+.. nbgallery::
+   :caption: Calculating Bulk Statistics
+
+   Calculating Bulk Statistics during Feature Detection <bulk_statistics/notebooks/compute_statistics_during_feature_detection>
+   Calculating Bulk Statistics during Segmentation <bulk_statistics/notebooks/compute_statistics_during_segmentation>
+   Calculating Bulk Statistics as a Postprocessing Step <bulk_statistics/notebooks/compute_statistics_postprocessing_example>
+
+
+.. nbgallery::
+   :caption: Examples of Using *tobac* with Large Datasets and in Parallel
+
+   Simple Dask Parallel Tutorial <big_datasets_examples/notebooks/parallel_processing_tobac>
+
+
 The notebooks can be found in the **examples** folder as part of the python package. The necessary input data for these examples is avaliable on zenodo and can be downloaded automatically by the Jupyter notebooks.
\ No newline at end of file
diff --git a/doc/feature_detection/notebooks/feature_detection_filtering.ipynb b/doc/feature_detection/notebooks/feature_detection_filtering.ipynb
index 1dcc998d..44c3d3d4 100644
--- a/doc/feature_detection/notebooks/feature_detection_filtering.ipynb
+++ b/doc/feature_detection/notebooks/feature_detection_filtering.ipynb
@@ -82,13 +82,13 @@
    ],
    "source": [
     "# Dimensions here are time, y, x.\n",
-    "input_field_arr = np.zeros((1,100,200))\n",
-    "input_field_arr[0, 15:85, 10:185]=50\n",
-    "input_field_arr[0, 20:80, 20:80]=100\n",
+    "input_field_arr = np.zeros((1, 100, 200))\n",
+    "input_field_arr[0, 15:85, 10:185] = 50\n",
+    "input_field_arr[0, 20:80, 20:80] = 100\n",
     "input_field_arr[0, 40:60, 125:170] = 100\n",
-    "input_field_arr[0, 30:40, 30:40]=200\n",
-    "input_field_arr[0, 50:75, 50:75]=200\n",
-    "input_field_arr[0, 55:70, 55:70]=300\n",
+    "input_field_arr[0, 30:40, 30:40] = 200\n",
+    "input_field_arr[0, 50:75, 50:75] = 200\n",
+    "input_field_arr[0, 55:70, 55:70] = 300\n",
     "\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
@@ -118,10 +118,14 @@
     }
    ],
    "source": [
-    "# We now need to generate an Iris DataCube out of this dataset to run tobac feature detection. \n",
-    "# One can use xarray to generate a DataArray and then convert it to Iris, as done here. \n",
-    "input_field_iris = xr.DataArray(input_field_arr, dims=['time', 'Y', 'X'], coords={'time': [np.datetime64('2019-01-01T00:00:00')]}).to_iris()\n",
-    "# Version 2.0 of tobac (currently in development) will allow the use of xarray directly with tobac. "
+    "# We now need to generate an Iris DataCube out of this dataset to run tobac feature detection.\n",
+    "# One can use xarray to generate a DataArray and then convert it to Iris, as done here.\n",
+    "input_field_iris = xr.DataArray(\n",
+    "    input_field_arr,\n",
+    "    dims=[\"time\", \"Y\", \"X\"],\n",
+    "    coords={\"time\": [np.datetime64(\"2019-01-01T00:00:00\")]},\n",
+    ").to_iris()\n",
+    "# Version 2.0 of tobac (currently in development) will allow the use of xarray directly with tobac."
    ]
   },
   {
@@ -164,23 +168,36 @@
    ],
    "source": [
     "thresholds = [50, 100, 150, 200]\n",
-    "fig, axarr = plt.subplots(2,2, figsize=(10,6))\n",
+    "fig, axarr = plt.subplots(2, 2, figsize=(10, 6))\n",
     "sigma_values = [0, 1, 2, 5]\n",
     "for sigma_value, ax in zip(sigma_values, axarr.flatten()):\n",
-    "    single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', sigma_threshold=sigma_value)\n",
-    "    \n",
+    "    single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "        field_in=input_field_iris,\n",
+    "        dxy=1000,\n",
+    "        threshold=thresholds,\n",
+    "        target=\"maximum\",\n",
+    "        sigma_threshold=sigma_value,\n",
+    "    )\n",
+    "\n",
     "    # This is what tobac sees\n",
-    "    filtered_field = scipy.ndimage.gaussian_filter(input_field_arr[0], sigma=sigma_value)\n",
+    "    filtered_field = scipy.ndimage.gaussian_filter(\n",
+    "        input_field_arr[0], sigma=sigma_value\n",
+    "    )\n",
     "    color_mesh = ax.pcolormesh(filtered_field)\n",
     "    plt.colorbar(color_mesh, ax=ax)\n",
     "    # Plot all features detected\n",
-    "    ax.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "    ax.scatter(\n",
+    "        x=single_threshold_features[\"hdim_2\"].values,\n",
+    "        y=single_threshold_features[\"hdim_1\"].values,\n",
+    "        color=\"r\",\n",
+    "        label=\"Detected Features\",\n",
+    "    )\n",
     "    ax.legend()\n",
     "    if sigma_value == 0:\n",
     "        sigma_val_str = \"0 (off)\"\n",
     "    else:\n",
     "        sigma_val_str = \"{0}\".format(sigma_value)\n",
-    "    ax.set_title(\"sigma_threshold= \"+ sigma_val_str)\n",
+    "    ax.set_title(\"sigma_threshold= \" + sigma_val_str)\n",
     "plt.tight_layout()\n",
     "plt.show()"
    ]
@@ -225,31 +242,44 @@
    ],
    "source": [
     "thresholds = [100]\n",
-    "fig, axarr = plt.subplots(2,2, figsize=(10,6))\n",
+    "fig, axarr = plt.subplots(2, 2, figsize=(10, 6))\n",
     "erosion_values = [0, 5, 10, 15]\n",
     "for erosion, ax in zip(erosion_values, axarr.flatten()):\n",
-    "    single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', n_erosion_threshold=erosion)\n",
-    "    \n",
+    "    single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "        field_in=input_field_iris,\n",
+    "        dxy=1000,\n",
+    "        threshold=thresholds,\n",
+    "        target=\"maximum\",\n",
+    "        n_erosion_threshold=erosion,\n",
+    "    )\n",
+    "\n",
     "    # Create our mask- this is what tobac does internally for each threshold.\n",
-    "    tobac_mask = 1*(input_field_arr[0] >= thresholds[0])\n",
-    "    \n",
+    "    tobac_mask = 1 * (input_field_arr[0] >= thresholds[0])\n",
+    "\n",
     "    if erosion > 0:\n",
-    "        # This is the parameter for erosion that gets passed to the scikit-image library. \n",
+    "        # This is the parameter for erosion that gets passed to the scikit-image library.\n",
     "        footprint = np.ones((erosion, erosion))\n",
-    "        # This is what tobac sees after erosion. \n",
-    "        filtered_mask = skimage.morphology.binary_erosion(tobac_mask, footprint).astype(np.int64)\n",
+    "        # This is what tobac sees after erosion.\n",
+    "        filtered_mask = skimage.morphology.binary_erosion(tobac_mask, footprint).astype(\n",
+    "            np.int64\n",
+    "        )\n",
     "    else:\n",
     "        filtered_mask = tobac_mask\n",
     "\n",
     "    color_mesh = ax.pcolormesh(filtered_mask)\n",
     "    # Plot all features detected\n",
-    "    ax.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "    ax.scatter(\n",
+    "        x=single_threshold_features[\"hdim_2\"].values,\n",
+    "        y=single_threshold_features[\"hdim_1\"].values,\n",
+    "        color=\"r\",\n",
+    "        label=\"Detected Features\",\n",
+    "    )\n",
     "    ax.legend()\n",
     "    if erosion == 0:\n",
     "        sigma_val_str = \"0 (off)\"\n",
     "    else:\n",
     "        sigma_val_str = \"{0}\".format(erosion)\n",
-    "    ax.set_title(\"n_erosion_threshold= \"+ sigma_val_str)\n",
+    "    ax.set_title(\"n_erosion_threshold= \" + sigma_val_str)\n",
     "plt.tight_layout()\n",
     "plt.show()"
    ]
diff --git a/doc/feature_detection/notebooks/multiple_thresholds_example.ipynb b/doc/feature_detection/notebooks/multiple_thresholds_example.ipynb
index 64f427a6..39bd33ed 100644
--- a/doc/feature_detection/notebooks/multiple_thresholds_example.ipynb
+++ b/doc/feature_detection/notebooks/multiple_thresholds_example.ipynb
@@ -73,13 +73,13 @@
    ],
    "source": [
     "# Dimensions here are time, y, x.\n",
-    "input_field_arr = np.zeros((1,100,200))\n",
-    "input_field_arr[0, 15:85, 10:185]=50\n",
-    "input_field_arr[0, 20:80, 20:80]=100\n",
+    "input_field_arr = np.zeros((1, 100, 200))\n",
+    "input_field_arr[0, 15:85, 10:185] = 50\n",
+    "input_field_arr[0, 20:80, 20:80] = 100\n",
     "input_field_arr[0, 40:60, 125:170] = 100\n",
-    "input_field_arr[0, 30:40, 30:40]=200\n",
-    "input_field_arr[0, 50:75, 50:75]=200\n",
-    "input_field_arr[0, 55:70, 55:70]=300\n",
+    "input_field_arr[0, 30:40, 30:40] = 200\n",
+    "input_field_arr[0, 50:75, 50:75] = 200\n",
+    "input_field_arr[0, 55:70, 55:70] = 300\n",
     "\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
@@ -118,7 +118,11 @@
     }
    ],
    "source": [
-    "input_field_iris = xr.DataArray(input_field_arr, dims=['time', 'Y', 'X'], coords={'time': [np.datetime64('2019-01-01T00:00:00')]}).to_iris()"
+    "input_field_iris = xr.DataArray(\n",
+    "    input_field_arr,\n",
+    "    dims=[\"time\", \"Y\", \"X\"],\n",
+    "    coords={\"time\": [np.datetime64(\"2019-01-01T00:00:00\")]},\n",
+    ").to_iris()"
    ]
   },
   {
@@ -160,13 +164,26 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
     "# Using 'center' here outputs the feature location as the arithmetic center of the detected feature\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold='center')\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=\"center\",\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
     "plt.title(\"Single Threshold of 50\")\n",
     "plt.show()"
@@ -203,13 +220,26 @@
     }
    ],
    "source": [
-    "thresholds = [150,]\n",
+    "thresholds = [\n",
+    "    150,\n",
+    "]\n",
     "# Using 'center' here outputs the feature location as the arithmetic center of the detected feature\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold='center')\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=\"center\",\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
     "plt.title(\"Single Threshold of 150\")\n",
     "plt.show()"
@@ -261,14 +291,27 @@
     }
    ],
    "source": [
-    "thresholds = [100, ]\n",
+    "thresholds = [\n",
+    "    100,\n",
+    "]\n",
     "# Using 'center' here outputs the feature location as the arithmetic center of the detected feature\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold='center')\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=\"center\",\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
     "plt.title(\"Single Threshold of 100\")\n",
     "plt.show()"
@@ -307,12 +350,23 @@
    "source": [
     "thresholds = [50, 100, 150, 200]\n",
     "# Using 'center' here outputs the feature location as the arithmetic center of the detected feature\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold='center')\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=\"center\",\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
     "plt.title(\"Thresholds: [50, 100, 150, 200]\")\n",
     "plt.show()"
@@ -344,17 +398,27 @@
    "source": [
     "thresholds = [50, 100, 150, 200, 250]\n",
     "# Using 'center' here outputs the feature location as the arithmetic center of the detected feature\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold='center')\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=\"center\",\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "\n",
     "# Plot all features detected\n",
     "for i, threshold in enumerate(thresholds):\n",
-    "    thresholded_points = single_threshold_features[single_threshold_features['threshold_value'] == threshold]\n",
-    "    plt.scatter(x=thresholded_points['hdim_2'].values, \n",
-    "                y=thresholded_points['hdim_1'].values, \n",
-    "                color='C'+str(i),\n",
-    "                label=\"Threshold: \"+str(threshold))\n",
+    "    thresholded_points = single_threshold_features[\n",
+    "        single_threshold_features[\"threshold_value\"] == threshold\n",
+    "    ]\n",
+    "    plt.scatter(\n",
+    "        x=thresholded_points[\"hdim_2\"].values,\n",
+    "        y=thresholded_points[\"hdim_1\"].values,\n",
+    "        color=\"C\" + str(i),\n",
+    "        label=\"Threshold: \" + str(threshold),\n",
+    "    )\n",
     "plt.legend()\n",
     "plt.title(\"Thresholds: [50, 100, 150, 200]\")\n",
     "plt.show()"
diff --git a/doc/feature_detection/notebooks/position_threshold_example.ipynb b/doc/feature_detection/notebooks/position_threshold_example.ipynb
index 37a1231b..016fb245 100644
--- a/doc/feature_detection/notebooks/position_threshold_example.ipynb
+++ b/doc/feature_detection/notebooks/position_threshold_example.ipynb
@@ -73,13 +73,13 @@
    ],
    "source": [
     "# Dimensions here are time, y, x.\n",
-    "input_field_arr = np.zeros((1,100,200))\n",
-    "input_field_arr[0, 15:85, 10:185]=50\n",
-    "input_field_arr[0, 20:80, 20:80]=100\n",
+    "input_field_arr = np.zeros((1, 100, 200))\n",
+    "input_field_arr[0, 15:85, 10:185] = 50\n",
+    "input_field_arr[0, 20:80, 20:80] = 100\n",
     "input_field_arr[0, 40:60, 125:170] = 100\n",
-    "input_field_arr[0, 30:40, 30:40]=200\n",
-    "input_field_arr[0, 50:75, 50:75]=200\n",
-    "input_field_arr[0, 55:70, 55:70]=300\n",
+    "input_field_arr[0, 30:40, 30:40] = 200\n",
+    "input_field_arr[0, 50:75, 50:75] = 200\n",
+    "input_field_arr[0, 55:70, 55:70] = 300\n",
     "\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
@@ -109,10 +109,14 @@
     }
    ],
    "source": [
-    "# We now need to generate an Iris DataCube out of this dataset to run tobac feature detection. \n",
-    "# One can use xarray to generate a DataArray and then convert it to Iris, as done here. \n",
-    "input_field_iris = xr.DataArray(input_field_arr, dims=['time', 'Y', 'X'], coords={'time': [np.datetime64('2019-01-01T00:00:00')]}).to_iris()\n",
-    "# Version 2.0 of tobac (currently in development) will allow the use of xarray directly with tobac. "
+    "# We now need to generate an Iris DataCube out of this dataset to run tobac feature detection.\n",
+    "# One can use xarray to generate a DataArray and then convert it to Iris, as done here.\n",
+    "input_field_iris = xr.DataArray(\n",
+    "    input_field_arr,\n",
+    "    dims=[\"time\", \"Y\", \"X\"],\n",
+    "    coords={\"time\": [np.datetime64(\"2019-01-01T00:00:00\")]},\n",
+    ").to_iris()\n",
+    "# Version 2.0 of tobac (currently in development) will allow the use of xarray directly with tobac."
    ]
   },
   {
@@ -154,15 +158,28 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
-    "position_threshold = 'center'\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold=position_threshold)\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
+    "position_threshold = \"center\"\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=position_threshold,\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
-    "plt.title(\"position_threshold \"+ position_threshold)\n",
+    "plt.title(\"position_threshold \" + position_threshold)\n",
     "plt.show()"
    ]
   },
@@ -205,15 +222,28 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
-    "position_threshold = 'extreme'\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold=position_threshold)\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
+    "position_threshold = \"extreme\"\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=position_threshold,\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
-    "plt.title(\"position_threshold \"+ position_threshold)\n",
+    "plt.title(\"position_threshold \" + position_threshold)\n",
     "plt.show()"
    ]
   },
@@ -256,15 +286,28 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
-    "position_threshold = 'weighted_diff'\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold=position_threshold)\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
+    "position_threshold = \"weighted_diff\"\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=position_threshold,\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
-    "plt.title(\"position_threshold \"+ position_threshold)\n",
+    "plt.title(\"position_threshold \" + position_threshold)\n",
     "plt.show()"
    ]
   },
@@ -307,15 +350,28 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
-    "position_threshold = 'weighted_abs'\n",
-    "single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold=position_threshold)\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
+    "position_threshold = \"weighted_abs\"\n",
+    "single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "    field_in=input_field_iris,\n",
+    "    dxy=1000,\n",
+    "    threshold=thresholds,\n",
+    "    target=\"maximum\",\n",
+    "    position_threshold=position_threshold,\n",
+    ")\n",
     "plt.pcolormesh(input_field_arr[0])\n",
     "plt.colorbar()\n",
     "# Plot all features detected\n",
-    "plt.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "plt.scatter(\n",
+    "    x=single_threshold_features[\"hdim_2\"].values,\n",
+    "    y=single_threshold_features[\"hdim_1\"].values,\n",
+    "    color=\"r\",\n",
+    "    label=\"Detected Features\",\n",
+    ")\n",
     "plt.legend()\n",
-    "plt.title(\"position_threshold \"+ position_threshold)\n",
+    "plt.title(\"position_threshold \" + position_threshold)\n",
     "plt.show()"
    ]
   },
@@ -357,18 +413,30 @@
     }
    ],
    "source": [
-    "thresholds = [50,]\n",
-    "fig, axarr = plt.subplots(2,2, figsize=(10,6))\n",
-    "testing_thresholds = ['center', 'extreme', 'weighted_diff', 'weighted_abs']\n",
+    "thresholds = [\n",
+    "    50,\n",
+    "]\n",
+    "fig, axarr = plt.subplots(2, 2, figsize=(10, 6))\n",
+    "testing_thresholds = [\"center\", \"extreme\", \"weighted_diff\", \"weighted_abs\"]\n",
     "for position_threshold, ax in zip(testing_thresholds, axarr.flatten()):\n",
-    "\n",
-    "    single_threshold_features = tobac.feature_detection_multithreshold(field_in = input_field_iris, dxy = 1000, threshold=thresholds, target='maximum', position_threshold=position_threshold)\n",
+    "    single_threshold_features = tobac.feature_detection_multithreshold(\n",
+    "        field_in=input_field_iris,\n",
+    "        dxy=1000,\n",
+    "        threshold=thresholds,\n",
+    "        target=\"maximum\",\n",
+    "        position_threshold=position_threshold,\n",
+    "    )\n",
     "    color_mesh = ax.pcolormesh(input_field_arr[0])\n",
     "    plt.colorbar(color_mesh, ax=ax)\n",
     "    # Plot all features detected\n",
-    "    ax.scatter(x=single_threshold_features['hdim_2'].values, y=single_threshold_features['hdim_1'].values, color='r', label=\"Detected Features\")\n",
+    "    ax.scatter(\n",
+    "        x=single_threshold_features[\"hdim_2\"].values,\n",
+    "        y=single_threshold_features[\"hdim_1\"].values,\n",
+    "        color=\"r\",\n",
+    "        label=\"Detected Features\",\n",
+    "    )\n",
     "    ax.legend()\n",
-    "    ax.set_title(\"position_threshold \"+ position_threshold)\n",
+    "    ax.set_title(\"position_threshold \" + position_threshold)\n",
     "plt.tight_layout()\n",
     "plt.show()"
    ]
diff --git a/doc/index.rst b/doc/index.rst
index d008e946..1b215f20 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -67,12 +67,20 @@ The project is currently being extended by several contributors to include addit
 
 
 .. toctree::
+   :caption: Developer guide
+   :maxdepth: 3
+
+   code_structure
+   contributing
+   code_reviews
+   mentoring
+
+.. toctree:: 
    :caption: Compute bulk statistics
    :maxdepth: 2
 
    bulk_statistics/index
 
-   
 .. toctree::
    :caption: API Reference
    :maxdepth: 3
diff --git a/doc/mentoring.rst b/doc/mentoring.rst
new file mode 100644
index 00000000..c956af5e
--- /dev/null
+++ b/doc/mentoring.rst
@@ -0,0 +1,25 @@
+Mentoring and Collaboration
+----------------------------
+
+============================
+Writing code collaboratively
+============================
+
+We firmly believe that code can only get better if more than two eyes and one brain work on it. Therefore, we aim to write code collaboratively, in particular, when comprehensive refactoring or enhancements of the code are done. In practice, this can be done by creating a **draft pull request**. This makes it really easy to iteratively improve a pull request with the feedback from others until the pull request is ready for review.
+
+
+When you work on a comprehensive feature with multiple developers, it is recommended to create a draft pull request on the :code:`dev_*` branch. As explained in :doc:`our branching strategy <contributing>`, this branch does not undergo any protection rules. It is meant to experiment with new code and all collaborators of the `tobac-project organization <https://github.com/tobac-project>`_ can directly push to this branch. Creating a draft pull request has the advantage of facilitating the communication with other developers who contribute to the same new feature. You can directly see which changes they make, comment these and discuss ways to go forward.
+
+==============
+Get a mentor
+==============
+
+**Is this your first time contributing to an open-source project?**
+
+Reach out to the **tobac** developer group and get a mentor! One of our developers will help you getting started and explain how our workflow works. You are, of course, always free to post any questions to GitHub discussions, our Slack channel, or write an email. But sometimes it can also be nice to have a specific person to refer to when things seem overwhelming in the beginning. 
+
+===============
+Pair reviews 
+===============
+
+Another great way of collaboration are pair reviews which means that you are reviewing code together with another developer. You can, for example, reach out to us when you have submitted a pull request and would like to talk through the review points with one of the reviewers in order to collaboratively come up with creative solutions to remaining issues. If you are a reviewer, you can offer a pair review to the person who created the pull request and help them addressing certain review points. 
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 550ba3d6..99bb7e46 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -3,3 +3,4 @@ nbsphinx
 numpy
 sphinx_rtd_theme
 sphinx-gallery
+sphinx-toolbox
diff --git a/doc/testing_sphinx-based_rendering.rst b/doc/testing_sphinx-based_rendering.rst
new file mode 100644
index 00000000..f0a7443e
--- /dev/null
+++ b/doc/testing_sphinx-based_rendering.rst
@@ -0,0 +1,154 @@
+How to check the Sphinx-based rendering
+---------------------------------------
+
+
+The workflow has been tested in a linux system. We aim to build a static
+website out of the documentation material present in ``tobac``.
+
+==================================
+1. Preparing the Local Environment
+==================================
+
+-  **choose a separate place for your testing**
+
+   I will use the temporary directory ``/tmp/website-testing`` which I
+   need to create. You can use a dedicated place of your choice …
+
+   .. code:: bash
+
+      > mkdir /tmp/website-testing
+      > cd /tmp/website-testing
+
+   I will indicate my position now with the ``/tmp/website-testing>``
+   prompt.
+
+-  **get the official repository**
+
+   .. code:: bash
+
+      /tmp/website-testing> git clone https://github.com/tobac-project/tobac
+
+   You might like to test a certain remote branch ``<BRANCH>`` then do:
+
+   .. code:: bash
+
+      /tmp/website-testing/tobac> git fetch --all
+      /tmp/website-testing/tobac> git checkout -t origin/<BRANCH> 
+
+-  **Python environment**
+
+   -  create a python virtual env
+
+      .. code:: bash
+      
+         /tmp/website-testing> python -m venv .python3-venv
+
+
+   -  and install requirements
+      
+      .. code:: bash
+         
+         # deactivation conda is only necessary if your loaded conda before … 
+         /tmp/website-testing> conda deactivate
+
+         # activate the new env and upgrade ``pip`` 
+         /tmp/website-testing> source .python3-venv/bin/activate 
+         /tmp/website-testing> pip install –upgrade pip
+
+         # now everything is installed into the local python env!
+         /tmp/website-testing> pip install -r tobac/doc/requirements.txt
+
+         # and also install RTD scheme 
+         /tmp/website-testing> pip install sphinx_rtd_theme
+
+      `pip`-based installation takes a bit of time, but is much faster than `conda`.
+   
+
+If the installation runs without problems, you are ready to build the website.
+
+
+==================================
+1. Building the Website
+==================================
+
+Actually, only few steps are needed to build the website, i.e.
+
+-  **running sphinx for rendering**
+
+   .. code:: bash
+
+      /tmp/website-testing> cd tobac
+
+      /tmp/website-testing/tobac> sphinx-build -b html doc doc/_build/html
+
+   If no severe error appeared
+
+-  **view the HTML content**
+
+   .. code:: bash
+
+      /tmp/website-testing/tobac> firefox doc/_build/html/index.html
+
+==================================
+3. Parsing Your Local Changes
+==================================
+
+Now, we connect to your locally hosted ``tobac`` repository and your
+development branch.
+
+-  **connect to your local repo**: Assume your repo is located at
+   ``/tmp/tobac-testing/tobac``, then add a new remote alias and fetch
+   all content with
+
+   .. code:: bash
+
+      /tmp/website-testing/tobac> git remote add local-repo /tmp/tobac-testing/tobac
+      /tmp/website-testing/tobac> git fetch --all
+
+-  **check your development branch out**: Now, assume the your
+   development branch is called ``my-devel``, then do
+
+   .. code:: bash
+
+      # to get a first overview on available branches
+      /tmp/website-testing/tobac> git branch --all
+
+      # and then actually get your development branch
+      /tmp/website-testing/tobac> git checkout -b my-devel local-repo/my-devel
+
+   You should see your developments, now …
+
+-  **build and view website again**
+
+   .. code:: bash
+
+      /tmp/website-testing/tobac> sphinx-build -M clean doc doc/_build
+      /tmp/website-testing/tobac> sphinx-build -b html doc doc/_build/html
+      /tmp/website-testing/tobac> firefox _build/html/index.html
+
+
+==========================================
+Option: Check Rendering of a Pull requests
+==========================================
+
+-  **check the pull request out**: Now, assume the PR has the ID ``<ID>`` and you define the branch name ``BRANCH_NAME`` as you like
+
+   .. code:: bash
+
+      # to get PR shown as dedicated branch
+      /tmp/website-testing/tobac> git fetch upstream pull/ID/head:BRANCH_NAME
+
+      # and then actually get this PR as branch
+      /tmp/website-testing/tobac> git checkout BRANCH_NAME
+
+   You should see the PR now ...
+
+-  **build and view website again**
+
+   .. code:: bash
+
+      /tmp/website-testing/tobac> sphinx-build -M clean doc doc/_build
+      /tmp/website-testing/tobac> sphinx-build -b html doc doc/_build/html
+      /tmp/website-testing/tobac> firefox _build/html/index.html
+
+
diff --git a/tobac/feature_detection.py b/tobac/feature_detection.py
index b937a0c0..fc491ff0 100644
--- a/tobac/feature_detection.py
+++ b/tobac/feature_detection.py
@@ -18,7 +18,7 @@
 """
 
 from __future__ import annotations
-from typing import Union, Callable
+from typing import Optional, Union, Callable
 import warnings
 import logging
 
@@ -199,8 +199,8 @@ def feature_position(
 
     elif position_threshold == "weighted_diff":
         # get position as centre of identified region, weighted by difference from the threshold:
-        weights = abs(track_data_region[region_small] - threshold_i)
-        if sum(weights) == 0:
+        weights = np.abs(track_data_region[region_small] - threshold_i)
+        if np.sum(weights) == 0:
             weights = None
         hdim1_weights = weights
         hdim2_weights = weights
@@ -211,8 +211,8 @@ def feature_position(
 
     elif position_threshold == "weighted_abs":
         # get position as centre of identified region, weighted by absolute values if the field:
-        weights = abs(track_data_region[region_small])
-        if sum(weights) == 0:
+        weights = np.abs(track_data_region[region_small])
+        if np.sum(weights) == 0:
             weights = None
         hdim1_weights = weights
         hdim2_weights = weights
@@ -230,14 +230,18 @@ def feature_position(
             hdim1_index = pbc_utils.weighted_circmean(
                 hdim1_indices, weights=hdim1_weights, high=hdim1_max + 1, low=hdim1_min
             )
+            hdim1_index = np.clip(hdim1_index, 0, hdim1_max + 1)
         else:
             hdim1_index = np.average(hdim1_indices, weights=hdim1_weights)
+            hdim1_index = np.clip(hdim1_index, 0, hdim1_max)
         if PBC_flag in ("hdim_2", "both"):
             hdim2_index = pbc_utils.weighted_circmean(
                 hdim2_indices, weights=hdim2_weights, high=hdim2_max + 1, low=hdim2_min
             )
+            hdim2_index = np.clip(hdim2_index, 0, hdim2_max + 1)
         else:
             hdim2_index = np.average(hdim2_indices, weights=hdim2_weights)
+            hdim2_index = np.clip(hdim2_index, 0, hdim2_max)
         if is_3D:
             vdim_index = np.average(vdim_indices, weights=vdim_weights)
 
@@ -390,7 +394,7 @@ def feature_detection_threshold(
     ] = "center",
     sigma_threshold: float = 0.5,
     n_erosion_threshold: int = 0,
-    n_min_threshold: int = 0,
+    n_min_threshold: Union[int, dict[float, int], list[int]] = 0,
     min_distance: float = 0,
     idx_start: int = 0,
     PBC_flag: Literal["none", "hdim_1", "hdim_2", "both"] = "none",
@@ -428,8 +432,10 @@ def feature_detection_threshold(
         Number of pixels by which to erode the identified features.
         Default is 0.
 
-    n_min_threshold : int, optional
+    n_min_threshold : int, dict of float to int, or list of int, optional
         Minimum number of identified contiguous pixels for a feature to be detected. Default is 0.
+        If given as a list, the number of elements must match number of thresholds.
+        If given as a dict, the keys need to match the thresholds and the values are the minimum number of identified contiguous pixels for a feature using that specific threshold.
 
     min_distance : float, optional
         Minimum distance between detected features (in meters). Default is 0.
@@ -897,7 +903,7 @@ def feature_detection_multithreshold_timestep(
     ] = "center",
     sigma_threshold: float = 0.5,
     n_erosion_threshold: int = 0,
-    n_min_threshold: int = 0,
+    n_min_threshold: Union[int, dict[float, int], list[int]] = 0,
     min_distance: float = 0,
     feature_number_start: int = 1,
     PBC_flag: Literal["none", "hdim_1", "hdim_2", "both"] = "none",
@@ -944,8 +950,10 @@ def feature_detection_multithreshold_timestep(
         Number of pixels by which to erode the identified features.
         Default is 0.
 
-    n_min_threshold : int, optional
+    n_min_threshold :  int, dict of float to int, or list of int, optional
         Minimum number of identified contiguous pixels for a feature to be detected. Default is 0.
+        If given as a list, the number of elements must match number of thresholds.
+        If given as a dict, the keys need to match the thresholds and the values are the minimum number of identified contiguous pixels for a feature using that specific threshold.
 
     min_distance : float, optional
         Minimum distance between detected features (in meters). Default is 0.
@@ -1026,7 +1034,7 @@ def feature_detection_multithreshold_timestep(
             ):
                 raise ValueError(
                     "Ambiguous input for threshold values. If n_min_threshold is given as a dict,"
-                    " the keys not to correspond to the values in threshold."
+                    " the keys must correspond to the values in threshold."
                 )
             # sort dictionary by keys (threshold values) so that they match sorted thresholds and
             # get values for n_min_threshold
@@ -1139,14 +1147,14 @@ def feature_detection_multithreshold(
     ] = "center",
     sigma_threshold: float = 0.5,
     n_erosion_threshold: int = 0,
-    n_min_threshold: int = 0,
+    n_min_threshold: Union[int, dict[float, int], list[int]] = 0,
     min_distance: float = 0,
     feature_number_start: int = 1,
     PBC_flag: Literal["none", "hdim_1", "hdim_2", "both"] = "none",
-    vertical_coord: str = None,
-    vertical_axis: int = None,
-    detect_subset: dict = None,
-    wavelength_filtering: tuple = None,
+    vertical_coord: Optional[str] = None,
+    vertical_axis: Optional[int] = None,
+    detect_subset: Optional[dict] = None,
+    wavelength_filtering: Optional[tuple] = None,
     dz: Union[float, None] = None,
     strict_thresholding: bool = False,
     statistic: Union[dict[str, Union[Callable, tuple[Callable, dict]]], None] = None,
@@ -1184,8 +1192,10 @@ def feature_detection_multithreshold(
         Number of pixels by which to erode the identified features.
         Default is 0.
 
-    n_min_threshold : int, optional
+    n_min_threshold :  int, dict of float to int, or list of int, optional
         Minimum number of identified contiguous pixels for a feature to be detected. Default is 0.
+        If given as a list, the number of elements must match number of thresholds.
+        If given as a dict, the keys need to match the thresholds and the values are the minimum number of identified contiguous pixels for a feature using that specific threshold.
 
     min_distance : float, optional
         Minimum distance between detected features (in meters). Default is 0.
diff --git a/tobac/utils/bulk_statistics.py b/tobac/utils/bulk_statistics.py
index cb9ce0cf..2dbcf6f2 100644
--- a/tobac/utils/bulk_statistics.py
+++ b/tobac/utils/bulk_statistics.py
@@ -122,7 +122,7 @@ def get_statistics(
     # mask must contain positive values to calculate statistics
     if np.any(labels > 0):
         if index is None:
-            index = features.feature.to_numpy()
+            index = features[id_column].to_numpy().astype(int)
         else:
             # get the statistics only for specified feature objects
             if np.max(index) > np.max(labels):
@@ -266,10 +266,16 @@ def get_statistics_from_mask(
          Updated feature dataframe with bulk statistics for each feature saved in a new column
     """
     # warning when feature labels are not unique in dataframe
-    if not features.feature.is_unique:
-        raise logging.warning(
+    if not features[id_column].is_unique:
+        logging.warning(
             "Feature labels are not unique which may cause unexpected results for the computation of bulk statistics."
         )
+    # extra warning when feature labels are not unique in timestep
+    uniques = features.groupby("time")[id_column].value_counts().values
+    if not uniques[uniques > 1].size == 0:
+        logging.warning(
+            "Note that non-unique feature labels occur also in the same timestep. This likely causes unexpected results for the computation of bulk statistics."
+        )
 
     if collapse_dim is not None:
         if isinstance(collapse_dim, str):
@@ -299,7 +305,7 @@ def get_statistics_from_mask(
 
         # make sure that the labels in the segmentation mask exist in feature dataframe
         if (
-            np.intersect1d(np.unique(segmentation_mask_t), features_t.feature).size
+            np.intersect1d(np.unique(segmentation_mask_t), features_t[id_column]).size
             > np.unique(segmentation_mask_t).size
         ):
             raise ValueError(
diff --git a/tobac/utils/decorators.py b/tobac/utils/decorators.py
index 8bc6657f..90e600b5 100644
--- a/tobac/utils/decorators.py
+++ b/tobac/utils/decorators.py
@@ -72,11 +72,13 @@ def _conv_kwargs_irispandas_to_xarray(conv_kwargs: dict):
 
     """
     return {
-        key: convert_cube_to_dataarray(arg)
-        if isinstance(arg, iris.cube.Cube)
-        else arg.to_xarray()
-        if isinstance(arg, pd.DataFrame)
-        else arg
+        key: (
+            convert_cube_to_dataarray(arg)
+            if isinstance(arg, iris.cube.Cube)
+            else arg.to_xarray()
+            if isinstance(arg, pd.DataFrame)
+            else arg
+        )
         for key, arg in zip(conv_kwargs.keys(), conv_kwargs.values())
     }
 
@@ -118,11 +120,13 @@ def _conv_kwargs_xarray_to_irispandas(conv_kwargs: dict):
         iris cubes
     """
     return {
-        key: xr.DataArray.to_iris(arg)
-        if isinstance(arg, xr.DataArray)
-        else arg.to_dataframe()
-        if isinstance(arg, xr.Dataset)
-        else arg
+        key: (
+            xr.DataArray.to_iris(arg)
+            if isinstance(arg, xr.DataArray)
+            else arg.to_dataframe()
+            if isinstance(arg, xr.Dataset)
+            else arg
+        )
         for key, arg in zip(conv_kwargs.keys(), conv_kwargs.values())
     }
 
@@ -166,9 +170,11 @@ def wrapper(*args, **kwargs):
                 # print("converting iris to xarray and back")
                 args = tuple(
                     [
-                        convert_cube_to_dataarray(arg)
-                        if type(arg) == iris.cube.Cube
-                        else arg
+                        (
+                            convert_cube_to_dataarray(arg)
+                            if type(arg) == iris.cube.Cube
+                            else arg
+                        )
                         for arg in args
                     ]
                 )
@@ -179,9 +185,11 @@ def wrapper(*args, **kwargs):
                 if type(output) == tuple:
                     output = tuple(
                         [
-                            xarray.DataArray.to_iris(output_item)
-                            if type(output_item) == xarray.DataArray
-                            else output_item
+                            (
+                                xarray.DataArray.to_iris(output_item)
+                                if type(output_item) == xarray.DataArray
+                                else output_item
+                            )
                             for output_item in output
                         ]
                     )
@@ -241,9 +249,11 @@ def wrapper(*args, **kwargs):
                 # print("converting xarray to iris and back")
                 args = tuple(
                     [
-                        xarray.DataArray.to_iris(arg)
-                        if type(arg) == xarray.DataArray
-                        else arg
+                        (
+                            xarray.DataArray.to_iris(arg)
+                            if type(arg) == xarray.DataArray
+                            else arg
+                        )
                         for arg in args
                     ]
                 )
@@ -257,9 +267,11 @@ def wrapper(*args, **kwargs):
                 if type(output) == tuple:
                     output = tuple(
                         [
-                            xarray.DataArray.from_iris(output_item)
-                            if type(output_item) == iris.cube.Cube
-                            else output_item
+                            (
+                                xarray.DataArray.from_iris(output_item)
+                                if type(output_item) == iris.cube.Cube
+                                else output_item
+                            )
                             for output_item in output
                         ]
                     )
@@ -325,11 +337,13 @@ def wrapper(*args, **kwargs):
                 # print("converting iris to xarray and back")
                 args = tuple(
                     [
-                        convert_cube_to_dataarray(arg)
-                        if type(arg) == iris.cube.Cube
-                        else arg.to_xarray()
-                        if type(arg) == pd.DataFrame
-                        else arg
+                        (
+                            convert_cube_to_dataarray(arg)
+                            if type(arg) == iris.cube.Cube
+                            else arg.to_xarray()
+                            if type(arg) == pd.DataFrame
+                            else arg
+                        )
                         for arg in args
                     ]
                 )
@@ -339,11 +353,15 @@ def wrapper(*args, **kwargs):
                 if type(output) == tuple:
                     output = tuple(
                         [
-                            xarray.DataArray.to_iris(output_item)
-                            if type(output_item) == xarray.DataArray
-                            else output_item.to_dataframe()
-                            if type(output_item) == xarray.Dataset
-                            else output_item
+                            (
+                                xarray.DataArray.to_iris(output_item)
+                                if type(output_item) == xarray.DataArray
+                                else (
+                                    output_item.to_dataframe()
+                                    if type(output_item) == xarray.Dataset
+                                    else output_item
+                                )
+                            )
                             for output_item in output
                         ]
                     )
@@ -415,11 +433,15 @@ def wrapper(*args, **kwargs):
                 # print("converting xarray to iris and back")
                 args = tuple(
                     [
-                        xarray.DataArray.to_iris(arg)
-                        if type(arg) == xarray.DataArray
-                        else arg.to_dataframe()
-                        if type(arg) == xarray.Dataset
-                        else arg
+                        (
+                            xarray.DataArray.to_iris(arg)
+                            if type(arg) == xarray.DataArray
+                            else (
+                                arg.to_dataframe()
+                                if type(arg) == xarray.Dataset
+                                else arg
+                            )
+                        )
                         for arg in args
                     ]
                 )
@@ -433,11 +455,15 @@ def wrapper(*args, **kwargs):
                 if type(output) == tuple:
                     output = tuple(
                         [
-                            xarray.DataArray.from_iris(output_item)
-                            if type(output_item) == iris.cube.Cube
-                            else output_item.to_xarray()
-                            if type(output_item) == pd.DataFrame
-                            else output_item
+                            (
+                                xarray.DataArray.from_iris(output_item)
+                                if type(output_item) == iris.cube.Cube
+                                else (
+                                    output_item.to_xarray()
+                                    if type(output_item) == pd.DataFrame
+                                    else output_item
+                                )
+                            )
                             for output_item in output
                         ]
                     )
diff --git a/tobac/utils/general.py b/tobac/utils/general.py
index 0d633edf..394034c4 100644
--- a/tobac/utils/general.py
+++ b/tobac/utils/general.py
@@ -4,6 +4,8 @@
 
 import copy
 import logging
+from typing_extensions import Literal
+import iris
 import pandas as pd
 
 from . import internal as internal_utils
@@ -16,13 +18,16 @@
 import warnings
 
 
-def add_coordinates(t, variable_cube):
+def add_coordinates(
+    features: pd.DataFrame,
+    variable_cube: iris.cube.Cube,
+):
     """Add coordinates from the input cube of the feature detection
     to the trajectories/features.
 
     Parameters
     ----------
-    t : pandas.DataFrame
+    features : pandas.DataFrame
         Trajectories/features from feature detection or linking step.
 
     variable_cube : iris.cube.Cube
@@ -32,7 +37,7 @@ def add_coordinates(t, variable_cube):
 
     Returns
     -------
-    t : pandas.DataFrame
+    pandas.DataFrame
         Trajectories with added coordinates.
 
     """
@@ -42,17 +47,17 @@ def add_coordinates(t, variable_cube):
     logging.debug("start adding coordinates from cube")
 
     # pull time as datetime object and timestr from input data and add it to DataFrame:
-    t["time"] = None
-    t["timestr"] = None
+    features["time"] = None
+    features["timestr"] = None
 
     logging.debug("adding time coordinate")
 
     time_in = variable_cube.coord("time")
     time_in_datetime = time_in.units.num2date(time_in.points)
 
-    t["time"] = time_in_datetime[t["frame"]]
-    t["timestr"] = [
-        x.strftime("%Y-%m-%d %H:%M:%S") for x in time_in_datetime[t["frame"]]
+    features["time"] = time_in_datetime[features["frame"]]
+    features["timestr"] = [
+        x.strftime("%Y-%m-%d %H:%M:%S") for x in time_in_datetime[features["frame"]]
     ]
 
     # Get list of all coordinates in input cube except for time (already treated):
@@ -88,7 +93,7 @@ def add_coordinates(t, variable_cube):
                     variable_cube.coord(coord).points,
                     fill_value="extrapolate",
                 )
-                coordinate_points = f(t["hdim_1"])
+                coordinate_points = f(features["hdim_1"])
 
             if variable_cube.coord_dims(coord) == (hdim_2,):
                 f = interp1d(
@@ -96,21 +101,21 @@ def add_coordinates(t, variable_cube):
                     variable_cube.coord(coord).points,
                     fill_value="extrapolate",
                 )
-                coordinate_points = f(t["hdim_2"])
+                coordinate_points = f(features["hdim_2"])
 
         # interpolate 2D coordinates:
         elif variable_cube.coord(coord).ndim == 2:
             if variable_cube.coord_dims(coord) == (hdim_1, hdim_2):
                 points = (dimvec_1, dimvec_2)
                 values = variable_cube.coord(coord).points
-                xi = np.column_stack((t["hdim_1"], t["hdim_2"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (hdim_2, hdim_1):
                 points = (dimvec_2, dimvec_1)
                 values = variable_cube.coord(coord).points
-                xi = np.column_stack((t["hdim_2"], t["hdim_1"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
         # interpolate 3D coordinates:
         # mainly workaround for wrf latitude and longitude (to be fixed in future)
@@ -119,44 +124,44 @@ def add_coordinates(t, variable_cube):
             if variable_cube.coord_dims(coord) == (ndim_time, hdim_1, hdim_2):
                 points = (dimvec_1, dimvec_2)
                 values = variable_cube[0, :, :].coord(coord).points
-                xi = np.column_stack((t["hdim_1"], t["hdim_2"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (ndim_time, hdim_2, hdim_1):
                 points = (dimvec_2, dimvec_1)
                 values = variable_cube[0, :, :].coord(coord).points
-                xi = np.column_stack((t["hdim_2"], t["hdim_1"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (hdim_1, ndim_time, hdim_2):
                 points = (dimvec_1, dimvec_2)
                 values = variable_cube[:, 0, :].coord(coord).points
-                xi = np.column_stack((t["hdim_1"], t["hdim_2"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (hdim_1, hdim_2, ndim_time):
                 points = (dimvec_1, dimvec_2)
                 values = variable_cube[:, :, 0].coord(coord).points
-                xi = np.column_stack((t["hdim_1"], t["hdim_2"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (hdim_2, ndim_time, hdim_1):
                 points = (dimvec_2, dimvec_1)
                 values = variable_cube[:, 0, :].coord(coord).points
-                xi = np.column_stack((t["hdim_2"], t["hdim_1"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
             if variable_cube.coord_dims(coord) == (hdim_2, hdim_1, ndim_time):
                 points = (dimvec_2, dimvec_1)
                 values = variable_cube[:, :, 0].coord(coord).points
-                xi = np.column_stack((t["hdim_2"], t["hdim_1"]))
-                coordinate_points = interpn(points, values, xi)
+                xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
+                coordinate_points = interpn(points, values, xi, bounds_error=False)
 
         # write resulting array or list into DataFrame:
-        t[coord] = coordinate_points
+        features[coord] = coordinate_points
 
         logging.debug("added coord: " + coord)
-    return t
+    return features
 
 
 def add_coordinates_3D(
diff --git a/tobac/utils/periodic_boundaries.py b/tobac/utils/periodic_boundaries.py
index 71ecb38e..e230aca5 100644
--- a/tobac/utils/periodic_boundaries.py
+++ b/tobac/utils/periodic_boundaries.py
@@ -1,5 +1,6 @@
 """Utilities for handling indexing and distance calculation with periodic boundaries
 """
+
 from __future__ import annotations
 import functools
 


Test patch:
diff --git a/tobac/tests/test_decorators.py b/tobac/tests/test_decorators.py
index 01a3a0ad..039babde 100644
--- a/tobac/tests/test_decorators.py
+++ b/tobac/tests/test_decorators.py
@@ -1,6 +1,7 @@
 """
 Tests for tobac.utils.decorators
 """
+
 import numpy as np
 import pandas as pd
 import xarray as xr
diff --git a/tobac/tests/test_utils.py b/tobac/tests/test_utils.py
index 439dd79b..e9953d7f 100644
--- a/tobac/tests/test_utils.py
+++ b/tobac/tests/test_utils.py
@@ -584,3 +584,63 @@ def test_transform_feature_points_3D():
     assert np.all(new_feat_df["hdim_1"] == [25, 30])
     assert np.all(new_feat_df["hdim_2"] == [5, 15])
     assert np.all(new_feat_df["vdim"] == [5, 10])
+
+
+def test_get_spacings():
+    """Tests tobac.utils.get_spacings."""
+
+    x_values = np.linspace(100, 500, 5)
+    y_values = np.linspace(400, 200, 5)
+    t_values = np.array([0, 1, 2])
+
+    in_xr = xr.DataArray(
+        np.zeros((len(t_values), len(y_values), len(x_values))),
+        dims=["time", "y", "x"],
+        coords={
+            "x": (
+                "x",
+                x_values,
+                {"units": "meters", "standard_name": "projection_x_coordinate"},
+            ),
+            "y": (
+                "y",
+                y_values,
+                {"units": "meters", "standard_name": "projection_y_coordinate"},
+            ),
+            "time": (
+                "time",
+                t_values,
+                {"units": "hours since 1970-01-01 00:00:00", "standard_name": "time"},
+            ),
+        },
+    )
+
+    # Test with arithmetic average and different dx and dy
+    dxy, dt = tb_utils.get_spacings(in_xr)
+    assert dxy == (100 + 50) / 2
+    assert dt == 3600
+
+    # Test with geometric average and different dx and dy
+    dxy, _ = tb_utils.get_spacings(in_xr, average_method="geometric")
+    assert dxy == np.sqrt(100 * 50)
+
+    # Test with specified grid spacing and time spacing
+    dxy, dt = tb_utils.get_spacings(in_xr, grid_spacing=15, time_spacing=1800)
+    assert dxy == 15
+    assert dt == 1800
+
+    in_xr = xr.DataArray(
+        np.zeros((len(t_values), len(y_values), len(x_values))),
+        dims=["time", "y", "x"],
+        coords={
+            "time": (
+                "time",
+                t_values,
+                {"units": "hours since 1970-01-01 00:00:00", "standard_name": "time"},
+            )
+        },
+    )
+
+    # Test with missing data
+    with pytest.raises(ValueError):
+        tb_utils.get_spacings(in_xr)
diff --git a/tobac/tests/test_utils_bulk_statistics.py b/tobac/tests/test_utils_bulk_statistics.py
index 4967a779..1db036b0 100644
--- a/tobac/tests/test_utils_bulk_statistics.py
+++ b/tobac/tests/test_utils_bulk_statistics.py
@@ -8,7 +8,10 @@
 import tobac.testing as tb_test
 
 
-def test_bulk_statistics():
+@pytest.mark.parametrize(
+    "id_column, index", [("feature", [1]), ("feature_id", [1]), ("cell", [1])]
+)
+def test_bulk_statistics(id_column, index):
     """
     Test to assure that bulk statistics for identified features are computed as expected.
 
@@ -46,10 +49,10 @@ def test_bulk_statistics():
     )
 
     #### 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
+        out_df, out_seg_mask, test_data_iris, statistic=stats, id_column=id_column
     )
     assert out_segmentation.equals(out_df)
 
@@ -86,11 +89,12 @@ def test_bulk_statistics():
     )
 
     ##### 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
+        out_df, out_seg_mask, test_data_iris, statistic=stats, id_column=id_column
     )
+
     assert out_segmentation.equals(out_df)
 
     # assure that column names in new dataframe correspond to keys in statistics dictionary