4-feedzones.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "vN6IjcnKtKEqn2LtN923",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# Geothermal Well Test Analysis with Python\n",
    "### Notebook 4: Fluid velocity analysis and feed zone interpretation\n",
    "#### Irene Wallis and Katie McLean \n",
    "#### Software Underground, Transform 2021\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "BJqasOgylrn8J4vcWIx0",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### Google Colab Setup\n",
    "\n",
    "If you are using Google Colab to run this notebook, we assume you have already followed the Google Colab setup steps outlined [here](https://github.com/ICWallis/T21-Tutorial-WellTestAnalysis).\n",
    "\n",
    "Because we are importing data, we need to \"mount your Google Drive\", which is where we tell this notebook to look for the data files. You will need to mount the Google Drive into each notebook.  \n",
    "\n",
    "1. Run the cell below if you are in Google Colab. If you are not in Google Colab, running the cell below will just return an error that says \"No module named 'google'\". If you get a Google Colab error that says \"Unrecognised runtime 'geothrm'; defaulting to 'python3' Notebook settings\", just ignore it. \n",
    "\n",
    "2. Follow the link generated by running this code. That link will ask you to sign in to your google account (use the one where you have saved these tutorial materials in) and to allow this notebook access to your google drive. \n",
    "\n",
    "3. Completing step 2 above will generate a code. Copy this code, paste below where it says \"Enter your authorization code:\", and press ENTER. \n",
    "\n",
    "Congratulations, this notebook can now import data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "EGWAperoLhU00srKqHaS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "from google.colab import drive\n",
    "drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "YjpLIefYNQnR4vFsJXuc",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "\n",
    "***\n",
    "\n",
    "## 15. Feed zone interpretation in geothermal wells\n",
    "\n",
    "Geothermal wells typically have very long sections of perforated liner (often longer than 1km). This entire length is potentially open to fluid flow from (or into) the reservoir, though whether there is actually any flow at a particular depth depends on whether there is permeability in the reservoir at that depth. In typical geothermal wells in NZ (and elsewhere) there are multiple distinct depths at which there is permeability - called \"permeable feed zones\", or just \"feed zones\". See the figure below for a schematic of a geothermal well showing feed zones. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "u36ZPj2D8XTMN5Gbn8ot",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image('https://raw.githubusercontent.com/ICWallis/T21-Tutorial-WellTestAnalysis/main/Figures/Figure2.jpg',width = 500,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "IdEvIP3LRkQrwtTwCYI7",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "It is of interest to know the depths and relative sizes of the feed zones. This allows the geothermal reservoir engineer to do various things, including: \n",
    "- Correlate feed zones to geological formations to improve future well targeting.\n",
    "- Accurately model the well and its likely future performance.\n",
    "- Maintain the well into the future, for example if targeted well stimulations are required.\n",
    "\n",
    "Feed zones are interpreted from various data which can be collected by PTS (pressure-temperature-spinner) wireline tools. This typically includes: \n",
    "- Temperature profiles during injection, at different flow rates.\n",
    "- Spinner profiles during injection, used to calculate fluid velocity profiles for different flow rates.\n",
    "- Pressure and temperature profiles during progressive heat-up after injection stops. \n",
    "- Fracture datasets from borehole imaging (if available). \n",
    "\n",
    "Feed zones are initially interpreted from the data types above, which are captured during completion testing and heat-up. Later during output testing, PTS data are captured as the well is producing, and this data confirms which of the feed zones are active under those conditions. As more testing is done and more data is available, the feed zone interpretation can evolve. Also the feed zones in a well can change over time, particularly if there is scaling blocking up particular feed zones. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "59NHIngKgVAST7OKgWJF",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 16. Import, munge and check data\n",
    "\n",
    "## 16.1 Use bespoke functions to import and munge data\n",
    "\n",
    "Install packages required for this notebook. If you do not already have iapws in your environment, then you will need to pip install it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "R66U43TkvXp7jmZAJulw",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "!pip install iapws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "2mzWoPijFfgQL819jIz6",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "import iapws # steam tables\n",
    "import openpyxl\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "from IPython.display import Image\n",
    "from ipywidgets import interactive, Layout, FloatSlider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "QG6Px2ZVLv48U5Y5Ix4t",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "def timedelta_seconds(dataframe_col, test_start):\n",
    "    '''\n",
    "    Make a float in seconds since the start of the test\n",
    "\n",
    "    args:   dataframe_col: dataframe column containing datetime objects\n",
    "            test_start: test start time formatted '2020-12-11 09:00:00'\n",
    "\n",
    "    returns: float in seconds since the start of the test\n",
    "    '''\n",
    "    test_start_datetime = pd.to_datetime(test_start)\n",
    "    list = []\n",
    "    for datetime in dataframe_col:\n",
    "        time_delta = datetime - test_start_datetime\n",
    "        seconds = time_delta.total_seconds()\n",
    "        list.append(seconds)\n",
    "    return list\n",
    "\n",
    "\n",
    "\n",
    "def read_flowrate(filename):\n",
    "    ''' \n",
    "    Read PTS-2-injection-rate.xlsx in as a pandas dataframe and munge for analysis\n",
    "\n",
    "    args: filename is r'PTS-2-injection-rate.xlsx'\n",
    "\n",
    "    returns: pandas dataframe with local NZ datetime and flowrate in t/hr\n",
    "    '''\n",
    "    df = pd.read_excel(filename, header=1) \n",
    "    df.columns = ['raw_datetime','flow_Lpm']\n",
    "\n",
    "    list = []\n",
    "    for date in df['raw_datetime']:\n",
    "        newdate = datetime.fromisoformat(date)\n",
    "        list.append(newdate)\n",
    "    df['ISO_datetime'] = list \n",
    "\n",
    "    list = []\n",
    "    for date in df.ISO_datetime:\n",
    "        newdate = pd.to_datetime(datetime.strftime(date,'%Y-%m-%d %H:%M:%S'))\n",
    "        list.append(newdate)\n",
    "    df['datetime'] = list\n",
    "\n",
    "    df['flow_tph'] = df.flow_Lpm * 0.060\n",
    "\n",
    "    df['timedelta_sec'] = timedelta_seconds(df.datetime, '2020-12-11 09:26:44.448')\n",
    "\n",
    "    df.drop(columns = ['raw_datetime', 'flow_Lpm', 'ISO_datetime'], inplace = True)\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "\n",
    "def read_pts(filename):\n",
    "    '''\n",
    "    Read PTS-2.xlsx in as a Pandas dataframe and munge for analysis\n",
    "\n",
    "    args: filename is r'PTS-2.xlsx'\n",
    "\n",
    "    returns: Pandas dataframe with datetime (local) and key coloumns of PTS data with the correct dtype\n",
    "    '''\n",
    "    df = pd.read_excel(filename)\n",
    "\n",
    "    dict = {\n",
    "        'DEPTH':'depth_m',\n",
    "        'SPEED': 'speed_mps',\n",
    "        'Cable Weight': 'cweight_kg',\n",
    "        'WHP': 'whp_barg',\n",
    "        'Temperature': 'temp_degC',\n",
    "        'Pressure': 'pressure_bara',\n",
    "        'Frequency': 'frequency_hz'\n",
    "    }\n",
    "    df.rename(columns=dict, inplace=True)\n",
    "\n",
    "    df.drop(0, inplace=True)\n",
    "    df.reset_index(drop=True, inplace=True)\n",
    "\n",
    "    list = []\n",
    "    for date in df.Timestamp:\n",
    "        newdate = openpyxl.utils.datetime.from_excel(date)\n",
    "        list.append(newdate)\n",
    "    df['datetime'] = list\n",
    "\n",
    "    df.drop(columns = ['Date', 'Time', 'Timestamp','Reed 0',\n",
    "       'Reed 1', 'Reed 2', 'Reed 3', 'Battery Voltage', \n",
    "       'PRT Ref Voltage','SGS Voltage', 'Internal Temp 1', \n",
    "       'Internal Temp 2', 'Internal Temp 3','Cal Temp', \n",
    "       'Error Code 1', 'Error Code 2', 'Error Code 3',\n",
    "       'Records Saved', 'Bad Pages',], inplace = True)\n",
    "    \n",
    "    df[\n",
    "        ['depth_m', 'speed_mps','cweight_kg','whp_barg','temp_degC','pressure_bara','frequency_hz']\n",
    "    ] = df[\n",
    "        ['depth_m','speed_mps','cweight_kg','whp_barg','temp_degC','pressure_bara','frequency_hz']\n",
    "        ].apply(pd.to_numeric)\n",
    "    \n",
    "    df['timedelta_sec'] = timedelta_seconds(df.datetime, '2020-12-11 09:26:44.448')\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "\n",
    "def append_flowrate_to_pts(flowrate_df, pts_df):\n",
    "    '''\n",
    "    Add surface flowrate to pts data\n",
    "\n",
    "    Note that the flowrate data is recorded at a courser time resolution than the pts data\n",
    "    The function makes a linear interpolation to fill the data gaps\n",
    "    Refer to bonus-combine-data.ipynb to review this method and adapt it for your own data\n",
    "\n",
    "    Args:   flowrate and pts dataframes generated by the read_flowrate and read_pts functions\n",
    "\n",
    "    Returns: pts dataframe with flowrate tph added\n",
    "    \n",
    "    '''\n",
    "    flowrate_df = flowrate_df.set_index('timedelta_sec')\n",
    "    pts_df = pts_df.set_index('timedelta_sec')\n",
    "    combined_df = pts_df.join(flowrate_df, how = 'outer',  lsuffix = '_pts', rsuffix = '_fr')\n",
    "    combined_df.drop(columns = ['datetime_fr'], inplace = True)\n",
    "    combined_df.columns = ['depth_m', 'speed_mps', 'cweight_kg', 'whp_barg', 'temp_degC',\n",
    "       'pressure_bara', 'frequency_hz', 'datetime', 'flow_tph']\n",
    "    combined_df['interpolated_flow_tph'] = combined_df['flow_tph'].interpolate(method='linear')\n",
    "    trimmed_df = combined_df[combined_df['depth_m'].notna()]\n",
    "    trimmed_df.reset_index(inplace=True)\n",
    "    return trimmed_df\n",
    "\n",
    "\n",
    "\n",
    "def find_index(value, df, colname):\n",
    "    '''\n",
    "    Find the dataframe index for the exact matching value or nearest two values\n",
    "\n",
    "    args:   value: (float or int) the search term\n",
    "            df: (obj) the name of the dataframe that is searched\n",
    "            colname: (str) the name of the coloum this is searched\n",
    "\n",
    "    returns:  dataframe index(s) for the matching value or the two adjacent values\n",
    "              rows can be called from a df using df.iloc[[index_number,index_number]]\n",
    "    '''\n",
    "    exactmatch = df[df[colname] == value]\n",
    "    if not exactmatch.empty:\n",
    "        return exactmatch.index\n",
    "    else:\n",
    "        lowerneighbour_index = df[df[colname] < value][colname].idxmax()\n",
    "        upperneighbour_index = df[df[colname] > value][colname].idxmin()\n",
    "        return [lowerneighbour_index, upperneighbour_index] \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "dMcLxlA0SCBMof4ipPar",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "The cells below will take a little while to run because it includes all steps required to import and munge the data (i.e., everything we did in notebook 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "lC8kOpwoXfq75gKlLUcz",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "flowrate = read_flowrate(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-FlowRate.xlsx')\n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#flowrate = read_flowrate(r'Data-FlowRate.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "uQB8GyMHdlfHv4pJM60z",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "pts = read_pts(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-PTS.xlsx')\n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#pts = read_pts(r'Data-PTS.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "JUbFD9HyCnsHG4u4EkWB",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "9yETDgOAswJn78hIt1Fd",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "dt9susHmwiZLaCJofcDQ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 16.2 Add flowrate values to the pts dataframe\n",
    "\n",
    "Our fluid velocity analysis requires that we know the pump rate and spinner frequency. There are several ways we could approach this:\n",
    "\n",
    "1. We could assume that the pump rate was held perfectly steady at the planned pump rate and set a single value\n",
    "2. We could use the actual flowrate data if that is available\n",
    "\n",
    "As we have good quality pump data, we use the bespoke function below to append the flowrate data to the pts dataframe. As the flowrate data is recorded at a coarser time resolution than the pts data, we used linear interpolation to fill the gaps. \n",
    "\n",
    "If you are using this workflow on your own data, you need to adjust the column names in the function. This method is documented in bonus-combine-date.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "BvsJGUKSq92P7DRfOfp4",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts = append_flowrate_to_pts(flowrate, pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "cF8u3S3WohPSUUpVCO2k",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "BYfP0CoOlAymq3VaeNDh",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "oJ3ypu37pFyUyYhHD5cz",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 16.3 Check the data\n",
    "\n",
    "It is good practice to check your data after import. \n",
    "\n",
    "You can use the Pandas methods listed in Section 2.1.1 (1-intro-and-data.ipynb) to check your data and the plots below. \n",
    "\n",
    "### 16.3.1 Visualise spinner by depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "fy1zgJ6MDpKSPm5NZKEL",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(24,8),sharey=True)\n",
    "\n",
    "ax1.scatter(pts.frequency_hz, pts.depth_m, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "ax2.scatter(pts.datetime, pts.depth_m, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax3 = ax2.twinx()\n",
    "ax3.plot(flowrate.datetime, flowrate.flow_tph, \n",
    "    c='k', linestyle = '-', linewidth = 3, alpha = 0.3, \n",
    "    label='Surface pump flowrate')\n",
    "\n",
    "ax1.set_ylim(1000,0)\n",
    "ax1.set_xlim(-30,30)\n",
    "\n",
    "ax1.set_ylabel('Depth [m]')\n",
    "ax1.set_xlabel('Spinner frequency [hz]')\n",
    "\n",
    "ax2.set_xlabel('Time [hh:mm]')\n",
    "ax2.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))\n",
    "\n",
    "ax3.set_ylabel('Flowrate [t/hr]')\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "oPC0zYGPilLa0Ut1hTSQ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 16.3.2 Visualise spinner by time\n",
    "\n",
    "When we plot spinner frequency by time we can see the sequence of up and down runs of the tool inside the well. A steady tool speed is maintained within each of these runs. However, as the tool approaches the bottom and top of the logged interval, it slows down before it stops. As we zoom into the data by changing the time interval plotted, we see where the tool is slowing before it stops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "U9ZHQzzejuyebvUB59lz",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(21,6))\n",
    "\n",
    "ax1.scatter(pts.datetime, pts.frequency_hz, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "ax2.scatter(pts.datetime, pts.depth_m, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax3 = ax2.twinx()\n",
    "ax3.plot(flowrate.datetime, flowrate.flow_tph, \n",
    "    c='k', linestyle = '-', linewidth = 1, alpha = 0.3, \n",
    "    label='Surface pump flowrate')\n",
    "\n",
    "ax4 = ax1.twinx()\n",
    "ax4.plot(pts.datetime, pts.depth_m, \n",
    "    c='k', linestyle = '-', linewidth = 1, alpha = 0.3, # edit linewidth to make visible\n",
    "    label='Tool depth [m]')\n",
    "\n",
    "ax4.set_ylim(1000,-1000)\n",
    "ax4.set_ylabel('Tool depth [m]')\n",
    "\n",
    "ax1.set_ylim(-30,30)\n",
    "ax1.set_ylabel('Spinner frequency [hz]')\n",
    "\n",
    "ax2.set_ylim(1000,0)\n",
    "ax2.set_ylabel('Tool depth [m]')\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.set_xlabel('Time [hh:mm]')\n",
    "    ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))\n",
    "\n",
    "ax3.set_ylabel('Flowrate [t/hr]')\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    ax.grid()\n",
    "\n",
    "ax4.set_ylim(1000,400)\n",
    "\n",
    "\n",
    "# Uncomment the code below\n",
    "# Edit the times to limit the plot to the desired time period\n",
    "\n",
    "#start_time = pd.to_datetime('2020-12-11 09:30:00')\n",
    "#end_time = pd.to_datetime('2020-12-11 10:30:00')\n",
    "\n",
    "#ax1.set_xlim(start_time,end_time)\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "e5OPJfx6hoUdmwJSENy4",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 16.4 Clean data\n",
    "\n",
    "We will remove data acquired when the tool is stationary or slowing and the data acquired while in the cased sections. \n",
    "\n",
    "To understand the method used here and how we decided which data to filter, refer to bonus-filter-by-toolspeed.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "nOlEWMCFXpQKTlhgF7Al",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "moving_pts = pts[\n",
    "    (pts.speed_mps > 0.9 ) & (pts.speed_mps < pts.speed_mps.max()) | \n",
    "    (pts.speed_mps > pts.speed_mps.min() ) & (pts.speed_mps < -0.9)\n",
    "    ]\n",
    "\n",
    "production_shoe = 462.5  \n",
    "\n",
    "clean_pts = moving_pts[(moving_pts.depth_m < moving_pts.depth_m.max()) & (moving_pts.depth_m > production_shoe)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "PaW2aLkQnYjle3rUSFAM",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "We now have a new working dataframe called clean_pts that will be used in the analysis. \n",
    "\n",
    "If we repeat the same plot we made in Section 14.3.2, we can see which data has been removed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "lWLdJW6nQHM9psaUEVy2",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(28,8))\n",
    "\n",
    "ax1.scatter(clean_pts.datetime, clean_pts.frequency_hz, c = clean_pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "ax2.scatter(clean_pts.datetime, clean_pts.depth_m, c = clean_pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax3 = ax2.twinx()\n",
    "ax3.plot(flowrate.datetime, flowrate.flow_tph, \n",
    "    c='k', linestyle = '-', linewidth = 3, alpha = 0.3, \n",
    "    label='Surface pump flowrate')\n",
    "\n",
    "ax4 = ax1.twinx()\n",
    "ax4.plot(pts.datetime, pts.depth_m, \n",
    "    c='k', linestyle = '-', linewidth = 0, alpha = 0.3,  # edit linewidth to make visible\n",
    "    label='Tool depth [m]')\n",
    "ax4.set_ylim(1000,400)\n",
    "ax4.set_ylabel('Tool depth [m]')\n",
    "\n",
    "ax1.set_ylim(-30,30)\n",
    "ax1.set_ylabel('Spinner frequency [hz]')\n",
    "\n",
    "ax2.set_ylim(1000,0)\n",
    "ax2.set_ylabel('Tool depth [m]')\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.set_xlabel('Time [hh:mm]')\n",
    "    ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "\n",
    "ax3.set_ylabel('Flowrate [t/hr]')\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    ax.grid()\n",
    "\n",
    "# Uncomment the code below\n",
    "# Edit the times to limit the plot to the desired time period\n",
    "\n",
    "#start_time = pd.to_datetime('2020-12-11 14:30:00')\n",
    "#end_time = pd.to_datetime('2020-12-11 15:30:00')\n",
    "\n",
    "#ax1.set_xlim(start_time,end_time)\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "oO2kUyD16T1XyjgoaQCq",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "## 17. Select data by flow rate\n",
    "\n",
    "The cross-plot analysis is done using PTS tool passes conducted at a single surface pump flow rate. \n",
    "\n",
    "In this section we generate a dataframe of PTS data for each flow rate using an interactive plotting tool.\n",
    "\n",
    "### 17.1 Interactive plot\n",
    "\n",
    "Use the sliders on the interactive plot below to find a time before (start, green) and after (end, red) each set of PTS tool passes that were conducted at a single flow rate. Because we have removed the stationary and slowing tool data, the slider values need to be close to the start and end of our desired data interval, but they do not need to be exact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "JM2p5nsUYV9hMsIIIOQy",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "min_timestamp = pts.timedelta_sec.iloc[0]\n",
    "max_timestamp = pts.timedelta_sec.iloc[-1]\n",
    "\n",
    "def subselect_plot(start_value, stop_value):\n",
    "    f,ax = plt.subplots(1,1, figsize = (20,6))\n",
    "    ax.scatter(clean_pts.timedelta_sec, clean_pts.depth_m,\n",
    "        c = 'k', s = 1, linewidths = 0, label = 'Tool depth')\n",
    "    ax1 = ax.twinx()\n",
    "    ax1.plot(flowrate.timedelta_sec, flowrate.flow_tph, \n",
    "        ':', c='k', label='Surface pump flowrate')\n",
    "    ymin = pts.depth_m.min()\n",
    "    ymax = pts.depth_m.max() + 100\n",
    "    ax.vlines(start_value, ymin, ymax, color='tab:green')\n",
    "    ax.vlines(stop_value, ymin, ymax, color='tab:red')\n",
    "    ax.set_ylim(pts.depth_m.max() + 100, 0)\n",
    "    ax.set_xlabel('Time elapsed since the test started [sec]')\n",
    "    ax.set_ylabel('Tool depth [m]')\n",
    "    ax1.set_ylabel('Flowrate [t/hr]')\n",
    "\n",
    "result = interactive(subselect_plot,\n",
    "         \n",
    "         start_value = FloatSlider\n",
    "         (\n",
    "             value = (max_timestamp - min_timestamp)/3 + min_timestamp,\n",
    "             description = 'start',\n",
    "             min = min_timestamp, \n",
    "             max = max_timestamp, \n",
    "             step = 10, \n",
    "             continuous_update=False,\n",
    "             layout = Layout(width='80%'),\n",
    "             ),\n",
    "          \n",
    "          stop_value = FloatSlider\n",
    "          (\n",
    "             value = (max_timestamp - min_timestamp)/2 + min_timestamp, \n",
    "             description = 'stop',\n",
    "             min = min_timestamp, \n",
    "             max = max_timestamp, \n",
    "             step = 10, \n",
    "             continuous_update=False,\n",
    "             layout = Layout(width='80%')\n",
    "             )\n",
    ")\n",
    "\n",
    "display(result);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "vl1HbMNzVToeAtKKshcK",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "print(\n",
    "    'start =',result.children[0].value, \n",
    "    '\\nstop =', result.children[1].value,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "xwvz0af5gEBj6exnoWpj",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 17.2 Record your analysis\n",
    "\n",
    "We want to make our completion test analysis repeatable and easy to come back to and check. Subsequently, we take the range selected above and manually define objects for making a PTS dataframe with only one flow rate.  \n",
    "\n",
    "We could have defined the start and stop objects used to generate our single rate dataframes as the result.childern\\[0\\].value and result.childern\\[1\\].value objects. But if we did this, you would lose your work because these ipywidget results change every time the sliders are moved or the notebook is re-run. \n",
    "\n",
    "Copy-paste the timestamps printed by the cell above into the markdown cell below to preserve our analysis. We can also take this opportunity to record any metadata that will help others (or our future self) understand our analysis.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "hYLYyNxpD7vUG50asnLC",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "**First pump flow rate (lowest rate)**\n",
    "\n",
    "    Insert your results here\n",
    "\n",
    "**Second pump flow rate (highest rate)**\n",
    "\n",
    "    Insert your results here\n",
    "\n",
    "**Third pump flow rate (middle rate)**\n",
    "\n",
    "    Insert your results here\n",
    "\n",
    "The how the pumps were shut off early in the third rate so we selected data from before they were shut off.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "cT0oL1PN7lWKmWVJHVH1",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 17.3 Make a PTS dataframe for each flow rate\n",
    "\n",
    "Select the data from the clean_pts dataframe for each of the three flow rates using the timestamps generated with the interactive plot\n",
    "\n",
    "#### 17.3.1 First flow rate (lowest pump rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "booNYdR2KKygdYN9jSYi",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# First flow rate\n",
    "\n",
    "start = 240.0\n",
    "stop = 3740.0 \n",
    "\n",
    "pts_first_rate = clean_pts[\n",
    "    (clean_pts.timedelta_sec > start) \n",
    "    & (clean_pts.timedelta_sec < stop)\n",
    "]\n",
    "\n",
    "pts_first_rate.tail(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "PIUK67MQfQu0y0hI3yN3",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### 17.3.2 Second flow rate (highest pump rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "kaxE45LVs2oqLM3b3fzS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Second flowrate\n",
    "\n",
    "start = 9460.0 \n",
    "stop = 12780.0 \n",
    "\n",
    "pts_second_rate = clean_pts[\n",
    "    (clean_pts.timedelta_sec > start) \n",
    "    & (clean_pts.timedelta_sec < stop)\n",
    "]\n",
    "\n",
    "pts_second_rate.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "hTUDtEfyl0bTGUQ5evzP",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### 17.3.3 Third flow rate (middle pump rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "CsZWPyKZw8RtJTbgO2gp",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Third flowrate\n",
    "\n",
    "start = 18660.0 \n",
    "stop = 21020.0 \n",
    "\n",
    "\n",
    "pts_third_rate = clean_pts[\n",
    "    (clean_pts.timedelta_sec > start) \n",
    "    & (clean_pts.timedelta_sec < stop)\n",
    "]\n",
    "\n",
    "pts_third_rate.tail(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "TsKS1P7O2b0Y0F609imO",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### 17.3.4 Plot data from one flow rate\n",
    "\n",
    "The plot below offers us the opportunity to look at the raw data we have just selected. \n",
    "\n",
    "Edit the pts_df and flow_df objects to switch between the three flow rate dataframes we generated above.\n",
    "\n",
    "Note how the flow rate is approximately the same but does vary with time. The right-hand plot nicely shows how the spinner frequency (hertz) various as the tool is run at different speeds up (positive values) and down (negative values) inside the well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "b57RMyBwSVHpchWNtcB2",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# define which dataframe you want to plot\n",
    "pts_df = pts_first_rate\n",
    "\n",
    "# test plot the data\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(20,6))\n",
    "\n",
    "ax1.scatter(pts_df.datetime, pts_df.frequency_hz, \n",
    "    c = pts_df.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax2.scatter(pts_df.frequency_hz, pts_df.depth_m, \n",
    "    c = pts_df.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax1.set_ylabel('Spinner frequency [hz]')\n",
    "ax1.set_xlabel('Time [hh:mm]')\n",
    "ax1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))\n",
    "\n",
    "ax2.set_xlabel('Spinner frequency [hz]')\n",
    "ax2.set_ylabel('Depth [m]')\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "N2ULzy6izNiOMYxjMS4X",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 18. Fluid velocity analysis\n",
    "\n",
    "We find the fluid velocity inside the well at any depth by determining the speed at which the PTS tool has to be travelling to match the speed of the fluid inside the well. When this is true the spinner will not turn and so the frequency will be zero. \n",
    "\n",
    "In summary, the cross-plot method includes the following steps:\n",
    "1. Define the cross-plot interval for analysis (a usual default is all data within 1 meter)\n",
    "2. Select the PTS data inside that interval\n",
    "3. Generate a linear interpolation of frequency (x) and tool speed (y)\n",
    "4. Return the y-intercept, which is the zero spin or fluid velocity\n",
    "5. Return data that helps us to QA/QC and clean the analysis results (R2, data used in the model fit, number of data points)\n",
    "6. QA/QC result\n",
    "7. Clean result to remove suspect values\n",
    "\n",
    "This is done for each of the three pump flow rate and at the end the results are interpreted along with the temperature profiles to identify feed zones.\n",
    "\n",
    "## 18.1 Illustrate the cross-plot method\n",
    "\n",
    "In this section, we use one meter of data to illustrate the cross-plot method implemented in this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "F8nrpItwO8WgGxX9ZZbG",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "interval_top = 700 # shallowest depth in the cross-plot analysis interval\n",
    "interval_bottom = 701 # deepest depth in the cross-plot analysis interval\n",
    "\n",
    "selected_data = pts_second_rate[\n",
    "    (pts_second_rate.depth_m > interval_top ) & (pts_second_rate.depth_m < interval_bottom)\n",
    "    ]\n",
    "\n",
    "selected_data = selected_data[selected_data['frequency_hz'].notna()]\n",
    "selected_data = selected_data[selected_data['pressure_bara'].notna()]\n",
    "\n",
    "linear_model = stats.linregress(selected_data.frequency_hz, selected_data.speed_mps)\n",
    "\n",
    "test_slope = linear_model[0]\n",
    "test_intercept = linear_model[1] # this is the tool speed that matches the fluid velocity\n",
    "test_rvalue = linear_model[2] # this is how well the model fits the data and will be a filter\n",
    "\n",
    "print(test_rvalue, test_rvalue**2, test_intercept) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "o439IsKQPnZItOdsFcA2",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "**Linear regression**\n",
    "\n",
    "We find the speed associated with zero-frequency spinner by using linear interpolation. We plan to add a bi-linear interpolation method to this notebook in future. There are many Python packages that can be used to do a linear regression. We selected the stats.linregress because it is fast and easy to use. \n",
    "\n",
    "The stats.linregress method returns the R value, which is a number between 1 and -1 that describes the relationship between independent (x) and dependant (y) variable:\n",
    "- -1 indicates that an **increase** in x has an associated **decrease** in y\n",
    "- +1 indicates that an **increase** in x has an associated **increase** in y\n",
    "- 0 indicates there is **no relationship** between x and y\n",
    "\n",
    "The $R^2$ value generally tells us how related the two variables are: More specifically, it describes the proportion of variation in the dependant variable (y) that can be predicted from the independent variable (x). $R^2$ varies between 0 and 1 (a percentage). If $R^2 = 0.99$ then the model explains 99% of the variation can be explained. While an $R^2 = 0.3$ indicates that the model explains only 30% of the variation. However, these metrics may not well describe the shape of our data (check out [Anscombe's quartet](https://en.wikipedia.org/wiki/Anscombe%27s_quartet)) and a percentage may not be easy to think about in relation to our data. \n",
    "\n",
    "Root mean squared error would probably be a better metric of fit quality (and is on our to-do list), but because we are working with a small number of data in each interval we can use $R^2$ as a nice, computationally cheap approach to evaluating fit quality. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "4BXdw2DeYrYuCmw8wh5L",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax) = plt.subplots(1, 1,figsize=(6,6))\n",
    "   \n",
    "ax.scatter(selected_data.frequency_hz, selected_data.speed_mps,\n",
    "    color = 'none', edgecolors = 'k', marker='o', s = 50, linewidths = 1, label='Data')\n",
    "\n",
    "model_y_vals = [] # y values using our model\n",
    "for n in selected_data.frequency_hz:\n",
    "    model_y_vals.append(test_slope * n + test_intercept)\n",
    "\n",
    "ax.plot(selected_data.frequency_hz, model_y_vals, \n",
    "    color='tab:orange', linestyle='-', linewidth=3, alpha=0.5, label='Linear fit')\n",
    "\n",
    "ax.scatter(0,test_intercept,\n",
    "    color='tab:orange', s = 100, label='Intercept - 0 spin')\n",
    "\n",
    "ax.hlines(0, -40, 40, color = 'k', linewidth = 0.5)\n",
    "ax.vlines(0, -2, 2, color = 'k', linewidth = 0.5)\n",
    "\n",
    "ax.set_xlim(-40,40)\n",
    "ax.set_ylim(-2,2)\n",
    "\n",
    "ax.set_xlabel('Spinner frequency [hz]')\n",
    "ax.set_ylabel('Tool speed [mps]')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "wnbuZS0iQW8NzVcJZ4H0",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "In the above cross-plot, the y axis intercept from our linear_model is the fluid velocity: it is the point at which the logging tool is moving at the same speed as the fluid inside the well. \n",
    "\n",
    "We will do this many times down the well for a specified depth interval that depends on the resolution and quality of our data. In this method any depth interval can be used, it does not have to be 1 m. A depth interval of 0.5 m has been chosen for the example below due to high data quality.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "jaZoJVoAkQOwBFQjieyJ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 18.2 Functions for cross-plot analysis\n",
    "\n",
    "The cross-plot analysis method in this notebook is a series of functions that are written assuming the dataframe column headers in this tutorial. They are not yet generalised. We have only just started working on this method and plan to refine it with time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "d2ASQ7bPkNV7XrMJXqUg",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# To do: write an error for validation if the bottom is less than the top. Have it return an informative error. \n",
    "\n",
    "def analysis_steps(analysis_top, analysis_bottom, step_size):\n",
    "    '''Make lists that define top and bottom of each analysis interval\n",
    "\n",
    "    The cross_plot_analysis function requires that we pass in two numbers that \n",
    "    define the top and bottom of each data interval that we will do the\n",
    "    cross-plot analysis on. \n",
    "\n",
    "    Args:   analysis_top: The shallowest depth of the spinner analysis interval\n",
    "            analysis_bottom: The deepest depth of the spinner analysis interval\n",
    "            step_size: the length of each cross-plot analysis interval within the spinner analysis interval\n",
    "\n",
    "    Returns:    list_tops: List of top of each cross-plot analysis interval\n",
    "                list_bots: List of the bottom of each cross-plot analysis interval   \n",
    "    \n",
    "    '''\n",
    "\n",
    "    list_tops = np.arange(\n",
    "        start = analysis_top, \n",
    "        stop = analysis_bottom - step_size, \n",
    "        step = step_size)\n",
    "\n",
    "    list_bots = np.arange(\n",
    "        start = analysis_top + step_size, \n",
    "        stop = analysis_bottom, \n",
    "        step = step_size)\n",
    "\n",
    "    return list_tops, list_bots\n",
    "\n",
    "# TODO: generalise this function so it does not matter what the column headers are\n",
    "\n",
    "def cross_plot_analysis(dataframe, interval_top, interval_bottom):  \n",
    "    '''\n",
    "    Cross plot analysis of spinner frequency and tool speed data to find fluid velocity\n",
    "\n",
    "    The method selects from the PTS dataframe between the interval_top and interval_bottom, \n",
    "    tests if there are any data in that interval and then calculates a linear model if there is. \n",
    "    The linear interpolation operates on spinner frequency (x) and tool speed (y) to find the fluid velocity (y-intercept).\n",
    "    The method also returns the model slope, R squared (goodness of fit), and the number of data points used.\n",
    "\n",
    "    Args:           dataframe, \n",
    "                    interval_top, \n",
    "                    interval_bottom\n",
    "\n",
    "    Returns:        freq_data: (list) frequency data used in that cross-plot interval \n",
    "                    speed_data: (list) tool speed data used in that cross-plot interval \n",
    "                    linear_model[1]: liner model intercept which is equivalent to the fluid velocity\n",
    "                    linear_model[0]: linear model slope \n",
    "                    r_squared: goodness of fit  \n",
    "                    number_of_df_rows: number of dataframe rows used in that cross-plot interval\n",
    "\n",
    "    '''  \n",
    "    # select the cross-plot interval\n",
    "    df = dataframe[(dataframe.depth_m > interval_top) & (dataframe.depth_m < interval_bottom)]\n",
    "    # remove nan values\n",
    "    df = df[df['frequency_hz'].notna()]\n",
    "    df = df[df['pressure_bara'].notna()]\n",
    "    # define data values for linear interpolation\n",
    "    freq_data = df.frequency_hz.tolist()\n",
    "    speed_data = df.speed_mps.tolist()\n",
    "    # test if there is any kind of data in that cross-plot interval\n",
    "    number_of_df_rows = df.shape[0]\n",
    "    if number_of_df_rows > 1:\n",
    "        # if there is data, do the linear regression\n",
    "        linear_model = stats.linregress(df.frequency_hz, df.speed_mps)\n",
    "        r_squared = linear_model[2]**2\n",
    "    else: \n",
    "        # if there is no data, return None\n",
    "        linear_model = [None,None]\n",
    "        r_squared = None\n",
    "    return freq_data, speed_data, linear_model[1], linear_model[0], r_squared, number_of_df_rows\n",
    "\n",
    "# TODO: Need to test for the presence of data for the cross-plot intervals rather than the dataframe rows hack\n",
    "# Currently I am assuming that the number of rows in the dataframe returned by the cross_plot_analysis function\n",
    "# is equivalent to the number of values available for the cross-plot analysis. However, this may not be the case.\n",
    "\n",
    "def calc_fluid_velocity(single_flowrate_df, top, bottom, step):\n",
    "    '''\n",
    "    Calculate fluid velocity from PTS data using cross-plot method for each defined step. \n",
    "    \n",
    "    Note that this function has not yet been generalised, so it assumes the column headers generated by tutorial method.\n",
    "\n",
    "    Args:       single_flowrate_df: (Pandas dataframe) generated using the method in this tutoral \n",
    "                top: (int) shallowest data depth \n",
    "                bottom: (int) deepest data depth \n",
    "                step: (int) interval thickness for each cross-plot\n",
    "    \n",
    "    Returns:    Pandas dataframe containing depth, fluid velocity (aka model intercept), \n",
    "                model slope, model R2, data in each cross-plot interval (observation number),\n",
    "                spinner frequency and tool speed data contained in the cross-plot interval.\n",
    "    '''\n",
    "    # define the interval steps\n",
    "    list_tops, list_bots = analysis_steps(top, bottom, step)\n",
    "\n",
    "    # define the lists that results will be placed into during the for loop\n",
    "    depth = [] # half way between top and bottom of step\n",
    "    frequency_data = [] # spinner frequency data in that cross-plot interval\n",
    "    speed_data = [] # tool speed data in that cross-plot interval\n",
    "    fluid_velocity = [] # model intercept\n",
    "    slope = [] # model slope\n",
    "    rsquared = [] # goodness of fit\n",
    "    obs_num = [] # number of observations in the step\n",
    "\n",
    "    # calculate a linear model for each cross-plot interval using a for loop\n",
    "    for top, bot in zip(list_tops, list_bots):\n",
    "        d = (bot - top)/2 + top \n",
    "        depth.append(d)\n",
    "        fd, sd, v, s, r, obs = cross_plot_analysis(single_flowrate_df, top, bot)\n",
    "        frequency_data.append(fd)\n",
    "        speed_data.append(sd)\n",
    "        fluid_velocity.append(v)\n",
    "        slope.append(s)\n",
    "        rsquared.append(r)\n",
    "        obs_num.append(obs) \n",
    "\n",
    "    # turn the results lists into a Pandas dataframe\n",
    "    df = pd.DataFrame()\n",
    "    df['depth_m'] = depth\n",
    "    df['intercept_velocity_mps'] = fluid_velocity\n",
    "    df['slope'] = slope\n",
    "    df['r_squared'] = rsquared\n",
    "    df['obs_num'] = obs_num\n",
    "    df['frequency_data_hz'] = frequency_data\n",
    "    df['speed_data_mps'] = speed_data\n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "IJ97fVdtm3yC5XOCc8lR",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 18.3 Calculate fluid velocity\n",
    "\n",
    "Call the wrapper function calc_fluid_velocity to find the fluid velocity for each of our flow rates. \n",
    "\n",
    "It is informative to trial this method for various depth intervals and see what happens to the $R^2$ and number of values in each cross-plot interval. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "tEpnRe4oQaACJj0vLUuL",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fvelocity_first_rate = calc_fluid_velocity(pts_first_rate, top = 460, bottom = 925, step = 0.5)\n",
    "fvelocity_second_rate = calc_fluid_velocity(pts_second_rate, top = 460, bottom = 925, step = 0.5)\n",
    "fvelocity_third_rate = calc_fluid_velocity(pts_third_rate, top = 460, bottom = 925, step = 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "o2kar8gJ6ZohzrSCdqTg",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fvelocity_first_rate.tail(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "2T86LZuQK0JP29pGTfKv",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fvelocity_second_rate.tail(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "sbLZST6BL8Y6Cqj8mzjG",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fvelocity_third_rate.tail(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "o5EPxlXqPgtFU28RghyC",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 18.4 Visulise raw fluid velocity results\n",
    "\n",
    "Below we plot the results for one flow rate. Edit the object fluid_velocity_df to change the results set you would like to view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "dwsrIAEa4h0ED3zQv5wL",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts_first_rate.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "S5lGxMTbPkhdrpWFDAEi",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fluid_velocity_df = fvelocity_second_rate # the results dataframe to be plotted\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,8), sharey=True)\n",
    "\n",
    "fig.suptitle('Evaluate the Quaility of the Raw Cross-plot Analysis', fontsize=18)\n",
    "\n",
    "ax1.set_title('R-squared', fontsize=16)\n",
    "im1 = ax1.scatter(fluid_velocity_df.intercept_velocity_mps, fluid_velocity_df.depth_m,\n",
    "    c = fluid_velocity_df.r_squared, s = 20, linewidths = 0)\n",
    "fig.colorbar(im1,ax=ax1)\n",
    "\n",
    "ax2.set_title('Number of data points', fontsize=16)\n",
    "im2 = ax2.scatter(fluid_velocity_df.intercept_velocity_mps, fluid_velocity_df.depth_m,\n",
    "    c = fluid_velocity_df.obs_num, s = 20, linewidths = 0)\n",
    "fig.colorbar(im2,ax=ax2)\n",
    "\n",
    "ax1.set_ylabel('Depth [m]')\n",
    "ax1.set_ylim(950,400)\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.set_xlabel('Fluid velocity [m/s]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "fNA3RH9YGEqhoxxPhoK1",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 18.5 Clean the fluid velocity data\n",
    "\n",
    "We will use the number of values in the cross-plot interval and the $R^2$ value to remove from our results dataframe those values that are likely to be suspect. Visually inspect the plots above to set limits on these filters.\n",
    "\n",
    "#### 18.5.1 First pump flow rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "TLsIrCJCW33vVyJesAqC",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# first rate\n",
    "\n",
    "# Filter data based on R2 value\n",
    "fvelocity_first_rate_trimmed = fvelocity_first_rate[(fvelocity_first_rate.r_squared > 0.98 )] # check filter\n",
    "\n",
    "# filter data based on number of values\n",
    "fvelocity_first_rate_trimmed = fvelocity_first_rate_trimmed[(fvelocity_first_rate_trimmed.obs_num > 6 )] # check filter\n",
    "\n",
    "\n",
    "print('before filter =', fvelocity_first_rate.shape, 'after filter =', fvelocity_first_rate_trimmed.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "XbkoTjqyITS8AUMbjxsg",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### 18.5.2 Second pump flow rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "qv7inyuxNRemhvgKA8NZ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# second rate\n",
    "\n",
    "# Filter data based on R2 value\n",
    "fvelocity_second_rate_trimmed = fvelocity_second_rate[(fvelocity_second_rate.r_squared > 0.98 )]\n",
    "\n",
    "# filter data based on number of values\n",
    "fvelocity_second_rate_trimmed = fvelocity_second_rate_trimmed[(fvelocity_second_rate_trimmed.obs_num > 6 )]\n",
    "\n",
    "print('before filter =', fvelocity_second_rate.shape, 'after filter =', fvelocity_second_rate_trimmed.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "Rt3jkrOocl21Na337vVM",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### 18.5.3 Third pump flow rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "FD2DaaRhXdI2sgC2iaVE",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# third rate\n",
    "\n",
    "# Filter data based on R2 value\n",
    "fvelocity_third_rate_trimmed = fvelocity_third_rate[(fvelocity_third_rate.r_squared > 0.98 )]\n",
    "\n",
    "# filter data based on number of values\n",
    "fvelocity_third_rate_trimmed = fvelocity_third_rate_trimmed[(fvelocity_third_rate_trimmed.obs_num > 6 )]\n",
    "\n",
    "print('before filter =', fvelocity_third_rate.shape, 'after filter =', fvelocity_third_rate_trimmed.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "f8BTisVTxOm31a0eMGH3",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 18.4 Visualise cleaned fluid velocity results\n",
    "\n",
    "Below we generate the same plot as used in Section 18.4 but with the dataframes that have had suspect cross-plot analysis results removed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "88OnB0jhXa9yI51t9NOS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fluid_velocity_df = fvelocity_first_rate_trimmed # the results dataframe to be plotted\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,8), sharey=True)\n",
    "\n",
    "fig.suptitle('Evaluate the Quaility of the Cross-plot Analysis after Filtering', fontsize=18)\n",
    "\n",
    "ax1.set_title('R-squared', fontsize=16)\n",
    "im1 = ax1.scatter(fluid_velocity_df.intercept_velocity_mps, fluid_velocity_df.depth_m,\n",
    "    c = fluid_velocity_df.r_squared, s = 20, linewidths = 0)\n",
    "fig.colorbar(im1,ax=ax1)\n",
    "\n",
    "ax2.set_title('Number of data points', fontsize=16)\n",
    "im2 = ax2.scatter(fluid_velocity_df.intercept_velocity_mps, fluid_velocity_df.depth_m,\n",
    "    c = fluid_velocity_df.obs_num, s = 20, linewidths = 0)\n",
    "fig.colorbar(im2,ax=ax2)\n",
    "\n",
    "ax1.set_ylabel('Depth [m]')\n",
    "ax1.set_ylim(950,400)\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.set_xlabel('Fluid velocity [m/s]')   \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "f9RpkiyODtBO42auy0Ay",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### 18.5 Visulise each cross-plot analysis\n",
    "\n",
    "The function below enables us to plot each cross-plot and save it to a folder so we can visually check the model fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "mdy13CAchtmBrMqsGIL0",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "def get_df_name(df):\n",
    "    '''detect the dataframe name because sometimes df.name does not work'''\n",
    "    name =[x for x in globals() if globals()[x] is df][0]\n",
    "    return name\n",
    "\n",
    "def make_crossplot_figures(dataframe, filename):\n",
    "    '''Export all cross-plots for a dataframe into a folder in our present working directory\n",
    "    \n",
    "    WARNING: A folder will be made with the specified filename and\n",
    "    if a folder already exists with that filename, it will be overwritten'''\n",
    "    \n",
    "    # make or overwite folder\n",
    "    import os\n",
    "    import shutil\n",
    "    directory = filename\n",
    "    parent_dir = os.getcwd() # detect present working directory\n",
    "    path = os.path.join(parent_dir, directory)\n",
    "    if os.path.exists(path): # if the directory already exists, remove it\n",
    "        shutil.rmtree(path)\n",
    "    os.mkdir(path) # make a directory in our present working directiory with our filename\n",
    "    \n",
    "    # make cross-plot for each row in the dataframe\n",
    "    for i, row in dataframe.iterrows():\n",
    "        frequency_data_hz = dataframe.iloc[i]['frequency_data_hz'] #.values.tolist()\n",
    "        speed_data_mps = dataframe.iloc[i]['speed_data_mps'] #.values.tolist()\n",
    "        slope = dataframe.iloc[i]['slope']\n",
    "        intercept = dataframe.iloc[i]['intercept_velocity_mps']\n",
    "\n",
    "        # calculate y values using our model\n",
    "        model_y_vals = []\n",
    "        for n in frequency_data_hz:\n",
    "            model_y_vals.append(slope * n + intercept)\n",
    "\n",
    "        # generate test plot\n",
    "        fig, (ax) = plt.subplots(1, 1,figsize=(8,8))\n",
    "        ax.set_title('Dataframe = {df_name}, Depth = {depth}, \\n R2 = {rsquared:.4f}, Datapoint num = {datanumber}'.format(\n",
    "            df_name = get_df_name(dataframe),\n",
    "            depth = dataframe.iloc[i]['depth_m'], \n",
    "            rsquared = dataframe.iloc[i]['r_squared'], \n",
    "            datanumber = dataframe.iloc[i]['obs_num']))\n",
    "\n",
    "        ax.scatter(frequency_data_hz, speed_data_mps,\n",
    "            color = 'none', edgecolors = 'k', marker='o', s = 50, linewidths = 1)\n",
    "\n",
    "        ax.plot(frequency_data_hz, model_y_vals, \n",
    "            color='tab:orange', linestyle='-', linewidth=3, alpha=0.5, label='Linear fit')\n",
    "\n",
    "        ax.hlines(0, -40, 40, color = 'k', linewidth = 0.5)\n",
    "        ax.vlines(0, -2, 2, color = 'k', linewidth = 0.5)\n",
    "\n",
    "        ax.set_xlim(-40,40)\n",
    "        ax.set_ylim(-2,2)\n",
    "\n",
    "        ax.set_xlabel('Spinner frequency [hz]')\n",
    "        ax.set_ylabel('Tool speed [mps]')\n",
    "\n",
    "        ax.grid()\n",
    "        \n",
    "        plt.savefig(path + '/{depth}.png'.format(depth = dataframe.iloc[i]['depth_m']), \n",
    "                    dpi = 300,\n",
    "                    facecolor='white', transparent=False)\n",
    "        plt.close()\n",
    "\n",
    "    return print('Cross-plot figures are saved in {folder}'.format(folder = path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "IJ5sVL84PEPV55KHRYR0",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "**Take care with the filename term in the make_crossplot_figures function: If this filename already exists, it will be deleted and a new one made.**\n",
    "\n",
    "Uncomment the code in this cell to export the cross-plots. This code will take quite a while to run because we are making many figures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "DmMnNVFaFZmBrNR73R8W",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# I don't reccomend running this cell in Google Colab\n",
    "#make_crossplot_figures(fvelocity_second_rate, 'xplots_second_rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "1pfB2LnE1wFCEN0SAsBu",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 19. Combine the results and find the feed zones\n",
    "\n",
    "In this section, we combine all the data together using a quite richly formatted plot to demonstrate the completion test data visualisation that's possible with Python. \n",
    "\n",
    "For the purposes of this tutorial the feed zone interpretation is made using only temperature and fluid velocity profiles during injection, and the first heating profile. The full set of heat-up data and fracture data are not included.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "7T6EVxQ1CmcSTMXdQL4W",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "heating_37days = pd.read_csv(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-Temp-Heating37days.csv') \n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#heating_37days = pd.read_csv('Data-Temp-Heating37days.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "nsOgjesYyoolW8OLYH95",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Convert bar gauge to bar atmosphere\n",
    "\n",
    "heating_37days['pressure_bara'] = heating_37days.pres_barg - 1\n",
    "\n",
    "heating_37days.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "CWuUDueqDQEw3NI9ZrHP",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the BPD\n",
    "\n",
    "# note that iapws uses SI units so some unit conversion is required\n",
    "\n",
    "heating_37days['pressure_mpa'] = heating_37days.pressure_bara * 0.1  # convert pressure to MPa for ipaws\n",
    "\n",
    "pressure = heating_37days['pressure_mpa'].tolist()\n",
    "tsat = []\n",
    "for p in pressure:\n",
    "    saturation_temp = iapws.iapws97._TSat_P(p) - 273.15  # calculate saturation temp in Kelvin & convert to degC\n",
    "    tsat.append(saturation_temp)\n",
    "heating_37days['tsat_degC'] = tsat\n",
    "\n",
    "heating_37days.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "ruYrGznN4JBI47g0FCHh",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,10), sharey=True)\n",
    "\n",
    "# feedzone interpretaion\n",
    "feedzones = [\n",
    "    (560, 620), # 1\n",
    "    (660, 690), # 2\n",
    "    (705, 715), # 3\n",
    "    (720, 740), # 4\n",
    "    (745, 775), # 5\n",
    "    (800, 835), # 6\n",
    "    (850, 875), # 7\n",
    "    (920, 930)  # 8\n",
    "    ]\n",
    "\n",
    "for ax in [ax1, ax2]: # plot all FZ\n",
    "    for top, bottom in feedzones:\n",
    "        ax.axhspan(top, bottom, color='tab:blue', alpha=.1)\n",
    "\n",
    "biggest_feedzones = [\n",
    "    (745, 775), # 5\n",
    "    (800, 835), # 6\n",
    "    ]\n",
    "\n",
    "for ax in [ax1, ax2]: # highlight the largest FZ\n",
    "    for top, bottom in biggest_feedzones:\n",
    "        ax.axhspan(top, bottom, color='tab:blue', alpha=.4)\n",
    "\n",
    "label_depth = [] # find the half way point for label depth\n",
    "for top, bottom in feedzones:\n",
    "    l = (bottom - top)/2 + top\n",
    "    label_depth.append(l)\n",
    "\n",
    "labels = ['FZ1','FZ2','FZ3','FZ4','FZ5','FZ6','FZ7','FZ8']    \n",
    "\n",
    "for depth, label in zip(label_depth,labels): # plot FZ labels\n",
    "    ax1.text(-0.16,depth,label,verticalalignment='center')\n",
    "\n",
    "\n",
    "# fluid veolocity profiles for each flow rate generate by cross-plot analaysis\n",
    "ax1.plot(fvelocity_first_rate_trimmed.intercept_velocity_mps,fvelocity_first_rate_trimmed.depth_m, \n",
    "    #marker = '.', # uncomment this to view the data points\n",
    "    color = '#440154', linestyle = '-', linewidth = 2, #alpha = 0.8,\n",
    "    label = '{rate:.0f} t/hr injection - velocity'.format(rate = pts_first_rate.flow_tph.mean()))\n",
    "ax1.plot(fvelocity_third_rate_trimmed.intercept_velocity_mps,fvelocity_third_rate_trimmed.depth_m,\n",
    "    #marker = '.',\n",
    "    color = '#5ec962', linestyle = '-', linewidth = 2, #alpha = 0.8,\n",
    "    label = '{rate:.0f} t/hr injection - velocity'.format(rate = pts_third_rate.flow_tph.mean())\n",
    "    )\n",
    "ax1.plot(fvelocity_second_rate_trimmed.intercept_velocity_mps,fvelocity_second_rate_trimmed.depth_m,\n",
    "    #marker = '.',\n",
    "    color = '#21918c', linestyle = '-', linewidth = 2, #alpha = 0.8,\n",
    "    label = '{rate:.0f} t/hr injection - veolcity'.format(rate = pts_second_rate.flow_tph.mean()))\n",
    "\n",
    "# completion test temp data\n",
    "ax2.scatter(clean_pts.temp_degC, clean_pts.depth_m,     \n",
    "    c = clean_pts.timedelta_sec, s = 5, linewidths = 0, alpha = 0.5)\n",
    "\n",
    "# false plots to generate the legand\n",
    "ax2.plot(0, 0, color = '#440154', linewidth = 2, # purple\n",
    "    label = '{rate:.0f} t/hr injection - temp'.format(rate = pts_first_rate.flow_tph.mean()))\n",
    "ax2.plot(0, 0, color = '#5ec962', linewidth = 2, # blue\n",
    "    label = '{rate:.0f} t/hr injection - temp'.format(rate = pts_third_rate.flow_tph.mean()))\n",
    "ax2.plot(0, 0, color = '#21918c', linewidth = 2, # green\n",
    "    label = '{rate:.0f} t/hr injection - temp'.format(rate = pts_second_rate.flow_tph.mean()))\n",
    "ax2.plot(0, 0, color = '#fde725', linewidth = 2, # yellow\n",
    "    label = 'Day 0 shut - temp')\n",
    "\n",
    "\n",
    "# stable temp data\n",
    "ax2.plot(heating_37days.temp_degC, heating_37days.depth_m, \n",
    "    color = '#fd7b25', linewidth = 2,\n",
    "    label = 'Day 37 shut - temp')\n",
    "\n",
    "# saturation temp for the stable pressure profile assuming pure water\n",
    "ax2.plot(heating_37days.tsat_degC, heating_37days.depth_m, \n",
    "    linestyle = ':', color = 'k', linewidth = 2,\n",
    "    label = 'Day 37 shut - BPD')\n",
    "\n",
    "\n",
    "production_shoe = 462.5 # 13 3/8 production casing shoe in meters measured depth (mMD) from the casing head flange (CHF)\n",
    "top_of_liner = 425 # top of perforated 10 3/4 liner in meters measured depth (mMD) from CHF\n",
    "terminal_depth = 946 # deepest drilled depth \n",
    "# the perforated liner is squatted on bottom but didn't quite make it all the way down (bottom of liner is 931 mMD)\n",
    "\n",
    "# blank well casing\n",
    "ax1.plot([-0.2, -0.2],[0, production_shoe],\n",
    "    color = 'k', linewidth = 8, linestyle = '-')\n",
    "ax2.plot([1, 1],[0, production_shoe],\n",
    "    color = 'k', linewidth = 3, linestyle = '-')\n",
    "\n",
    "# perforated well casing\n",
    "ax1.plot([-0.18, -0.18],[top_of_liner, terminal_depth],     \n",
    "    color = 'k', linewidth = 1.5, linestyle = '--')\n",
    "ax2.plot([5, 5],[top_of_liner, terminal_depth],     \n",
    "    color = 'k', linewidth = 1.5, linestyle = '--')\n",
    "\n",
    "\n",
    "ax1.set_xlim(-0.2,1.2)\n",
    "ax1.set_xlabel('Fluid velocity [m/s]')\n",
    "\n",
    "ax2.set_xlim(0,300)\n",
    "ax2.set_xlabel('Temperature [degC]')\n",
    "\n",
    "ax1.set_ylim(950,300) #950,300 to show produciton zone\n",
    "ax1.set_ylabel('Depth [m]')\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.grid()\n",
    "    ax.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "yZjZp8CDtNCelFGMWUPS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### This is not an easy well to interpret! \n",
    "\n",
    "- There are multiple feed zones and they interact differently at the various injection rates \n",
    "- Some feed zones behave the same regardless of injection rate, while in other feed zones the flow direction switches when the injection rate is changed \n",
    "- The majority of the fluid is exiting at FZ5 and FZ6 (dark blue) \n",
    "- It is difficult to pinpoint whether there is a single major feed zone from this data alone (i.e., from the fluid velocity profiles and temperature profiles) \n",
    "- Pivoting of pressure profiles during progressive heat-up runs and combined analysis with borehole image log data would aid the interpretation, but this next analysis step is beyond the scope of our tutorial. \n",
    "\n",
    "\n",
    "|FZ number | Upper bound |  Lower bound | Features |\n",
    "| --- | --- | --- | --- |\n",
    "| 1 | 560 | 620 | Inflow at lowest injection rate and then outflow at the higher injection rates, as shown by: <br /> - Large increase fluid velocity at the low injection rate, which reverses for the other two higher rates, becoming a step down. <br /> - Increase in temperature gradient present for the first injection rate, which is not present for the two higher rates.|\n",
    "| 2 | 660 | 690 | Small and high-velocity inflow of two-phase fluid, as seen at the highest two flow rates from: <br /> - Spikes in fluid velocity which subside below the feed zone. <br /> - Increases in temperature gradient.  |\n",
    "| 3 | 705 | 715 | Small and high-velocity inflow of two-phase fluid, as seen at all three flow rates from: <br /> - Spikes in fluid velocity which subside below the feed zone. <br /> - Increases in temperature gradient. <br /> - Rapid heating at this depth after injection stops.  |\n",
    "| 4 | 720 | 740 | Small and high-velocity inflow of two-phase fluid, as seen at all three flow rates from: <br /> - Spikes in fluid velocity which subside below the feed zone. <br /> - Increases in temperature gradient. <br /> - Rapid heating at this depth after injection stops. |\n",
    "| 5 | 745 | 775 | Outflow of fluid from the wellbore, as seen at all three flow rates from: <br /> - Drop in fluid velocity. |\n",
    "| 6 | 800 | 835 | Outflow of fluid from the wellbore, as seen at all three flow rates from: <br /> - Drop in fluid velocity. |\n",
    "| 7 | 850 | 875 | Outflow of fluid from the wellbore, as seen at all three flow rates from: <br /> - Drop in fluid velocity. <br /> - Increase in temperature gradient. <br /> - Anomaly in heatup temperatures. |\n",
    "| 8 | 920 | 930 | Must be some minor permeability at or below this depth, due to separation of the temperature profiles at different injection rates. |\n",
    "\n",
    "***\n",
    "\n",
    "#### What's next for completion test analysis with Python?\n",
    "\n",
    "Future possible refinements of this fluid velocity method:\n",
    "- If an open-hole calliper log was run:\n",
    "    - We could limit the fluid velocity analysis to only those intervals that are in gauge.\n",
    "    - Or correct the fluid velocity profile for hole size (and compare to the spinner ratio method currently used to achieve the same thing).\n",
    "- Perhaps we could look for data where the spinner was stuck and remove these.\n",
    "- Set a condition on the cross-plot analysis that it does not return an intercept if there are only negative or positive values for spinner frequency.\n",
    "- Use root mean squared error instead of $R^2$ and consider other measures of the quality of our model fit.\n",
    "\n",
    "Do you have any suggestions?\n",
    "\n",
    "Send your comments and feedback to irene@cubicearth.nz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "UbCMQvEtSv8FCdwHDF5x",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "You have finished the T21 geothermal well completion test tutoral. Well done!\n",
    "\n",
    "***\n",
    "\n",
    "<p><center>© 2021 <a href=\"https://www.cubicearth.nz/\">Irene Wallis</a> and <a href=\"https://www.linkedin.com/in/katie-mclean-25994315/\">Katie McLean</a> <a href=\"https://creativecommons.org/licenses/by/4.0/\"</a></center></p>\n",
    "\n",
    "<p><center>Licensed under the Apache License, Version 2.0</center></p>\n",
    "\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "iooxa": {
   "id": {
    "block": "Yq1SQWsDsoJ2sJYwDcNB",
    "project": "mmReuqVTAa9JzPpNr22I",
    "version": 1
   }
  },
  "kernelspec": {
   "display_name": "geothrm",
   "language": "python",
   "name": "geothrm"
  },
  "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}