post_mask_data_reduction_single_pipeline_master.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PIPELINE FOR SINGLE RASTER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# §0 - Read-in and set-up\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.1 - file read-in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your file iris2model.sav contains an inversion of a single raster and will be returned as a numpy.record...\n",
      "...with the required structure\n"
     ]
    }
   ],
   "source": [
    "from sav_extract import sav_extract\n",
    "filename='iris2model.sav'\n",
    "iris2model=sav_extract(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.2 - Initial data read-in and trim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iris_cal has successfully been trimmed to (217, 183, 32)\n"
     ]
    }
   ],
   "source": [
    "from data_extract import tag_extract, kpeak_ind_extract, trimmed_rows, data_trim\n",
    "\n",
    "wavelength  = tag_extract(iris2model,'wl') # need this one first\n",
    "iris_cal    = tag_extract(iris2model,'iris_cal') # then this one to get trimmed rows from kpeak frame\n",
    "\n",
    "kpeak_ind   = kpeak_ind_extract(wavelength) # isolate the wavelength index closest to Mg II k peak\n",
    "frame       = iris_cal[kpeak_ind,...] # create a frame at this peak wavelength, where FL is most visible\n",
    "trimmed_rows= trimmed_rows(frame) # cuts FOV to only include observed data, not zeros at top and bottom (returns list of top and bottom row)\n",
    "\n",
    "del frame # it has served its purpose to get trimmed rows, now delete to save confusion pre-trim\n",
    "del wavelength\n",
    "\n",
    "iris_cal    = data_trim(iris_cal, trimmed_rows,'iris_cal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.3 - Initial frame creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame=iris_cal[kpeak_ind]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# §1 - Frame cleaning function list\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 - Slit row extraction function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "def slit_row_extract(frame):\n",
    "    '''\n",
    "    Extract the index of the row with the slit in it from monochromatic image.\n",
    "    Works by getting the row index of the min val in each col and taking the most common\n",
    "    row value from that list\n",
    "    '''\n",
    "    slit_rows=[]\n",
    "    for col in range(frame.shape[1]):\n",
    "        slit_rows.append(np.argmin(frame[:,col]))\n",
    "    slit_row=int(stats.mode(slit_rows)[0])\n",
    "    return slit_row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 - Threshold extraction function (using hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np # explicitly state package dependencies for each sep function\n",
    "from scipy.signal import find_peaks\n",
    "\n",
    "def thresh_extract(frame,slit_row):\n",
    "    '''\n",
    "    Extract flat threshold frame for filament sep from k3 monochromatic frame.\n",
    "    It does this by raising slit vals to a max value so the first trough of histogram \n",
    "    is roughly the threshold between filamenty and quiet sun.\n",
    "    (Only tested for quiet sun filament obs as of 21/09/2020)\n",
    "    '''\n",
    "    ### raise slit vals to max for hist peaks\n",
    "    frame[slit_row,:]=np.amax(frame) \n",
    "    frame[slit_row-1,:]=np.amax(frame)\n",
    "\n",
    "    ### then do find peaks for thresh here:\n",
    "    frame_linear=np.reshape(frame,-1)\n",
    "    bin_number=int(len(frame_linear)/80) # 80 decided empirically to give consistently good peaks\n",
    "    hist,bins = np.histogram(frame_linear, bins = bin_number) # iteratively decided on 70 bins\n",
    "    bin_centres = 0.5*(bins[1:]+bins[:-1])\n",
    "\n",
    "    peaks, _ = find_peaks(-hist) # find troughs\n",
    "    index=peaks[0] # since bin size is so big, first trough corresponds to post-FL dip.\n",
    "    thresh=bin_centres[index] \n",
    "    frame[frame>thresh]=0\n",
    "    return frame\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 - Binary mask extraction function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import medfilt\n",
    "\n",
    "def bin_mask_extract(frame,slit_row=None):\n",
    "    '''\n",
    "    Extracts binary mask array-like 1s for filament, 0s for quiet sun from post-thresh frame.\n",
    "    It does this using a basic median filter that cleans out rough edges of post-thresh frame.\n",
    "    For border creation, don't give slit_row keyword (COSMETIC)\n",
    "    For index sep and further data reduction, give slit_row keyword (DATA REDUCTION)\n",
    "    '''\n",
    "    bin_mask=np.where(frame>0,1,0)\n",
    "    if slit_row is None: # COSMETIC no slits nan'd to create line border\n",
    "        bin_mask=medfilt(bin_mask,5) \n",
    "        return bin_mask # COSMETIC\n",
    "    else:\n",
    "        # DATA REDUCTION slit rows are nan'd to stop slit row values (which are neither FL or QS) being wronlgy counted as such\n",
    "        # need to give slit row value as explicit argument when calling function\n",
    "        bin_mask[slit_row,:]=100 \n",
    "        bin_mask[slit_row-1,:]=100\n",
    "        bin_mask=np.where(bin_mask>1,np.nan,bin_mask)\n",
    "        bin_mask=medfilt(bin_mask,5) \n",
    "        return bin_mask # Further DATA REDUCTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 - Mask border extraction function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.collections import LineCollection\n",
    "\n",
    "def get_all_boundary_edges(bool_img):\n",
    "    '''\n",
    "    Get a list of all edges\n",
    "    (where the value changes from 'True' to 'False') in the 2D image.\n",
    "    Return the list as indices of the image.\n",
    "    '''\n",
    "    ij_boundary = []\n",
    "    ii, jj = np.nonzero(bool_img)\n",
    "    for i, j in zip(ii, jj):\n",
    "        # North\n",
    "        if j == bool_img.shape[1]-1 or not bool_img[i, j+1]:\n",
    "            ij_boundary.append(np.array([[i, j+1],\n",
    "                                         [i+1, j+1]]))\n",
    "        # East\n",
    "        if i == bool_img.shape[0]-1 or not bool_img[i+1, j]:\n",
    "            ij_boundary.append(np.array([[i+1, j],\n",
    "                                         [i+1, j+1]]))\n",
    "        # South\n",
    "        if j == 0 or not bool_img[i, j-1]:\n",
    "            ij_boundary.append(np.array([[i, j],\n",
    "                                         [i+1, j]]))\n",
    "        # West\n",
    "        if i == 0 or not bool_img[i-1, j]:\n",
    "            ij_boundary.append(np.array([[i, j],\n",
    "                                         [i, j+1]]))\n",
    "    if not ij_boundary:\n",
    "        return np.zeros((0, 2, 2))\n",
    "    else:\n",
    "        return np.array(ij_boundary)\n",
    "\n",
    "def close_loop_boundary_edges(xy_boundary, clean=True):\n",
    "    '''\n",
    "    Connect all edges defined by 'xy_boundary' to closed \n",
    "    boundary lines.\n",
    "    If not all edges are part of one surface return a list of closed \n",
    "    boundaries is returned (one for every object).\n",
    "    '''\n",
    "    boundary_loop_list = []\n",
    "    while xy_boundary.size != 0:\n",
    "        # Current loop\n",
    "        xy_cl = [xy_boundary[0, 0], xy_boundary[0, 1]]  # Start with first edge\n",
    "        xy_boundary = np.delete(xy_boundary, 0, axis=0)\n",
    "\n",
    "        while xy_boundary.size != 0:\n",
    "            # Get next boundary edge (edge with common node)\n",
    "            ij = np.nonzero((xy_boundary == xy_cl[-1]).all(axis=2))\n",
    "            if ij[0].size > 0:\n",
    "                i = ij[0][0]\n",
    "                j = ij[1][0]\n",
    "            else:\n",
    "                xy_cl.append(xy_cl[0])\n",
    "                break\n",
    "\n",
    "            xy_cl.append(xy_boundary[i, (j + 1) % 2, :])\n",
    "            xy_boundary = np.delete(xy_boundary, i, axis=0)\n",
    "\n",
    "        xy_cl = np.array(xy_cl)\n",
    "\n",
    "        boundary_loop_list.append(xy_cl)\n",
    "\n",
    "    return boundary_loop_list\n",
    "\n",
    "def plot_world_outlines(bool_img, ax=None, **kwargs):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    ij_boundary = get_all_boundary_edges(bool_img=bool_img)\n",
    "    xy_boundary = ij_boundary - 0.5\n",
    "    xy_boundary = close_loop_boundary_edges(xy_boundary=xy_boundary)\n",
    "    cl = LineCollection(xy_boundary, **kwargs)\n",
    "    ax.add_collection(cl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# §2 - Data reduction function list\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 - Pixel region assingment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def pix_region_assign(bin_mask):\n",
    "    '''\n",
    "    Takes input of binary mask with slit values as nans and returns two lists as index tuples \n",
    "    for pixels belonging to fl and qs regions (in that order) only if all pix are accounted for.\n",
    "    Slit pixels are added to a seperate list but not returned. This is created to ensure all \n",
    "    pixels of bin mask are accounted for but serve no further purpose.\n",
    "    '''\n",
    "    fl_pix=[]\n",
    "    qs_pix=[]\n",
    "    slit_pix=[]\n",
    "    for pix_ind, pix_val in np.ndenumerate(bin_mask):\n",
    "        if pix_val==0:\n",
    "            qs_pix.append(pix_ind)\n",
    "        if pix_val==1:\n",
    "            fl_pix.append(pix_ind)\n",
    "        if np.isnan(pix_val)==True:\n",
    "            slit_pix.append(pix_ind)\n",
    "\n",
    "    if len(fl_pix+qs_pix+slit_pix)==len(np.reshape(bin_mask,-1)):\n",
    "        print('All binary mask pixels successfully accounted for...')\n",
    "        return fl_pix,qs_pix        \n",
    "    else:\n",
    "        print('Not all binary pixels are accounted for...')\n",
    "        return "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 - Data and unc mean (removing bad fits - NEW METHOD) - *IMPORANT*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def data_mean_extract(model_dat_cube,pixels):\n",
    "    '''\n",
    "    Creates arrays for model means and their uncertainties for QS/Fl seperately\n",
    "    (model is either temp or ne, pixels is eiher qs or fl list)\n",
    "    '''\n",
    "    model_sum=np.zeros(model_dat_cube.shape[0]) # makes zero array like ltau\n",
    "    for tuple_index in pixels:\n",
    "        model_ltau=model_dat_cube[:,tuple_index[0],tuple_index[1]] #for all values of ltau but single pix\n",
    "        print(type(model_ltau),model_ltau.shape)\n",
    "        np.isfinite()\n",
    "        # want to skip pixels that have nans or infs\n",
    "        \n",
    "        model_ltau=np.nan_to_num(model_ltau,nan=0,posinf=0) #some uncs were 'inf' which messed up the sum, so for these I have removed them.\n",
    "        model_sum=model_sum+model_ltau\n",
    "    model_mean=model_sum/len(pixels)    \n",
    "    return model_mean"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "iris2_venv",
   "language": "python",
   "name": "iris2_venv"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}