NDVI_smoothing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NDVI SMOOTHING ANALYSIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import array\n",
    "import datetime\n",
    "from datetime import timedelta\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import pickle\n",
    "from termcolor import colored\n",
    "import time\n",
    "import xarray as xr\n",
    "\n",
    "%reload_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cython functions, modified `ws2dvopt` and `ws2dvoptp` in order to get vcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from cpython.array cimport array, clone\n",
    "from libc.math cimport log, pow, sqrt\n",
    "cimport numpy as np\n",
    "import numpy as np\n",
    "\n",
    "tFloat = np.double\n",
    "ctypedef np.double_t dtype_t\n",
    "\n",
    "\n",
    "cpdef lag1corr(np.ndarray[dtype_t] data1, np.ndarray[dtype_t] data2, double nd):\n",
    "    \"\"\"Calculates Lag-1 autocorrelation.\n",
    "\n",
    "    Adapted from https://stackoverflow.com/a/29194624/5997555\n",
    "\n",
    "    Args:\n",
    "        data1: fist data series\n",
    "        data2: second data series\n",
    "        nd: no-data value (will be exluded from calulation)\n",
    "\n",
    "    Returns:\n",
    "        Lag-1 autocorrelation value\n",
    "    \"\"\"\n",
    "\n",
    "    cdef int M, sub\n",
    "    cdef double sum1, sum2, var_sum1, var_sum2, cross_sum, std1, std2, cross_mean\n",
    "\n",
    "    M = data1.size\n",
    "\n",
    "    sum1 = 0.\n",
    "    sum2 = 0.\n",
    "    sub = 0\n",
    "    for i in range(M):\n",
    "        if data1[i] != nd and data2[i] != nd:\n",
    "            sum1 += data1[i]\n",
    "            sum2 += data2[i]\n",
    "        else:\n",
    "            sub += 1\n",
    "    mean1 = sum1 / (M-sub)\n",
    "    mean2 = sum2 / (M-sub)\n",
    "\n",
    "    var_sum1 = 0.\n",
    "    var_sum2 = 0.\n",
    "    cross_sum = 0.\n",
    "    for i in range(M):\n",
    "        if data1[i] != nd and data2[i] != nd:\n",
    "            var_sum1 += (data1[i] - mean1) ** 2\n",
    "            var_sum2 += (data2[i] - mean2) ** 2\n",
    "            cross_sum += (data1[i] * data2[i])\n",
    "\n",
    "    std1 = (var_sum1 / (M-sub)) ** .5\n",
    "    std2 = (var_sum2 / (M-sub)) ** .5\n",
    "    cross_mean = cross_sum / (M-sub)\n",
    "    return (cross_mean - mean1 * mean2) / (std1 * std2)\n",
    "\n",
    "cpdef ws2d(np.ndarray[dtype_t] y, double lmda, np.ndarray[dtype_t] w):\n",
    "    cdef array dbl_array_template = array('d', [])\n",
    "    cdef int i, i1, i2, m, n\n",
    "    cdef array z, d, c, e\n",
    "\n",
    "    n = y.shape[0]\n",
    "    m = n - 1\n",
    "\n",
    "    z = clone(dbl_array_template, n, zero=False)\n",
    "    d = clone(dbl_array_template, n, zero=False)\n",
    "    c = clone(dbl_array_template, n, zero=False)\n",
    "    e = clone(dbl_array_template, n, zero=False)\n",
    "\n",
    "    d.data.as_doubles[0] = w[0] + lmda\n",
    "    c.data.as_doubles[0] = (-2 * lmda) / d.data.as_doubles[0]\n",
    "    e.data.as_doubles[0] = lmda /d.data.as_doubles[0]\n",
    "    z.data.as_doubles[0] = w[0] * y[0]\n",
    "    d.data.as_doubles[1] = w[1] + 5 * lmda - d.data.as_doubles[0] * (c.data.as_doubles[0] * c.data.as_doubles[0])\n",
    "    c.data.as_doubles[1] = (-4 * lmda - d.data.as_doubles[0] * c.data.as_doubles[0] * e.data.as_doubles[0]) / d.data.as_doubles[1]\n",
    "    e.data.as_doubles[1] =  lmda / d.data.as_doubles[1]\n",
    "    z.data.as_doubles[1] = w[1] * y[1] - c.data.as_doubles[0] * z.data.as_doubles[0]\n",
    "    for i in range(2, m-1):\n",
    "        i1 = i - 1\n",
    "        i2 = i - 2\n",
    "        d.data.as_doubles[i]= w[i] + 6 *  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "        c.data.as_doubles[i] = (-4 *  lmda - d.data.as_doubles[i1] * c.data.as_doubles[i1] * e.data.as_doubles[i1])/ d.data.as_doubles[i]\n",
    "        e.data.as_doubles[i] =  lmda / d.data.as_doubles[i]\n",
    "        z.data.as_doubles[i] = w[i] * y[i] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]\n",
    "    i1 = m - 2\n",
    "    i2 = m - 3\n",
    "    d.data.as_doubles[m - 1] = w[m - 1] + 5 *  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "    c.data.as_doubles[m - 1] = (-2 *  lmda - d.data.as_doubles[i1] * c.data.as_doubles[i1] * e.data.as_doubles[i1]) / d.data.as_doubles[m - 1]\n",
    "    z.data.as_doubles[m - 1] = w[m - 1] * y[m - 1] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]\n",
    "    i1 = m - 1\n",
    "    i2 = m - 2\n",
    "    d.data.as_doubles[m] = w[m] +  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "    z.data.as_doubles[m] = (w[m] * y[m] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]) / d.data.as_doubles[m]\n",
    "    z.data.as_doubles[m - 1] = z.data.as_doubles[m - 1] / d.data.as_doubles[m - 1] - c.data.as_doubles[m - 1] * z.data.as_doubles[m]\n",
    "    for i in range(m-2, -1, -1):\n",
    "        z.data.as_doubles[i] = z.data.as_doubles[i] / d.data.as_doubles[i] - c.data.as_doubles[i] * z.data.as_doubles[i + 1] - e.data.as_doubles[i] * z.data.as_doubles[i + 2]\n",
    "    return z\n",
    "\n",
    "cpdef ws2dp(np.ndarray[dtype_t] y, double lmda, np.ndarray[dtype_t] w, double p):\n",
    "  \"\"\"Whittaker smoother with asymmetric smoothing and fixed lambda (S).\n",
    "\n",
    "  Args:\n",
    "      y: time-series numpy array\n",
    "      l: smoothing parameter lambda (S)\n",
    "      w: weights numpy array\n",
    "      p: \"Envelope\" value\n",
    "\n",
    "  Returns:\n",
    "      Smoothed time-series array z\n",
    "  \"\"\"\n",
    "  cdef array template = array('d', [])\n",
    "  cdef int m, i, j\n",
    "  cdef double y_tmp, z_tmp, p1\n",
    "\n",
    "  m = y.shape[0]\n",
    "  i = 0\n",
    "  j = 0\n",
    "  p1 = 1-p\n",
    "\n",
    "  template = array('d', [])\n",
    "  z = clone(template, m, True)\n",
    "  znew = clone(template, m, True)\n",
    "  wa = clone(template, m, False)\n",
    "  ww = clone(template, m, False)\n",
    "\n",
    "  # Calculate weights\n",
    "\n",
    "  for i in range(10):\n",
    "    for j in range(m):\n",
    "      y_tmp = y[j]\n",
    "      z_tmp = z.data.as_doubles[j]\n",
    "\n",
    "      if y_tmp > z_tmp:\n",
    "        wa.data.as_doubles[j] = p\n",
    "      else:\n",
    "        wa.data.as_doubles[j] = p1\n",
    "      ww.data.as_doubles[j] = w[j] * wa.data.as_doubles[j]\n",
    "\n",
    "    znew[0:m] = _ws2d(y, lmda, ww)\n",
    "    z_tmp = 0.0\n",
    "    j = 0\n",
    "    for j in range(m):\n",
    "      z_tmp += abs(znew.data.as_doubles[j] - z.data.as_doubles[j])\n",
    "\n",
    "    if z_tmp == 0.0:\n",
    "      break\n",
    "\n",
    "    z[0:m]= znew[0:m]\n",
    "\n",
    "  z[0:m] = _ws2d(y, lmda, ww)\n",
    "  return z\n",
    "\n",
    "\n",
    "cdef _ws2d(np.ndarray[dtype_t] y, double lmda, array[double] w):\n",
    "    \"\"\"Internal whittaker function for use in asymmetric smoothing.\n",
    "    Args:\n",
    "      y: time-series numpy array\n",
    "      lmbda: lambda (s) value\n",
    "      w: weights numpy array\n",
    "    Returns:\n",
    "        smoothed time-series array z\n",
    "    \"\"\"\n",
    "\n",
    "    cdef array dbl_array_template = array('d', [])\n",
    "    cdef int i, i1, i2, m, n\n",
    "    cdef array z, d, c, e\n",
    "\n",
    "    n = y.shape[0]\n",
    "    m = n - 1\n",
    "\n",
    "    z = clone(dbl_array_template, n, zero=False)\n",
    "    d = clone(dbl_array_template, n, zero=False)\n",
    "    c = clone(dbl_array_template, n, zero=False)\n",
    "    e = clone(dbl_array_template, n, zero=False)\n",
    "\n",
    "    d.data.as_doubles[0] = w.data.as_doubles[0] + lmda\n",
    "    c.data.as_doubles[0] = (-2 * lmda) / d.data.as_doubles[0]\n",
    "    e.data.as_doubles[0] = lmda /d.data.as_doubles[0]\n",
    "    z.data.as_doubles[0] = w.data.as_doubles[0] * y[0]\n",
    "    d.data.as_doubles[1] = w.data.as_doubles[1] + 5 * lmda - d.data.as_doubles[0] * (c.data.as_doubles[0] * c.data.as_doubles[0])\n",
    "    c.data.as_doubles[1] = (-4 * lmda - d.data.as_doubles[0] * c.data.as_doubles[0] * e.data.as_doubles[0]) / d.data.as_doubles[1]\n",
    "    e.data.as_doubles[1] =  lmda / d.data.as_doubles[1]\n",
    "    z.data.as_doubles[1] = w.data.as_doubles[1] * y[1] - c.data.as_doubles[0] * z.data.as_doubles[0]\n",
    "    for i in range(2, m-1):\n",
    "        i1 = i - 1\n",
    "        i2 = i - 2\n",
    "        d.data.as_doubles[i]= w.data.as_doubles[i] + 6 *  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "        c.data.as_doubles[i] = (-4 *  lmda - d.data.as_doubles[i1] * c.data.as_doubles[i1] * e.data.as_doubles[i1])/ d.data.as_doubles[i]\n",
    "        e.data.as_doubles[i] =  lmda / d.data.as_doubles[i]\n",
    "        z.data.as_doubles[i] = w.data.as_doubles[i] * y[i] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]\n",
    "    i1 = m - 2\n",
    "    i2 = m - 3\n",
    "    d.data.as_doubles[m - 1] = w.data.as_doubles[m - 1] + 5 *  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "    c.data.as_doubles[m - 1] = (-2 *  lmda - d.data.as_doubles[i1] * c.data.as_doubles[i1] * e.data.as_doubles[i1]) / d.data.as_doubles[m - 1]\n",
    "    z.data.as_doubles[m - 1] = w.data.as_doubles[m - 1] * y[m - 1] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]\n",
    "    i1 = m - 1\n",
    "    i2 = m - 2\n",
    "    d.data.as_doubles[m] = w.data.as_doubles[m] +  lmda - (c.data.as_doubles[i1] * c.data.as_doubles[i1]) * d.data.as_doubles[i1] - (e.data.as_doubles[i2] * e.data.as_doubles[i2]) * d.data.as_doubles[i2]\n",
    "    z.data.as_doubles[m] = (w.data.as_doubles[m] * y[m] - c.data.as_doubles[i1] * z.data.as_doubles[i1] - e.data.as_doubles[i2] * z.data.as_doubles[i2]) / d.data.as_doubles[m]\n",
    "    z.data.as_doubles[m - 1] = z.data.as_doubles[m - 1] / d.data.as_doubles[m - 1] - c.data.as_doubles[m - 1] * z.data.as_doubles[m]\n",
    "    for i in range(m-2, -1, -1):\n",
    "        z.data.as_doubles[i] = z.data.as_doubles[i] / d.data.as_doubles[i] - c.data.as_doubles[i] * z.data.as_doubles[i + 1] - e.data.as_doubles[i] * z.data.as_doubles[i + 2]\n",
    "    return z\n",
    "\n",
    "\n",
    "cpdef ws2doptv(np.ndarray[dtype_t] y, np.ndarray[dtype_t] w, array[double] llas):\n",
    "    \"\"\"Whittaker smoother with normal V-curve optimization of lambda (S).\n",
    "    Args:\n",
    "        y: time-series numpy array\n",
    "        w: weights numpy array\n",
    "        llas: array with lambda values to iterate (S-range)\n",
    "    Returns:\n",
    "        Smoothed time-series array z and optimized lambda (S) value lopt\n",
    "    \"\"\"\n",
    "    cdef array template = array('d', [])\n",
    "    cdef array fits, pens, diff1, lamids, v, z\n",
    "    cdef int m, m1, m2, nl, nl1, lix, i, k\n",
    "    cdef double w_tmp, y_tmp, z_tmp, z2, llastep, f1, f2, p1, p2, l, l1, l2, vmin, lopt\n",
    "\n",
    "    m = y.shape[0]\n",
    "    m1 = m - 1\n",
    "    m2 = m - 2\n",
    "    nl = len(llas)\n",
    "    nl1 = nl - 1\n",
    "    i = 0\n",
    "    k = 0\n",
    "\n",
    "    template = array('d', [])\n",
    "\n",
    "    fits = clone(template, nl, True)\n",
    "    pens = clone(template, nl, True)\n",
    "    z = clone(template, m, False)\n",
    "    diff1 = clone(template, m1, True)\n",
    "    lamids = clone(template, nl1, False)\n",
    "    v = clone(template, nl1, False)\n",
    "\n",
    "    # Compute v-curve\n",
    "    for lix in range(nl):\n",
    "        l = pow(10,llas.data.as_doubles[lix])\n",
    "        z[0:m] = ws2d(y, l, w)\n",
    "        for i in range(m):\n",
    "            w_tmp = w[i]\n",
    "            y_tmp = y[i]\n",
    "            z_tmp = z.data.as_doubles[i]\n",
    "            fits.data.as_doubles[lix] += pow(w_tmp * (y_tmp - z_tmp),2) \n",
    "        fits.data.as_doubles[lix] = log(fits.data.as_doubles[lix])\n",
    "\n",
    "        for i in range(m1):\n",
    "            z_tmp = z.data.as_doubles[i]\n",
    "            z2 = z.data.as_doubles[i+1]\n",
    "            diff1.data.as_doubles[i] = z2 - z_tmp\n",
    "        for i in range(m2):\n",
    "            z_tmp = diff1.data.as_doubles[i]\n",
    "            z2 = diff1.data.as_doubles[i+1]\n",
    "            pens.data.as_doubles[lix] += pow(z2 - z_tmp,2)\n",
    "        pens.data.as_doubles[lix] = log(pens.data.as_doubles[lix])\n",
    "\n",
    "    # Construct v-curve\n",
    "    llastep = llas[1] - llas[0]\n",
    "\n",
    "    for i in range(nl1):\n",
    "        l1 = llas.data.as_doubles[i]\n",
    "        l2 = llas.data.as_doubles[i+1]\n",
    "        f1 = fits.data.as_doubles[i]\n",
    "        f2 = fits.data.as_doubles[i+1]\n",
    "        p1 = pens.data.as_doubles[i]\n",
    "        p2 = pens.data.as_doubles[i+1]\n",
    "        v.data.as_doubles[i] = sqrt(pow(f2 - f1,2) + pow(p2 - p1,2)) / (log(10) * llastep)\n",
    "        lamids.data.as_doubles[i] = (l1+l2) / 2\n",
    "\n",
    "    vmin = v.data.as_doubles[k]\n",
    "    for i in range(1, nl1):\n",
    "        if v.data.as_doubles[i] < vmin:\n",
    "            vmin = v.data.as_doubles[i]\n",
    "            k = i\n",
    "\n",
    "    lopt = pow(10, lamids.data.as_doubles[k])\n",
    "\n",
    "    z[0:m] = ws2d(y, lopt, w)\n",
    "\n",
    "    return z, lopt, v, lamids\n",
    "\n",
    "\n",
    "\n",
    "cpdef ws2doptvp(np.ndarray[dtype_t] y, np.ndarray[dtype_t] w, array[double] llas, double p):\n",
    "    \"\"\"Whittaker smoother with asymmetric V-curve optimization of lambda (S).\n",
    "    Args:\n",
    "        y: time-series numpy array\n",
    "        w: weights numpy array\n",
    "        llas: array with lambda values to iterate (S-range)\n",
    "        p: \"Envelope\" value\n",
    "    Returns:\n",
    "        Smoothed time-series array z and optimized lambda (S) value lopt\n",
    "    \"\"\"\n",
    "    cdef array template = array('d', [])\n",
    "    cdef array fits, pens, diff1, lamids, v, z\n",
    "    cdef int m, m1, m2, nl, nl1, lix, i, j, k\n",
    "    cdef double w_tmp, y_tmp, z_tmp, z2, llastep, fit1, fit2, pen1, pen2, l, l1, l2, vmin, lopt, p1\n",
    "\n",
    "    m = y.shape[0]\n",
    "    m1 = m - 1\n",
    "    m2 = m - 2\n",
    "    nl = len(llas)\n",
    "    nl1 = nl - 1\n",
    "    i = 0\n",
    "    k = 0\n",
    "    j = 0\n",
    "    p1 = 1-p\n",
    "\n",
    "    template = array('d', [])\n",
    "    fits = clone(template, nl, True)\n",
    "    pens = clone(template, nl, True)\n",
    "    z = clone(template, m, True)\n",
    "    znew = clone(template, m, True)\n",
    "    diff1 = clone(template, m1, True)\n",
    "    lamids = clone(template, nl1, False)\n",
    "    v = clone(template, nl1, False)\n",
    "    wa = clone(template, m, False)\n",
    "    ww = clone(template, m, False)\n",
    "\n",
    "    # Compute v-curve\n",
    "    for lix in range(nl):\n",
    "        l = pow(10,llas.data.as_doubles[lix])\n",
    "\n",
    "        for i in range(10):\n",
    "          for j in range(m):\n",
    "            y_tmp = y[j]\n",
    "            z_tmp = z.data.as_doubles[j]\n",
    "            if y_tmp > z_tmp:\n",
    "              wa.data.as_doubles[j] = p\n",
    "            else:\n",
    "              wa.data.as_doubles[j] = p1\n",
    "            ww.data.as_doubles[j] = w[j] * wa.data.as_doubles[j]\n",
    "\n",
    "          znew[0:m] = _ws2d(y, l, ww)\n",
    "          z_tmp = 0.0\n",
    "          j = 0\n",
    "          for j in range(m):\n",
    "            z_tmp += abs(znew.data.as_doubles[j] - z.data.as_doubles[j])\n",
    "\n",
    "          if z_tmp == 0.0:\n",
    "            break\n",
    "\n",
    "          z[0:m]= znew[0:m]\n",
    "\n",
    "        for i in range(m):\n",
    "            w_tmp = w[i]\n",
    "            y_tmp = y[i]\n",
    "            z_tmp = z.data.as_doubles[i]\n",
    "            fits.data.as_doubles[lix] += pow(w_tmp * (y_tmp - z_tmp),2)\n",
    "        fits.data.as_doubles[lix] = log(fits.data.as_doubles[lix])\n",
    "\n",
    "        for i in range(m1):\n",
    "            z_tmp = z.data.as_doubles[i]\n",
    "            z2 = z.data.as_doubles[i+1]\n",
    "            diff1.data.as_doubles[i] = z2 - z_tmp\n",
    "        for i in range(m2):\n",
    "            z_tmp = diff1.data.as_doubles[i]\n",
    "            z2 = diff1.data.as_doubles[i+1]\n",
    "            pens.data.as_doubles[lix] += pow(z2 - z_tmp,2)\n",
    "        pens.data.as_doubles[lix] = log(pens.data.as_doubles[lix])\n",
    "\n",
    "    # Construct v-curve\n",
    "    llastep = llas[1] - llas[0]\n",
    "\n",
    "    for i in range(nl1):\n",
    "        l1 = llas.data.as_doubles[i]\n",
    "        l2 = llas.data.as_doubles[i+1]\n",
    "        fit1 = fits.data.as_doubles[i]\n",
    "        fit2 = fits.data.as_doubles[i+1]\n",
    "        pen1 = pens.data.as_doubles[i]\n",
    "        pen2 = pens.data.as_doubles[i+1]\n",
    "        v.data.as_doubles[i] = sqrt(pow(fit2 - fit1,2) + pow(pen2 - pen1,2)) / (log(10) * llastep)\n",
    "        lamids.data.as_doubles[i] = (l1+l2) / 2\n",
    "\n",
    "    vmin = v.data.as_doubles[k]\n",
    "    for i in range(1, nl1):\n",
    "        if v.data.as_doubles[i] < vmin:\n",
    "            vmin = v.data.as_doubles[i]\n",
    "            k = i\n",
    "\n",
    "    lopt = pow(10, lamids.data.as_doubles[k])\n",
    "\n",
    "    del z\n",
    "    z = clone(template, m, True)\n",
    "\n",
    "    for i in range(10):\n",
    "      for j in range(m):\n",
    "        y_tmp = y[j]\n",
    "        z_tmp = z.data.as_doubles[j]\n",
    "\n",
    "        if y_tmp > z_tmp:\n",
    "          wa.data.as_doubles[j] = p\n",
    "        else:\n",
    "          wa.data.as_doubles[j] = p1\n",
    "        ww.data.as_doubles[j] = w[j] * wa.data.as_doubles[j]\n",
    "\n",
    "      znew[0:m] = _ws2d(y, lopt, ww)\n",
    "      z_tmp = 0.0\n",
    "      j = 0\n",
    "      for j in range(m):\n",
    "        z_tmp += abs(znew.data.as_doubles[j] - z.data.as_doubles[j])\n",
    "\n",
    "      if z_tmp == 0.0:\n",
    "        break\n",
    "\n",
    "      z[0:m]= znew[0:m]\n",
    "\n",
    "    z[0:m] = _ws2d(y, lopt, ww)\n",
    "    \n",
    "    return z, lopt, v, lamids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def daily_interpolation(rawdates: np.ndarray):\n",
    "    \n",
    "    #create daily date range (we add 7 days so that we don't have any problems when shifting the values later)\n",
    "    pd_daily_dates = pd.date_range(start=rawdates[0],end=rawdates[-1] + datetime.timedelta(16))\n",
    "    \n",
    "    #convert pd data range to numpy array \n",
    "    np_daily_dates = np.array(pd_daily_dates.to_pydatetime())\n",
    "    \n",
    "    #convert datetime to date\n",
    "    np_daily_dates = np.array([x.date() for x in np_daily_dates])\n",
    "    \n",
    "    return np_daily_dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def fromstring(x):\n",
    "    \n",
    "    '''Converts string to datetime object'''\n",
    "    \n",
    "    try:\n",
    "        d = datetime.datetime.strptime(x, '%d/%m/%Y').date()\n",
    "    except:\n",
    "        d = datetime.datetime.strptime(x, '%Y-%m-%d').date()\n",
    "        \n",
    "    return d\n",
    "\n",
    "def fromjulian(x):\n",
    "    \n",
    "    '''Converts julian to datetime object'''\n",
    "\n",
    "    return datetime.datetime.strptime(x, '%Y%j').date()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract time series functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def ndvi_select_satellite(sat: str, MD: tuple):\n",
    "    \n",
    "    (MOD, MYD, MXD) = MD\n",
    "    \n",
    "    if (sat == 'MYD'):\n",
    "        df = MYD\n",
    "    elif (sat == 'MOD'):\n",
    "        df = MOD\n",
    "    else:\n",
    "        df = MXD\n",
    "        \n",
    "    return df\n",
    "\n",
    "def ndvi_extract_ts(df: pd.DataFrame, location: str, date_begin: int, date_end: int):\n",
    "    \n",
    "    y = df['NDVI'].loc[location].values\n",
    "\n",
    "    dts = df['Date'].loc[location].values\n",
    "\n",
    "    c_dts = df['Composite_date'].loc[location].values\n",
    "\n",
    "    same_c = df['Same_composite'].loc[location].values\n",
    "    \n",
    "    #Crop to date range\n",
    "    date_range = np.all([dts>=datetime.date(date_begin,1,1), dts<=datetime.date(date_end,12,31)], axis=0)\n",
    "    y = y[date_range]\n",
    "    dts = dts[date_range]\n",
    "    c_dts = c_dts[date_range]\n",
    "    same_c = same_c[date_range]\n",
    "    \n",
    "    return (y, dts, c_dts, same_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smoothing functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def ndvi_smoothing_A0(y: np.ndarray, dts: np.ndarray, pvalue: float, nopval: bool, lrange: np.ndarray, nd: int):\n",
    "    \n",
    "     # create weights\n",
    "    w = np.array((y!=nd)*1,dtype='double')\n",
    "\n",
    "    # apply whittaker filter with V-curve\n",
    "    if (nopval):\n",
    "        z, lopt, vcurve, l = ws2doptv(y, w, array.array('d',lrange))\n",
    "    else:\n",
    "        z, lopt, vcurve, l = ws2doptvp(y, w, array.array('d',lrange), pvalue)\n",
    "    \n",
    "    return (z, lopt, vcurve, l)\n",
    "\n",
    "\n",
    "    \n",
    "def ndvi_smoothing_A1(y: np.ndarray, dts: np.ndarray, pvalue: float, nopval: bool, lrange: np.ndarray, nd: int):\n",
    "\n",
    "    # create weights\n",
    "    w = np.array((y!=nd)*1,dtype='double')\n",
    "\n",
    "    # apply whittaker filter with V-curve\n",
    "    if (nopval):\n",
    "        z, lopt, vcurve, l = ws2doptv(y, w, array.array('d',lrange))\n",
    "    else:\n",
    "        z, lopt, vcurve, l = ws2doptvp(y, w, array.array('d',lrange), pvalue)\n",
    "\n",
    "    # Temporal interpolation\n",
    "    daily = daily_interpolation(dts)\n",
    "    dvec = np.full(len(daily), nd, dtype='double')\n",
    "\n",
    "    # shift observations to midpoint of acquisition (these positions are set to 0 instead of nodata)\n",
    "    mid_point = 16/2\n",
    "    for d in dts:\n",
    "        dl = daily.tolist()\n",
    "        dvec[dl.index(d + datetime.timedelta(mid_point))] = 0\n",
    "\n",
    "    # place des filtered values in the midpoints\n",
    "    dvec[ dvec != nd ] = z\n",
    "\n",
    "    # recreate weights\n",
    "    w = np.array((dvec != nd) * 1,dtype='double')\n",
    "\n",
    "    #refilter with low lanmbda\n",
    "    z =  ws2d(dvec, 0.0001, w)\n",
    "    \n",
    "    return (z, lopt, vcurve, l, daily)\n",
    "\n",
    "\n",
    "def ndvi_smoothing_A2(y: np.ndarray, dts: np.ndarray, pvalue: float, nopval: bool, lrange: np.ndarray, nd: int):\n",
    "    \n",
    "    # Temporal interpolation\n",
    "    daily = daily_interpolation(dts)\n",
    "    dvec = np.full(len(daily), nd, dtype='double')\n",
    "\n",
    "    mid_point = 16/2\n",
    "    for d in dts:\n",
    "        dl = daily.tolist()\n",
    "        dvec[dl.index(d + datetime.timedelta(mid_point))] = 0\n",
    "\n",
    "    # place des values in the right place\n",
    "    dvec[ dvec != nd ] = y\n",
    "\n",
    "    # create weights\n",
    "    w = np.array((dvec!=nd)*1,dtype='double')\n",
    "\n",
    "    # apply whittaker filter with V-curve\n",
    "    if (nopval):\n",
    "        z, lopt, vcurve, l = ws2doptv(dvec, w, array.array('d',lrange))\n",
    "    else:\n",
    "        z, lopt, vcurve, l = ws2doptvp(dvec, w, array.array('d',lrange), pvalue)        \n",
    "        \n",
    "    return (z, lopt, vcurve, l)\n",
    "\n",
    "\n",
    "def ndvi_smoothing_A3(y: np.ndarray, dts: np.ndarray, same_c: np.ndarray, c_dts: np.ndarray, pvalue: float, nopval: bool, lrange, nd: int):\n",
    "    \n",
    "    #remove the values with the same day composite\n",
    "    y_unique = np.delete(y, np.argwhere(same_c))\n",
    "    c_dts = np.delete(c_dts, np.argwhere(same_c))\n",
    "\n",
    "    # Temporal interpolation\n",
    "    daily = daily_interpolation(dts)\n",
    "    dvec = np.full(len(daily), nd, dtype='double')\n",
    "\n",
    "\n",
    "    for d in c_dts:\n",
    "        dl = daily.tolist()\n",
    "        dvec[dl.index(d)] = 0\n",
    "\n",
    "    # place des values in the right place\n",
    "    dvec[ dvec != nd ] = y_unique\n",
    "\n",
    "    # create weights\n",
    "    w = np.array((dvec!=nd)*1,dtype='double')\n",
    "\n",
    "    # apply whittaker filter with V-curve\n",
    "    if (nopval):\n",
    "        z, lopt, vcurve, l = ws2doptv(dvec, w, array.array('d',lrange))\n",
    "    else:\n",
    "        z, lopt, vcurve, l = ws2doptvp(dvec, w, array.array('d',lrange), pvalue)\n",
    "        \n",
    "    return (z, lopt, vcurve, l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print and Plot functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def ndvi_print_info(location: str, latlon: dict, lagCorr1: float):\n",
    "    \n",
    "    print('\\033[1m' + 'Selected Point: ', location, '\\033[0m')\n",
    "    print('(Lat, Lon) =', latlon[location] )\n",
    "    print('\\n')\n",
    "    \n",
    "    print('LagCorr: ', round(lagCorr1,3))\n",
    "    print('\\n')\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "def ndvi_plot_all(A0: bool, A1: bool, A2: bool, A3: bool,\n",
    "                  dts: np.ndarray, y: np.ndarray, daily: np.ndarray, \n",
    "                  z0: np.ndarray, z1: np.ndarray, z2: np.ndarray, z3: np.ndarray, \n",
    "                  nd: int, \n",
    "                  yauto: bool, ylimits: tuple):\n",
    "    \n",
    "    #replace nd by nan\n",
    "    ynan = y.copy()\n",
    "    ynan[ynan == nd] = np.nan\n",
    "    \n",
    "    plt.figure(figsize=(20,5))\n",
    "    \n",
    "    A = [A0, A1, A2, A3]\n",
    "    xA = [dts, daily, daily, daily]\n",
    "    z = [z0, z1, z2, z3]\n",
    "    col = ['b', 'g', 'r', 'orange']\n",
    "    label = ['A0', 'A1', 'A2', 'A3']\n",
    "    leg = []\n",
    "    \n",
    "    for i,a in enumerate(A):\n",
    "        if a:\n",
    "            plt.plot(xA[i], z[i], color = col[i], alpha = 0.5)\n",
    "            leg.append(label[i])\n",
    "\n",
    "    leg.append('raw values')\n",
    "    plt.plot(dts, ynan, color = 'grey', marker = 'o', alpha = 0.5)\n",
    "    \n",
    "    if not(yauto):\n",
    "        plt.ylim(ylimits)\n",
    "    \n",
    "    plt.xlabel('Date', fontsize=15)\n",
    "    plt.ylabel('NDVI', fontsize=15)\n",
    "    plt.legend(leg, fontsize=17, loc = 'lower right')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "    \n",
    "def ndvi_plot_vcurve(A0: bool, A1: bool, A2: bool, A3: bool,\n",
    "                     l0: np.ndarray, l1: np.ndarray, l2: np.ndarray, l3: np.ndarray, \n",
    "                     vcurve0: np.ndarray, vcurve1: np.ndarray, vcurve2: np.ndarray, vcurve3: np.ndarray, \n",
    "                     lopt0: float, lopt1: float, lopt2: float, lopt3: float):\n",
    "    \n",
    "    plt.figure(figsize=(20,10))\n",
    "    \n",
    "    A = [A0, A1, A2, A3]\n",
    "    xA = [l0, l1, l2, l3]\n",
    "    v = [vcurve0, vcurve1, vcurve2, vcurve3]\n",
    "    lopt = [lopt0, lopt1, lopt2, lopt3]\n",
    "    col = ['blue', 'green', 'red', 'orange']\n",
    "    leg = []\n",
    "    \n",
    "    for i,a in enumerate(A):\n",
    "        if a:\n",
    "            plt.plot(xA[i], v[i], color = col[i], alpha = 0.5, marker = 'o')\n",
    "            leg.append('lopt: ' + str(round(np.log10(lopt[i]),2)))\n",
    "            \n",
    "    for i,a in enumerate(A):\n",
    "        if a:\n",
    "            plt.axvline(x = np.log10(lopt[i]), ls = '--', color = col[i])\n",
    "\n",
    "    plt.xlabel('log10(l)', fontsize=15)\n",
    "    plt.ylabel('V', fontsize=15)\n",
    "    plt.title('V-curves', fontsize=23)\n",
    "    plt.legend(leg, fontsize=17, loc = 'lower right')\n",
    "    plt.show()\n",
    "            \n",
    "    \n",
    "def ndvi_plot_year(A0: bool, A1: bool, A2: bool, A3: bool,\n",
    "                   raw: bool, raw_midpoint: bool, raw_realpoint: bool,\n",
    "                   year: float, month: int, \n",
    "                   lta_mean: dict, lta_std: dict, \n",
    "                   daily: np.ndarray, \n",
    "                   z0: np.ndarray, z1: np.ndarray, z2: np.ndarray, z3: np.ndarray, \n",
    "                   z0_lta: np.ndarray, z1_lta: np.ndarray, z2_lta: np.ndarray, \n",
    "                   dts: np.ndarray,  c_dts:np.ndarray, y: np.ndarray, \n",
    "                   nd: int, \n",
    "                   yauto: bool, ylimits: tuple,\n",
    "                   satellite: str):\n",
    "\n",
    "    \n",
    "    #replace nd by nan\n",
    "    ynan = y.copy()\n",
    "    ynan[ynan == nd] = np.nan\n",
    "    \n",
    "    fig, axs = plt.subplots(1,2, figsize=(20, 10))\n",
    "\n",
    "    # cropping daily data to year\n",
    "    year_index = np.all([daily>=datetime.date(year,month,1), daily<datetime.date(year+1,month,1)], axis=0)\n",
    "\n",
    "    year_daily = daily[year_index]\n",
    "    year_z1 = np.array(z1)[year_index]\n",
    "    year_z2 = np.array(z2)[year_index]\n",
    "    year_z3 = np.array(z3)[year_index]\n",
    "    year_z1_lta = np.array(z1_lta)[year_index]\n",
    "    year_z2_lta = np.array(z2_lta)[year_index]\n",
    "\n",
    "    # cropping raw data to year\n",
    "    year_index2 = np.all([dts>=datetime.date(year,month,1), dts<datetime.date(year+1,month,1)], axis=0)\n",
    "\n",
    "    year_dts = dts[year_index2]\n",
    "    year_ynan = ynan[year_index2]\n",
    "    year_z0 = np.array(z0)[year_index2]\n",
    "    year_z0_lta = np.array(z0_lta)[year_index2]\n",
    "    year_lta_mean = []\n",
    "    year_lta_std = []\n",
    "    for dt in year_dts:\n",
    "        year_lta_mean.append(lta_mean[(dt.day, dt.month)])\n",
    "        year_lta_std.append(lta_std[(dt.day, dt.month)])\n",
    "        \n",
    "    #creating year_dts_midpoint and year_dts_realpoint\n",
    "    midpoint = 16/2\n",
    "    year_dts_midpoint = np.array([dt + timedelta(days=midpoint) for dt in year_dts])\n",
    "    \n",
    "    #cropping composite date\n",
    "    year_index3 = np.all([c_dts>=datetime.date(year,month,1), c_dts<datetime.date(year+1,month,1)], axis=0)\n",
    "    year_ynan_c = ynan[year_index3]\n",
    "    year_dts_realpoint = c_dts[year_index3]\n",
    "        \n",
    "    \n",
    "    \n",
    "    #PLOTTING THE YEARLY DATA    \n",
    "    A = [A0, A1, A2, A3]\n",
    "    xA = [year_dts, year_daily, year_daily, year_daily]\n",
    "    z = [year_z0, year_z1, year_z2, year_z3]\n",
    "    col = ['b', 'g', 'r', 'orange']\n",
    "    label = ['A0', 'A1', 'A2', 'A3']\n",
    "    leg = []\n",
    "    \n",
    "    for i,a in enumerate(A):\n",
    "        if a:\n",
    "            axs[0].plot(xA[i], z[i], color = col[i], alpha = 0.5)\n",
    "            leg.append(label[i])\n",
    "\n",
    "    #RAW values\n",
    "    R = [raw, raw_midpoint, raw_realpoint]\n",
    "    xR = [year_dts, year_dts_midpoint, year_dts_realpoint]\n",
    "    ylist = [year_ynan, year_ynan, year_ynan_c]\n",
    "    label = ['raw', 'raw midpoint', 'raw composite date']\n",
    "    \n",
    "    for i,r in enumerate(R):\n",
    "        if r:\n",
    "            axs[0].plot(xR[i], ylist[i], color = 'grey', marker = 'o', alpha = 0.5)\n",
    "            leg.append(label[i])\n",
    "            \n",
    "    \n",
    "    #Color MOD and MYD points if MXD data\n",
    "    if (satellite == 'MXD'):\n",
    "        len_date = len(dts[dts < datetime.date(year,month,1)])\n",
    "        if (len_date % 2) == 0:\n",
    "            first = 'MYD'\n",
    "            second = 'MOD'\n",
    "        else:\n",
    "            first = 'MOD'\n",
    "            second = 'MYD'\n",
    "            \n",
    "        for i,r in enumerate(R):\n",
    "            if r:      \n",
    "                axs[0].plot(xR[i][0::2], ylist[i][0::2], 'ks')\n",
    "                axs[0].plot(xR[i][1::2], ylist[i][1::2], 'k^')\n",
    "        leg.append(first)\n",
    "        leg.append(second)\n",
    "    \n",
    "    if not(yauto):\n",
    "        axs[0].set_ylim(ylimits)\n",
    "\n",
    "    axs[0].set_xlabel('Date', fontsize=15)\n",
    "    axs[0].set_ylabel('NDVI', fontsize=15)\n",
    "    axs[0].legend(leg, fontsize=15, loc = 'lower right')\n",
    "    axs[0].set_title('Year ' + str(year), fontsize = 20)\n",
    "    \n",
    "    \n",
    "    #PLOTTING THE LTA\n",
    "    \n",
    "    #shifting to midpoint\n",
    "    A = [A0, A1, A2]\n",
    "    xA = [year_dts, year_daily, year_daily]\n",
    "    z = [year_z0_lta, year_z1_lta, year_z2_lta]\n",
    "    col = ['b', 'g', 'r']\n",
    "    label = ['A0 (lta)', 'A1 (lta)', 'A2 (lta)']\n",
    "    leg = []\n",
    "    \n",
    "    for i,a in enumerate(A):\n",
    "        if a:\n",
    "            axs[1].plot(xA[i], z[i], color = col[i], alpha = 0.5)\n",
    "            leg.append(label[i])\n",
    "\n",
    "    axs[1].plot(year_dts_midpoint, year_lta_mean, color = 'grey', marker = 'o', alpha = 0.5)\n",
    "    axs[1].fill_between(year_dts_midpoint, np.array(year_lta_mean)-2*np.array(year_lta_std), np.array(year_lta_mean)+2*np.array(year_lta_std), color = 'whitesmoke')\n",
    "    axs[1].fill_between(year_dts_midpoint, np.array(year_lta_mean)-np.array(year_lta_std), np.array(year_lta_mean)+np.array(year_lta_std), color = 'lightgrey')\n",
    "    \n",
    "    leg.append('raw values')\n",
    "    leg.append('+- 2sd')\n",
    "    leg.append('+-1sd')\n",
    "    \n",
    "    if not(yauto):\n",
    "        axs[1].set_ylim(ylimits)\n",
    "\n",
    "    axs[1].set_xlabel('Date', fontsize=15)\n",
    "    axs[1].set_ylabel('LTA NDVI', fontsize=15)\n",
    "    axs[1].legend(leg, fontsize=15, loc = 'lower right')\n",
    "    axs[1].set_title('Long Term Averages ' + str(year), fontsize = 20)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LTA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def ndvi_lta_dict(y: np.ndarray, dts: np.ndarray, nd: float):\n",
    "    \n",
    "    dict_index = set((dt.day, dt.month) for dt in dts)\n",
    "\n",
    "    # initialize the dicts\n",
    "    lta_mean = {}\n",
    "    lta_std = {}\n",
    "\n",
    "    for dt in dict_index:\n",
    "\n",
    "        lta_date = [] \n",
    "\n",
    "        for ix, date_sel in enumerate(dts):\n",
    "            \n",
    "            if (date_sel.month == dt[1]) & (date_sel.day == dt[0]):     \n",
    "                if (y[ix] != nd):\n",
    "                    lta_date.append(y[ix])\n",
    "\n",
    "\n",
    "        # add to LTA dict\n",
    "        lta_mean[dt] = np.mean(lta_date)\n",
    "        lta_std[dt] = np.std(lta_date)\n",
    "        \n",
    "    return lta_mean, lta_std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Main function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def ndvi_main(A0: bool, A1: bool, A2: bool, A3: bool,\n",
    "              raw: bool, raw_midpoint: bool, raw_realpoint: bool,\n",
    "              MD: tuple, location: str, latlon: dict, satellite: str, \n",
    "              pvalue:float, nopval: bool, \n",
    "              lrange1: tuple, lrange2: tuple, step: float, \n",
    "              nd: int, \n",
    "              date_begin: int, date_end: int, \n",
    "              year: int, month: int,\n",
    "              yauto: bool, ylimits: tuple):\n",
    "    \n",
    "    lrange1 = np.arange(lrange1[0],lrange1[1], step)\n",
    "    lrange2 = np.arange(lrange2[0],lrange2[1], step)\n",
    "    \n",
    "    #extract time series\n",
    "    df = ndvi_select_satellite(satellite, MD)\n",
    "    (y, dts, c_dts, same_c) = ndvi_extract_ts(df, location, date_begin, date_end)\n",
    "    \n",
    "    #LTA\n",
    "    lta_mean, lta_std = ndvi_lta_dict(y, dts, nd)\n",
    "    \n",
    "    #smoothing\n",
    "    (z0, lopt0, vcurve0, l0) = ndvi_smoothing_A0(y, dts, pvalue, nopval, lrange1, nd)\n",
    "    (z1, lopt1, vcurve1, l1, daily) = ndvi_smoothing_A1(y, dts, pvalue, nopval, lrange1, nd)\n",
    "    (z2, lopt2, vcurve2, l2) = ndvi_smoothing_A2(y, dts, pvalue, nopval, lrange2, nd)\n",
    "    (z3, lopt3, vcurve3, l3) = ndvi_smoothing_A3(y, dts, same_c, c_dts, pvalue, nopval, lrange2, nd)\n",
    "    \n",
    "    #smoothing LTA data\n",
    "    # no daily date for LTA\n",
    "    lta_mean_list = []\n",
    "    for dt in dts:\n",
    "        lta_mean_list.append(lta_mean[(dt.day, dt.month)])\n",
    "    (z0_lta, lopt0_lta, vcurve0_lta, l0_lta) = ndvi_smoothing_A0(np.array(lta_mean_list), dts, pvalue, nopval, lrange1, nd)\n",
    "    (z1_lta, lopt1_lta, vcurve1_lta, l1_lta, daily_lta) = ndvi_smoothing_A1(np.array(lta_mean_list), dts, pvalue, nopval, lrange1, nd)\n",
    "    (z2_lta, lopt2_lta, vcurve2_lta, l2_lta) = ndvi_smoothing_A2(np.array(lta_mean_list), dts, pvalue, nopval, lrange2, nd)\n",
    "    \n",
    "    #lagCorr\n",
    "    lagCorr1 = lag1corr(np.array(y[0:len(y)-1]), np.array(y[1:]), nd)\n",
    "    \n",
    "    #Prints and plots\n",
    "    ndvi_print_info(location, latlon, lagCorr1)\n",
    "    ndvi_plot_all(A0, A1, A2, A3, dts, y, daily, z0, z1, z2, z3, nd, yauto, ylimits)\n",
    "    ndvi_plot_year(A0, A1, A2, A3, raw, raw_midpoint, raw_realpoint, year, month, lta_mean, lta_std, daily, z0, z1, z2, z3, z0_lta, z1_lta, z2_lta, dts,  c_dts, y, nd, yauto, ylimits, satellite)\n",
    "    ndvi_plot_vcurve(A0, A1, A2, A3, l0, l1, l2, l3, vcurve0, vcurve1, vcurve2, vcurve3, lopt0, lopt1, lopt2, lopt3)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_nd = -3000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# loading ndvi_MOD data from csv\n",
    "ndvi_MOD = pd.read_csv(\n",
    "    'data/MOD13A2-MOD13A2-006-results.csv', \n",
    "    index_col=0, \n",
    "    usecols = ['ID', \n",
    "               'Date', \n",
    "               'MOD13A2_006__1_km_16_days_NDVI', \n",
    "               'MOD13A2_006__1_km_16_days_composite_day_of_the_year'],\n",
    "    dtype = {'MOD13A2_006__1_km_16_days_composite_day_of_the_year': int})\n",
    "\n",
    "\n",
    "#renaming the columns\n",
    "ndvi_MOD = ndvi_MOD.rename(columns={\"MOD13A2_006__1_km_16_days_NDVI\": \"NDVI\", \n",
    "                        \"MOD13A2_006__1_km_16_days_composite_day_of_the_year\": \"Composite_date\"})\n",
    "\n",
    "\n",
    "# Convert string Date to datetime.date\n",
    "ndvi_MOD['Date'] = ndvi_MOD['Date'].apply(fromstring)\n",
    "\n",
    "\n",
    "# Convert composite_date from julian to datetime\n",
    "# Add a True/False column to keep track of the values with the same composite_date\n",
    "\n",
    "compo = []\n",
    "same_compo = []\n",
    "\n",
    "for i in range(len(ndvi_MOD)):\n",
    "    \n",
    "    d = ndvi_MOD['Date'][i]\n",
    "    c = ndvi_MOD['Composite_date'][i]\n",
    "    \n",
    "    if c != -1:    \n",
    "        if (d.month == 12) and (c<20):\n",
    "            compo.append(fromjulian(str(d.year+1)+str(c)))\n",
    "        else:\n",
    "            compo.append(fromjulian(str(d.year)+str(c)))\n",
    "    else:\n",
    "        # nodata so we don't care about the date\n",
    "        compo.append(d)\n",
    "    \n",
    "    if (i==0):\n",
    "        same_compo.append(False)\n",
    "    elif (compo[i]==compo[i-1]):\n",
    "        same_compo.append(True)\n",
    "    else:\n",
    "        same_compo.append(False)\n",
    "            \n",
    "ndvi_MOD['Composite_date'] = compo   \n",
    "ndvi_MOD['Same_composite'] = same_compo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# loading MYD data from csv\n",
    "ndvi_MYD = pd.read_csv(\n",
    "    'data/MYD13A2-MYD13A2-006-results.csv', \n",
    "    index_col=0, \n",
    "    usecols = ['ID', \n",
    "               'Date', \n",
    "               'MYD13A2_006__1_km_16_days_NDVI', \n",
    "               'MYD13A2_006__1_km_16_days_composite_day_of_the_year'],\n",
    "    dtype = {'MYD13A2_006__1_km_16_days_composite_day_of_the_year': int})\n",
    "\n",
    "\n",
    "#renaming the columns\n",
    "ndvi_MYD = ndvi_MYD.rename(columns={\"MYD13A2_006__1_km_16_days_NDVI\": \"NDVI\", \n",
    "                        \"MYD13A2_006__1_km_16_days_composite_day_of_the_year\": \"Composite_date\"})\n",
    "\n",
    "\n",
    "# Convert string Date to datetime.date\n",
    "ndvi_MYD['Date'] = ndvi_MYD['Date'].apply(fromstring)\n",
    "\n",
    "\n",
    "# Convert composite_date from julian to datetime\n",
    "# Add a True/False column to keep track of the values with the same composite_date\n",
    "\n",
    "compo = []\n",
    "same_compo = []\n",
    "\n",
    "for i in range(len(ndvi_MYD)):\n",
    "    \n",
    "    d = ndvi_MYD['Date'][i]\n",
    "    c = ndvi_MYD['Composite_date'][i]\n",
    "    \n",
    "    if c != -1:    \n",
    "        if (d.month == 12) and (c<20):\n",
    "            compo.append(fromjulian(str(d.year+1)+str(c)))\n",
    "        else:\n",
    "            compo.append(fromjulian(str(d.year)+str(c)))\n",
    "    else:\n",
    "        # nodata so we don't care about the date\n",
    "        compo.append(d)\n",
    "    \n",
    "    if (i==0):\n",
    "        same_compo.append(False)\n",
    "    elif (compo[i]==compo[i-1]):\n",
    "        same_compo.append(True)\n",
    "    else:\n",
    "        same_compo.append(False)\n",
    "            \n",
    "ndvi_MYD['Composite_date'] = compo   \n",
    "ndvi_MYD['Same_composite'] = same_compo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "#creating MXD\n",
    "\n",
    "#changing the index to numbers\n",
    "number_index1 = pd.Index(range(0,2*len(ndvi_MYD),2))\n",
    "number_index2 = pd.Index(range(1,2*len(ndvi_MYD)+1,2))\n",
    "\n",
    "ndvi_MYD_nb = ndvi_MYD.set_index(number_index1)\n",
    "ndvi_MYD_nb['ID'] = ndvi_MYD.index\n",
    "ndvi_MOD_nb = ndvi_MOD.set_index(number_index2)\n",
    "ndvi_MOD_nb['ID'] = ndvi_MOD.index\n",
    "\n",
    "\n",
    "#concatenating and resetting ID as index\n",
    "ndvi_MXD = pd.concat([ndvi_MYD_nb, ndvi_MOD_nb]).sort_index()\n",
    "ndvi_MXD = ndvi_MXD.set_index('ID')\n",
    "\n",
    "#re-run same_composite\n",
    "same_compo = []\n",
    "for i in range(len(ndvi_MXD)):\n",
    "    if (i==0 or i==1):\n",
    "        same_compo.append(False)\n",
    "    elif (ndvi_MXD['Composite_date'][i]==ndvi_MXD['Composite_date'][i-1] or ndvi_MXD['Composite_date'][i]==ndvi_MXD['Composite_date'][i-2]):\n",
    "        same_compo.append(True)\n",
    "    else:\n",
    "        same_compo.append(False)\n",
    "              \n",
    "ndvi_MXD['Same_composite'] = same_compo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "#Creating lat lon dictionnary\n",
    "\n",
    "latlon = {}\n",
    "\n",
    "pd_latlon = pd.read_csv(\n",
    "    'data/MYD13A2-MYD13A2-006-results.csv', \n",
    "    index_col=0, \n",
    "    usecols = ['ID', \n",
    "               'Latitude', \n",
    "               'Longitude'])\n",
    "\n",
    "dict_index = set(pd_latlon.index)\n",
    "\n",
    "for location in dict_index:\n",
    "    latlon[location] = (round(pd_latlon.loc[location]['Latitude'][0],3), round(pd_latlon.loc[location]['Longitude'][0],3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}