doc/bulk_statistics/notebooks/compute_statistics_postprocessing_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tobac example: Compute bulk statistics as a postprocessing step\n",
    "=== \n",
    "\n",
    "Instead of during the feature detection or segmentation process, you can also calculate bulk statistics of your detected/tracked objects as a postprocessing step, i.e. based on a segmentation mask. This makes it possible to combine different datasets and derive statistics for your detected features based on other input fields (e.g., get precipitation statistics under cloud features that were segmented based on brightness temperatures or outgoing longwave radiation). \n",
    "\n",
    "This notebook shows an example for how to compute bulk statistics for detected features as a postprocessing step, that is based on the segmentation mask that we have already created. We perform the feature detection and segmentation with data from [our example for precipitation tracking](https://github.com/tobac-project/tobac/blob/main/examples/Example_Precip_Tracking/Example_Precip_Tracking.ipynb). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:27.087494Z",
     "iopub.status.busy": "2024-03-22T17:17:27.087185Z",
     "iopub.status.idle": "2024-03-22T17:17:31.199088Z",
     "shell.execute_reply": "2024-03-22T17:17:31.194250Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import iris\n",
    "import numpy as np\n",
    "import pandas as pd\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",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:31.211608Z",
     "iopub.status.busy": "2024-03-22T17:17:31.210964Z",
     "iopub.status.idle": "2024-03-22T17:17:36.419147Z",
     "shell.execute_reply": "2024-03-22T17:17:36.416362Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using tobac version 1.5.3\n"
     ]
    }
   ],
   "source": [
    "# Import tobac itself\n",
    "import tobac\n",
    "print('using tobac version', str(tobac.__version__))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.455941Z",
     "iopub.status.busy": "2024-03-22T17:17:36.455446Z",
     "iopub.status.idle": "2024-03-22T17:17:36.460672Z",
     "shell.execute_reply": "2024-03-22T17:17:36.460009Z"
    }
   },
   "outputs": [],
   "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Feature detection** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.464189Z",
     "iopub.status.busy": "2024-03-22T17:17:36.463883Z",
     "iopub.status.idle": "2024-03-22T17:17:36.468048Z",
     "shell.execute_reply": "2024-03-22T17:17:36.467231Z"
    }
   },
   "outputs": [],
   "source": [
    "#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",
    "if not plot_dir.is_dir():\n",
    "    plot_dir.mkdir()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.471915Z",
     "iopub.status.busy": "2024-03-22T17:17:36.471614Z",
     "iopub.status.idle": "2024-03-22T17:17:36.479102Z",
     "shell.execute_reply": "2024-03-22T17:17:36.478486Z"
    }
   },
   "outputs": [],
   "source": [
    "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",
    "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",
    "    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'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.482460Z",
     "iopub.status.busy": "2024-03-22T17:17:36.482198Z",
     "iopub.status.idle": "2024-03-22T17:17:36.706441Z",
     "shell.execute_reply": "2024-03-22T17:17:36.705850Z"
    }
   },
   "outputs": [],
   "source": [
    "Precip=iris.load_cube(str(data_file[0]),'surface_precipitation_average')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.709868Z",
     "iopub.status.busy": "2024-03-22T17:17:36.709531Z",
     "iopub.status.idle": "2024-03-22T17:17:36.716902Z",
     "shell.execute_reply": "2024-03-22T17:17:36.715785Z"
    }
   },
   "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",
    "\n",
    "# get temporal and spation resolution of the data\n",
    "dxy,dt=tobac.get_spacings(Precip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:36.720988Z",
     "iopub.status.busy": "2024-03-22T17:17:36.720671Z",
     "iopub.status.idle": "2024-03-22T17:17:42.032429Z",
     "shell.execute_reply": "2024-03-22T17:17:42.031730Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting feature detection based on multiple thresholds\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "feature detection done\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features saved\n"
     ]
    }
   ],
   "source": [
    "# Feature detection based on based on surface precipitation field and a range of thresholds\n",
    "print('starting feature detection based on multiple thresholds')\n",
    "Features= tobac.feature_detection_multithreshold(Precip,dxy,**parameters_features) \n",
    "print('feature detection done')\n",
    "Features.to_hdf(savedir / 'Features.h5','table')\n",
    "print('features saved')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Segmentation** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:42.036502Z",
     "iopub.status.busy": "2024-03-22T17:17:42.035726Z",
     "iopub.status.idle": "2024-03-22T17:17:42.045228Z",
     "shell.execute_reply": "2024-03-22T17:17:42.042712Z"
    }
   },
   "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",
    "\n",
    "# get temporal and spation resolution of the data\n",
    "dxy,dt=tobac.get_spacings(Precip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:42.049101Z",
     "iopub.status.busy": "2024-03-22T17:17:42.048656Z",
     "iopub.status.idle": "2024-03-22T17:17:52.522098Z",
     "shell.execute_reply": "2024-03-22T17:17:52.514694Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting segmentation based on surface precipitation\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "segmentation based on surface precipitation performed, start saving results to files\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "segmentation surface precipitation performed and saved\n"
     ]
    }
   ],
   "source": [
    "# Perform Segmentation and save resulting mask to NetCDF file:\n",
    "print('Starting segmentation based on surface precipitation')\n",
    "Mask_Precip,Features_Precip=tobac.segmentation_2D(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')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get bulk statistics from segmentation mask file**\n",
    "\n",
    "You can decide which statistics to calculate by providing a dictionary with the name of the metric as keys (this will be the name of the column added to the dataframe) and functions as values. Note that it is also possible to provide input parameter to these functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:52.528489Z",
     "iopub.status.busy": "2024-03-22T17:17:52.527447Z",
     "iopub.status.idle": "2024-03-22T17:17:52.538706Z",
     "shell.execute_reply": "2024-03-22T17:17:52.535020Z"
    }
   },
   "outputs": [],
   "source": [
    "from tobac.utils import get_statistics_from_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-10T23:49:21.717102Z",
     "iopub.status.busy": "2023-07-10T23:49:21.716272Z",
     "iopub.status.idle": "2023-07-10T23:49:21.985540Z",
     "shell.execute_reply": "2023-07-10T23:49:21.984704Z"
    }
   },
   "source": [
    "#### Defining the dictionary for the statistics to be calculated "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:52.548583Z",
     "iopub.status.busy": "2024-03-22T17:17:52.548132Z",
     "iopub.status.idle": "2024-03-22T17:17:52.554124Z",
     "shell.execute_reply": "2024-03-22T17:17:52.553158Z"
    }
   },
   "outputs": [],
   "source": [
    "statistics = {}\n",
    "statistics['mean_precip'] = np.mean\n",
    "statistics['total_precip'] = np.sum\n",
    "statistics['max_precip'] = np.max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some functions, we need to provide additional input parameters, e.g. [np.percentile()](https://numpy.org/doc/stable/reference/generated/numpy.percentile.html). These can be provided as key word arguments in form of a dictionary. So instead of the function, you can provide a tuple with both the function and its respective input parameters: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:52.562198Z",
     "iopub.status.busy": "2024-03-22T17:17:52.561074Z",
     "iopub.status.idle": "2024-03-22T17:17:52.571050Z",
     "shell.execute_reply": "2024-03-22T17:17:52.569254Z"
    }
   },
   "outputs": [],
   "source": [
    "statistics['percentiles'] = (np.percentile, {'q': [95,99]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:52.575766Z",
     "iopub.status.busy": "2024-03-22T17:17:52.575276Z",
     "iopub.status.idle": "2024-03-22T17:17:57.003737Z",
     "shell.execute_reply": "2024-03-22T17:17:57.001134Z"
    }
   },
   "outputs": [],
   "source": [
    "features_with_stats = get_statistics_from_mask(Features_Precip, Mask_Precip, Precip, statistic=statistics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at the output: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:57.007064Z",
     "iopub.status.busy": "2024-03-22T17:17:57.006774Z",
     "iopub.status.idle": "2024-03-22T17:17:57.013033Z",
     "shell.execute_reply": "2024-03-22T17:17:57.012425Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1.629695\n",
       "1    1.409547\n",
       "2    2.441526\n",
       "3    1.938501\n",
       "4    2.486886\n",
       "Name: mean_precip, dtype: object"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features_with_stats.mean_precip.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:57.015793Z",
     "iopub.status.busy": "2024-03-22T17:17:57.015558Z",
     "iopub.status.idle": "2024-03-22T17:17:57.020851Z",
     "shell.execute_reply": "2024-03-22T17:17:57.020126Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    16.296951\n",
       "1    14.095468\n",
       "2    26.856783\n",
       "3    36.831512\n",
       "4    49.737709\n",
       "Name: total_precip, dtype: object"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features_with_stats.total_precip.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-22T17:17:57.023633Z",
     "iopub.status.busy": "2024-03-22T17:17:57.023395Z",
     "iopub.status.idle": "2024-03-22T17:17:57.028693Z",
     "shell.execute_reply": "2024-03-22T17:17:57.028166Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0      ([2.221776068210602, 2.276183712482452],)\n",
       "1    ([1.8030404090881347, 1.8164567756652832],)\n",
       "2      ([3.710712432861328, 3.759503173828125],)\n",
       "3      ([3.940941762924194, 4.042321195602417],)\n",
       "4     ([4.087516045570374, 4.3222578477859495],)\n",
       "Name: percentiles, dtype: object"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features_with_stats.percentiles.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}