{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. Creating a Dataset from numpy arrays\n",
    "The goal of this notebook is to demonstrate how to create a `xarray.Dataset` from scratch. This can be useful if your device is not supported or if you would like to integrate the `dtscalibration` library in your current routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:05.700886Z",
     "iopub.status.busy": "2022-04-06T08:11:05.700153Z",
     "iopub.status.idle": "2022-04-06T08:11:07.170523Z",
     "shell.execute_reply": "2022-04-06T08:11:07.169990Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "\n",
    "from dtscalibration import read_silixa_files\n",
    "\n",
    "# The following line introduces the .dts accessor for xarray datasets\n",
    "import dtscalibration  # noqa: E401  # noqa: E401\n",
    "from dtscalibration.variance_stokes import variance_stokes_constant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a `xarray.Dataset` object, a few things are needed:\n",
    "\n",
    "- timestamps\n",
    "\n",
    "- Stokes signal\n",
    "\n",
    "- anti-Stokes signal\n",
    "\n",
    "- x (length along fiber)\n",
    "\n",
    "Let's grab the data from an existing silixa dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.173171Z",
     "iopub.status.busy": "2022-04-06T08:11:07.172928Z",
     "iopub.status.idle": "2022-04-06T08:11:07.369257Z",
     "shell.execute_reply": "2022-04-06T08:11:07.368742Z"
    }
   },
   "outputs": [],
   "source": [
    "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"single_ended\")\n",
    "\n",
    "ds_silixa = read_silixa_files(directory=filepath, silent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will get all the numpy arrays from this `xarray.Dataset` to create a new one from 'scratch'.\n",
    "\n",
    "Let's start with the most basic data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.394645Z",
     "iopub.status.busy": "2022-04-06T08:11:07.394464Z",
     "iopub.status.idle": "2022-04-06T08:11:07.398245Z",
     "shell.execute_reply": "2022-04-06T08:11:07.397718Z"
    }
   },
   "outputs": [],
   "source": [
    "x = ds_silixa.x.values\n",
    "time = ds_silixa.time.values\n",
    "ST = ds_silixa.st.values\n",
    "AST = ds_silixa.ast.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this data has to be inserted into an xarray `Dataset`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.400661Z",
     "iopub.status.busy": "2022-04-06T08:11:07.400499Z",
     "iopub.status.idle": "2022-04-06T08:11:07.407067Z",
     "shell.execute_reply": "2022-04-06T08:11:07.406546Z"
    }
   },
   "outputs": [],
   "source": [
    "ds = xr.Dataset()\n",
    "ds[\"x\"] = (\"x\", x)\n",
    "ds[\"time\"] = (\"time\", time)\n",
    "ds[\"st\"] = ([\"x\", \"time\"], ST)\n",
    "ds[\"ast\"] = ([\"x\", \"time\"], AST)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.409666Z",
     "iopub.status.busy": "2022-04-06T08:11:07.409313Z",
     "iopub.status.idle": "2022-04-06T08:11:07.417393Z",
     "shell.execute_reply": "2022-04-06T08:11:07.416796Z"
    }
   },
   "outputs": [],
   "source": [
    "print(ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For calibration, a few more paramaters are needed:\n",
    "\n",
    "- acquisition time (for calculating residuals for WLS calibration)\n",
    "\n",
    "- reference temperatures\n",
    "\n",
    "- a double ended flag\n",
    "\n",
    "We'll put these into the custom `xarray.Dataset`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.419863Z",
     "iopub.status.busy": "2022-04-06T08:11:07.419639Z",
     "iopub.status.idle": "2022-04-06T08:11:07.427490Z",
     "shell.execute_reply": "2022-04-06T08:11:07.426850Z"
    }
   },
   "outputs": [],
   "source": [
    "ds[\"acquisitiontimeFW\"] = ds_silixa[\"acquisitiontimeFW\"].values\n",
    "ds[\"userAcquisitionTimeFW\"] = ds_silixa[\"acquisitiontimeFW\"].values\n",
    "ds[\"temp1\"] = ds_silixa[\"probe1Temperature\"]\n",
    "ds[\"temp2\"] = ds_silixa[\"probe2Temperature\"]\n",
    "\n",
    "ds.attrs[\"isDoubleEnded\"] = \"0\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can calibrate the data as usual (ordinary least squares in this example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T08:11:07.430326Z",
     "iopub.status.busy": "2022-04-06T08:11:07.429752Z",
     "iopub.status.idle": "2022-04-06T08:11:07.508805Z",
     "shell.execute_reply": "2022-04-06T08:11:07.508335Z"
    }
   },
   "outputs": [],
   "source": [
    "ds = ds.sel(x=slice(-30, 101))\n",
    "sections = {\n",
    "    \"temp1\": [slice(20, 25.5)],  # warm bath\n",
    "    \"temp2\": [slice(5.5, 15.5)],  # cold bath\n",
    "}\n",
    "\n",
    "st_var, resid = variance_stokes_constant(\n",
    "    ds.dts.st, sections, ds.dts.acquisitiontime_fw, reshape_residuals=True\n",
    ")\n",
    "ast_var, _ = variance_stokes_constant(\n",
    "    ds.dts.ast, sections, ds.dts.acquisitiontime_fw, reshape_residuals=False\n",
    ")\n",
    "out = ds.dts.calibrate_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)\n",
    "out.isel(time=0).tmpf.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}