PCA_baseline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.- Data preparation:\n",
    "\n",
    "### Aux function to grow nan areas in space for DataArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "\n",
    "def buffer_nans(da, kn):\n",
    "    k = np.zeros((kn,kn))\n",
    "    k[kn//2,kn//2] = 1\n",
    "\n",
    "    arr = da.values\n",
    "    mask = np.ones(arr.shape).astype(np.float32)\n",
    "\n",
    "    for i in range(arr.shape[0]):\n",
    "        mask[i,:,:] = signal.convolve2d(arr[i,:,:], k, boundary='fill', mode='same')\n",
    "\n",
    "    return da.where(~np.isnan(mask))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries and load dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "ds = xr.open_dataset(\"Murrumbidgee_near_Bundure__MUR_B3.nc\")\n",
    "\n",
    "ds = ds.isel(x=slice(400,800), y=slice(0,400))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute masking using the blue band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blue = ds.nbart_blue.astype(np.float32) / 1e4\n",
    "\n",
    "# 1.- Filter reflectances greater than 0.5 \n",
    "blue = blue.where(blue<.5)\n",
    "\n",
    "# 2.- Filter reflectances with difference to lower quartile larger than 0.05 => (0.07)\n",
    "blue = blue.where((blue - blue.quantile(0.25, dim='time'))<.07)\n",
    "\n",
    "# 3.- Grow a 5x5 buffer around NaN pixels \n",
    "blue = buffer_nans(blue, 5)\n",
    "\n",
    "# 4.- Discard frames with more than 25% missing pixels\n",
    "blue = blue.isel(time=(blue.count(dim=('x','y'))/(400*400))>.25)\n",
    "\n",
    "blue.isel(time=3).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply mask to all the other bands and save for PCA training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for band_name in ds:\n",
    "    band = ds[band_name].astype(np.float32) / 1e4\n",
    "    \n",
    "    # 1. Apply blue mask\n",
    "    band = band.sel(time=blue.time).where(~np.isnan(blue))\n",
    "\n",
    "    # 2.- Interpolate NaNs over time linearly\n",
    "    band = band.interpolate_na(dim='time')\n",
    "\n",
    "    # 3.- Interpolate NaNs at the start and end using nearest neighbor\n",
    "    band = band.interpolate_na(dim='time', method='nearest', fill_value='extrapolate')\n",
    "\n",
    "    # 4.- Apply median rolling filter along time (window=3)\n",
    "    band = band.rolling(time=3, min_periods=1).median()\n",
    "    \n",
    "    np.save(band_name, band)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA decomposition\n",
    "\n",
    "### Load libraries and create data stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(735, 400, 400)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "stack = np.empty((0,400,400))\n",
    "\n",
    "for fname in [\"nbart_red.npy\",\"nbart_green.npy\",\"nbart_blue.npy\",\n",
    "              \"nbart_nir_1.npy\",\"nbart_nir_2.npy\",\"nbart_swir_2.npy\",\n",
    "              \"nbart_swir_3.npy\"]:\n",
    "    \n",
    "    band = np.load(fname)\n",
    "    stack = np.append(stack, band, axis=0)\n",
    "    \n",
    "stack = stack.reshape(stack.shape[0], -1)\n",
    "    \n",
    "stack.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute PCA and plot explained variance with the number of components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8XGV97/HPd19yJSFAQowksANGSkQEmhNRqFLRCl6g1V7AKoJYSgt4bwu2p1Z7OdRqKz1wShGjohaOeKFoseCxBESoECBAuAQDBBMIJESTHZLs28zv/LGeSSY7c1nZMHvN3vm+X8xr1nrW7Tckmd88z7PW8ygiMDMza6aj6ADMzGxscMIwM7NcnDDMzCwXJwwzM8vFCcPMzHJxwjAzs1ycMMzMLBcnDDMzy8UJw8zMcukqOoCX0syZM6Onp6foMMzMxox77rnn+YiYlWffcZUwenp6WLZsWdFhmJmNGZKeyruvm6TMzCwXJwwzM8vFCcPMzHJxwjAzs1ycMMzMLJeWJQxJSyStl7SiznZJ+mdJqyQ9IOnYqm0nS1qZtl3UqhjNzCy/VtYwvgKc3GD7KcCC9DoX+BcASZ3A5Wn7QuAMSQtbGKeZmeXQsucwIuI2ST0NdjkNuDqyOWL/W9IMSXOAHmBVRDwBIOnatO/DrYrVxo+IIAJKEZQjKJehHEEpgihDkG0PsvJsOSvYtWznuSr7DD+OVFau3h7Z9iyWneXliHTO7Lhy5bxV5bss79gnK6NSVq6Kbcdn3vHpd1mv3h51t+06RfPO7cP2r3PcsEs3uE7j7Qw7387r1inf0/ibbN95fO3yHdtrF9fdv/pae3bMnl1oysQuznvjYfVP+BIp8sG9g4A1VetrU1mt8tfWO4mkc8lqKBx88MEvfZTjRKkcDAyV6R8q0T9U3mW5f6hM/2CZgVKZ/sHq7WUGS2WGysHQjvegVE7LaX2oXGufYLBUzt7L2TGDqbwc2ZdhOS2XysO/5LPtWXmk8kb7k86Z7Wc2nkm7l83cZ+K4Txg1PjbRoLymiLgSuBJg0aJF4+7rIiJ4oX+ITdsG+eW2gR3vm7cP8sutg/T2DbJtYIgX+kts6x9i68AQ2wZKbO3P3l/oH2L7QImhl/ibtLtTdHaI7o4OOjtFV0cHXR2iq1Ppfed6Z0cH3R3Z/hO7OujsEB0SHYLODiGJTomODlK5UjlZuURHx879K/vscnw6dse5RDpm120i+weXrae/bKm8ukyiav/djyOVdwhE4+M6OrJ90n9V16k+btflDlW+GHaWDz8uhU7aa9f1GuWVf1g7v3CGHbNju4at174Gw7fnPG7Y5etub3Y+dvusLzL+enEyfP/aW+rtX+uazc7VropMGGuBeVXrc4FngAl1yselzdsHeWLDC6zeuJV1m/t4bnMfz/X282xvH8/19rFhS3/DL/upEzqZOrErvTqZMqGLA6ZOYN7+U3Zsm9zdyaTuTiZ0dTCxqyO9d1Yt71o2MS1P6Oqgu5IMKgmisyN9mY2tv+hm9uIVmTBuAC5IfRSvBTZHxDpJG4AFkuYDTwOnA+8pMM6XRETw1MZt3L92E8vXbOKhp3t54vkXeP6FgV32mzapi9nTJ/Gy6ZM47LCZzJo2kf2ndjNjygT2mzKBGVO62W9Ktr7v5G66O31ntJmNjpYlDEnXACcCMyWtBT4FdANExBXAjcDbgFXANuDstG1I0gXATUAnsCQiHmpVnK20edsgP161gaUrN3DrYxvYsKUfgMndnSx8+XRO+pXZHDprKofO2of5M6cwZ9/JTJ04rsaDNLNxpJV3SZ3RZHsA59fZdiNZQhlzBobK/OiR5/jOfU+zdOV6BkvBvpO7+bUFMzn+FTM5et4MFhy4D12uGZjZGOOfsy+RgaEy31y2hn9Z+jhPb9rOgdMmctbrezj5yDkcPW8GnR1u8zezsc0J4yXwwNpN/Ml1D7DyuS0cPW8GnzntVZx4+IFOEmY2rjhhvAgRwVfuWM3f/McjzNxnAl88cxFvPuJA30FkZuOSE8YIlcvBn1//INfctYa3LJzN537nNew7ubvosMzMWsYJYwQidiaL83/9MD7+lsPpcPOTmY1zThgj8PmbH+Oau9Zwwa+/gk+89fCiwzEzGxW+t3MP3fLoei67ZRW/t2geH/+NVxYdjpnZqHHC2APre/v46DeXc8Sc6Xz6tFe5c9vM9ipOGHvg7258hG39JS57zzFM6u4sOhwzs1HlhJHTT5/YyPXLn+EP33goh83ap+hwzMxGnRNGDhHBZ77/MAfNmMwfn/iKosMxMyuEE0YOt696noee6eXDJy1g8gQ3RZnZ3skJI4cv/2Q1s6ZN5LRjXl50KGZmhXHCaOK53j6WrlzP7/zqXCZ2uXZhZnsvJ4wmvnPv05QDfvtX5xYdiplZoZwwmvj2vWtZdMh+HOo7o8xsL+eE0cBTG7eyav0LvP2oOUWHYmZWuFwJQ9IJks5Oy7PSfNvj3tKVGwD49cMPLDgSM7PiNU0Ykj4F/BlwcSrqBr7eyqDaxdKV65k/cyo9M6cWHYqZWeHy1DB+CzgV2AoQEc8A01oZVDsolYNlq3/J6w47oOhQzMzaQp6EMRARAQSApL3i5/Yj63rZ0j/Ea+fvX3QoZmZtIU/C+KakfwVmSPoD4P8BX2xtWMW768lfAPA/epwwzMwgxwRKEfE5SW8BeoHDgb+MiB+2PLKC3bdmE3P2ncTLZ0wuOhQzs7bQNGGkO6J+XEkSkiZL6omI1a0OrkgPPb2ZVx+0b9FhmJm1jTxNUtcB5ar1Uiobt7b0DfLE81s50gnDzGyHPAmjKyIGKitpeULrQireI+u2AHDkQdMLjsTMrH3kSRgbJJ1aWZF0GvB860Iq3iPregFYOMc1DDOziqZ9GMB5wDckXQYIWAOc2dKoCrZ641amTOhk9vSJRYdiZtY28twl9ThwnKR9AEXEltaHVayfb9zGwftPQVLRoZiZtY08d0lNBN4N9ABdlS/RiPhMSyMr0OqNW1lw4Lh/mN3MbI/k6cP4d+A0YIhseJDKa1wqlYM1v9jOITOnFB2KmVlbydOHMTciTm55JG3i2d4+Bkpleg7YK0ZAMTPLLU8N4w5Jr255JG3i5xu3AXDw/q5hmJlVy1PDOAE4S9KTQD/ZnVIREUe1NLKCPNu7HYA5+04qOBIzs/aSJ2GcMtKTSzoZuBToBK6KiEuGbd8PWAIcBvQBH4iIFWnbR4EPko2S+yBwdkT0jTSWvJ7ZlF1izr4eQ8rMrFrTJqmIeCoingK2k3157xjqvBFJncDlZAlnIXCGpIXDdvsksDzVVs4kSy5IOgj4ELAoIo4kSzin5/1QL8azm/uYMaWbyRM6R+NyZmZjRp4Z906V9DPgSeBWYDXwgxznXgysiogn0nAi15LdbVVtIfAjgIh4FOiRNDtt6wImS+oCpgDP5Ljmi7Zucx8vm+7mKDOz4fJ0ev81cBzwWETMB04CfpLjuIPIngqvWJvKqt0PvAtA0mLgELK7sp4GPgf8HFgHbI6Im3Nc80V7tne7+y/MzGrIkzAGI2Ij0CGpIyJuAY7OcVytx6SHN2VdAuwnaTlwIXAfMJT6Nk4D5gMvB6ZKem/Ni0jnSlomadmGDRtyhNXYuk19vMz9F2Zmu8nT6b0pDQtyG9mYUuvJHuJrZi0wr2p9LsOalSKiFzgbQNkj5E+m11uBJyNiQ9r2HeD1wNeHXyQirgSuBFi0aFHTvpVG+odKbNw64BqGmVkNeWoYp5F1eH8U+E/gceCdOY67G1ggab6kCWSd1jdU7yBpRtoG2R1Rt6Uk8nOy8aumpERyEvBIng/0YmzaNgjAAfuM69HbzcxGJM/gg9XDgHw174kjYkjSBcBNZHc5LYmIhySdl7ZfARwBXC2pBDwMnJO2/VTSt4B7yWoz95FqEa3Uuz1LGNMndbf6UmZmY07dhCHp9og4QdIWdu17qDy413R2oYi4EbhxWNkVVct3AgvqHPsp4FPNrvFS6u3LEsa+k50wzMyGq5swIuKE9L7XDNu6uVLDcMIwM9tNwz4MSR2SVoxWMEXr3Z715U+flOdeADOzvUvDhBERZeB+SQePUjyFqtQw3CRlZra7PD+l5wAPSbqLqnkwIuLU+oeMTZVO72nu9DYz202ehPHplkfRJjZvH2RydycTuvLcbWxmtnfJc1vtraMRSDvo7Rt0c5SZWR15Bh88TtLdkl6QNCCpJKl3NIIbbb3bh5g+2R3eZma15Gl7uQw4A/gZMJnsiezLWhlUUTZvH/RDe2ZmdeRqrI+IVUBnRJQi4svAiS2NqiBukjIzqy9P+8u2NN7TckmfJRtufGprwypGb98gCw7cp+gwzMzaUp4axvvSfheQ3VY7D3h3K4Mqyrb+ElMmug/DzKyWPN+OxwI3plFkx/Uttn2DJSZ3e2pWM7Na8tQwTgUek/Q1SW9PU6aOOxFB31CZSd1+BsPMrJam344RcTbwCuA64D3A45KuanVgo22wFJTKwaQu1zDMzGrJVVuIiEFJPyAb5nwy2aRKH2xlYKOtb6gEwOQJThhmZrXkeXDvZElfAVYBvw1cRTa+1LjSN5gljInuwzAzqylPDeMs4FrgDyOiv7XhFKdvoAzgTm8zszryjCV1+mgEUrRKk5Q7vc3MavO3Y1JpknKnt5lZbU4YyfYBd3qbmTXihJH0DWV9GG6SMjOrrW4fhqQHyW6jrSkijmpJRAXpT01SEzpdwzAzq6VRp/c70vv56f1r6f33gW0ti6ggQ+UsN3q2PTOz2uomjIh4CkDS8RFxfNWmiyT9BPhMq4MbTYOlrEmqq1MFR2Jm1p7y/JyeKumEyoqk1zMOhzcfSH0YEzpdwzAzqyXPg3vnAEsk7UvWp7EZ+EBLoypApUmq2wnDzKymPA/u3QO8RtJ0QBGxufVhjT43SZmZNZZnLKnZkr4E/N+I2CxpoaRzRiG2UVVpknINw8ystjzfjl8BbgJentYfAz7SqoCKsuMuKScMM7Oa8nw7zoyIbwJlgIgYAkotjaoAg0NukjIzayRPwtgq6QDSQ3ySjiPr+B5XdvRhdDhhmJnVkucuqY8BNwCHpecvZpHNizGuDJaDCZ0dSE4YZma15LlL6l5JbwQOBwSsjIjBlkc2ygaHym6OMjNrINcUrcBioCftf6wkIuLqlkVVgKFy+A4pM7MGmiYMSV8DDgOWs7OzO4BxlTAGSmW6XcMwM6srTw1jEbAwIuqOXFuPpJOBS4FO4KqIuGTY9v2AJWQJqQ/4QESsSNtmkM0ffiRZgvpARNy5pzHkNThUdg3DzKyBPN+QK4CX7emJJXUClwOnAAuBMyQtHLbbJ4Hlaaj0M8mSS8WlwH9GxK8ArwEe2dMY9oSbpMzMGstTw5gJPCzpLqC/UhgRpzY5bjGwKiKeAJB0LXAa8HDVPguB/5XO96ikHkmzge3AG4Cz0rYBYCDPBxqpgZI7vc3MGsmTMP5qhOc+CFhTtb4WeO2wfe4H3gXcLmkxcAgwl6yvZAPwZUmvAe4BPhwRW0cYS1ODQ2U/5W1m1kCe22pvHeG5a/1cH94PcglwqaTlwIPAfcAQ0A0cC1wYET+VdClwEfA/d7uIdC5wLsDBBx88wlDdJGVm1kzdb0hJt6f3LZJ6q15bJPXmOPdaYF7V+lzgmeodIqI3Is6OiKPJ+jBmAU+mY9dGxE/Trt8iSyC7iYgrI2JRRCyaNWtWjrBqG3STlJlZQ41m3DshvU8b4bnvBhZImg88DZwOvKd6h3Qn1LbUR/FB4LaI6AV6Ja2RdHhErAROYte+j5fcgO+SMjNrKO+De0g6EJhUWY+InzfaPyKGJF1ANtJtJ7AkIh6SdF7afgVwBHC1pBJZQqgeNv1C4BuSJgBPAGfnjXUkhsrB5O7OVl7CzGxMy/Pg3qnA58mGN19P1jH9CPCqZsdGxI3AjcPKrqhavhNYUOfY5WTPgIyKwVKZaZNy508zs71OnjaYvwaOAx6LiPlkzUM/aWlUBRgsudPbzKyRPN+QgxGxEeiQ1BERtwBHtziuUTdY8m21ZmaN5GmD2SRpH+A2sj6F9WS3vo4rvkvKzKyxPD+pTyN78vqjwH8CjwPvbGVQRRhyk5SZWUN5Htyrfrr6qy2MpVAerdbMrLG6CUPSFnZ9MltpXUBExPQWxzaqBkt+DsPMrJFGD+6N9IG9MclNUmZmjeV68EDSscAJZDWM2yPivpZGVQCPVmtm1ljTn9SS/pKs7+IAsqHOvyLpL1od2GjzbbVmZo3lqWGcARwTEX0Aki4B7gX+ppWBjaZSOYjATVJmZg3k+YZcTdUYUsBEsltrx43BUhnATVJmZg3kqWH0Aw9J+iFZH8ZbyCY8+meAiPhQC+MbFZWE4SYpM7P68iSM76ZXxdLWhFKcwVJ297CbpMzM6suTMH4QEeurC6rmqRgX3CRlZtZcnp/UP5b0u5UVSR9n1xrHmFdJGN0drmGYmdWTp4ZxInClpN8BZpPNhbG4lUGNtlI5a5JyDcPMrL6mP6kjYh3ZoIOvA3qAqyPihRbHNapSvqBDThhmZvXkmXHvh8A64EhgLrBE0m0R8YlWBzdaKjUM5wszs/ryNNpfHhFnRsSmiFgBvB7Y3OK4RlVEljA6O5wxzMzqydMkdb2kQyS9ORV1A19obVijy01SZmbN5RlL6g+AbwH/mormAte3MqjRVmmScgXDzKy+PE1S5wPHA70AEfEz4MBWBjXaylFJGM4YZmb15EkY/RExUFmR1MWuEyuNeU4YZmbN5UkYt0r6JDBZ0luA64DvtTas0VXpw3Cnt5lZfXkSxkXABuBB4A+BG4FxNR9GpYbhCoaZWX1Nn8OIiDLwxfQal8plN0mZmTXjwZNwk5SZWR5OGLhJyswsj9wJQ9LUVgZSJDdJmZk1l+fBvddLephslFokvUbS/2l5ZKPITVJmZs3lqWH8E/BWYCNARNwPvKGVQY22nc9hFByImVkby9UkFRFrhhWVWhBLYUo7+jCcMczM6skzgdIaSa8HQtIE4EOk5qnxYsdotU4YZmZ15alhnEc2ntRBwFrg6LQ+bqQZWt3pbWbWQJ4ahiLi91seSYF8W62ZWXN5ahh3SLpZ0jmSZuzJySWdLGmlpFWSLqqxfT9J35X0gKS7JB05bHunpPskfX9PrrunPIGSmVlzeSZQWkA2dtSrgHslfV/Se5sdJ6kTuBw4BVgInCFp4bDdPgksj4ijgDOBS4dt/zCj0F/iJikzs+by3iV1V0R8DFgM/AL4ao7DFgOrIuKJNDz6tcBpw/ZZCPwoXeNRoEfSbABJc4G3A1flifHFKO+oYbT6SmZmY1eeB/emS3q/pB8AdwDryJJBMwcB1bfjrk1l1e4H3pWusxg4hGxGP8imgf1ToNwkvnMlLZO0bMOGDTnC2l3Zt9WamTWV5zf1/WR3Rn0mIl4ZEX8WEffkOK7Wt+/wiZcuAfaTtBy4ELgPGJL0DmB9nutExJURsSgiFs2aNStHWLvzBEpmZs3luUvq0Kj0Cu+ZtcC8qvW5wDPVO0REL3A2gLKf90+m1+nAqZLeBkwCpkv6ekQ07TsZiXKqw/g5DDOz+uomDElfiIiPADdI2i1hRMSpTc59N7BA0nzgabIk8J5h15gBbEt9HB8EbktJ5OL0QtKJwCdalSzAt9WameXRqIbxtfT+uZGcOCKGJF0A3AR0Aksi4iFJ56XtVwBHAFdLKgEPA+eM5Fov1o4mKd9Wa2ZWV92EUdV/cHRE7HK7q6QPA7c2O3lE3Eg2pWt12RVVy3cCC5qcYymwtNm1Xowdo9W6imFmVleeTu/31yg76yWOo1AerdbMrLlGfRhnkPU5zJd0Q9WmaaShzseLygRKvq3WzKy+Rn0YlWcuZgKfryrfAjzQyqBGmydQMjNrrlEfxlPAU8DrRi+cYpTKbpIyM2smz5Pex0m6W9ILkgYklST1jkZwo8VPepuZNZen0/sy4AzgZ8Bksucl/ncrgxpt4SYpM7Om8jzpTUSsktQZESXgy5LuaHFco6rku6TMzJrKkzC2palZl0v6LFlH+NTWhjW6PJaUmVlzeZqk3kf2pPYFwFay8aHe3cqgRlulScoJw8ysvqY1jHS3FMB24NOtDacYvkvKzKy5Rg/uPcjuw5HvkGbJGxfKnqLVzKypRjWMd4xaFAWrPLjn22rNzOpr9uDeXqFcDjdHmZk10bQPQ9IWdjZNTQC6ga0RMb2VgY2mcoSbo8zMmsjT6T2tel3Sb5JvTu8xoxxujjIzaybPbbW7iIjrgTe1IJbClMNNUmZmzeRpknpX1WoHsIgGd0+NReVyePIkM7Mm8jzp/c6q5SFgNXBaS6IpSCnCD+2ZmTWRpw/j7NEIpEgR4HxhZtZYniap+cCFQE/1/hFxauvCGl2+S8rMrLk8TVLXA18CvgeUWxtOMUplN0mZmTWTJ2H0RcQ/tzySApUDOlzDMDNrKE/CuFTSp4Cbgf5KYUTc27KoRln4tlozs6byJIxXkw1x/iZ2NkkF4+hZDDdJmZk1lydh/BZwaEQMtDqYopTDc2GYmTWT50nv+4EZrQ6kSBFBxx4/825mtnfJU8OYDTwq6W527cMYN7fV+sE9M7Pm8iSMT7U8ioKVAw8NYmbWRJ4nvW8djUCKVI7wk95mZk14PgwqEyg5Y5iZNeL5MPDQIGZmeXg+DKBU9gRKZmbNeD4MsttqO31brZlZQ54Pg8qMe65hmJk10tL5MCSdDFwKdAJXRcQlw7bvBywBDgP6gA9ExApJ84CrgZeRDUdyZURcOtI4mil5Tm8zs6aaNsRI+qqkGVXr+0lakuO4TuBy4BRgIXCGpIXDdvsksDwijgLOJEsukNVkPh4RRwDHAefXOPYlExF0Ol+YmTWUp+X+qIjYVFmJiF8Cx+Q4bjGwKiKeSONQXcvuTVkLgR+l8z4K9EiaHRHrKqPhRsQW4BHgoBzXHBE3SZmZNZcnYXSkpiMAJO1Pvr6Pg4A1Vetr2f1L/37gXem8i4FDgLnVO0jqIUtQP81xzRHxaLVmZs3l+eL/PHCHpG+R3R31u8Df5jiu1jfw8LurLiGbb2M58CBwH1lzVHYCaR/g28BHIqK35kWkc4FzAQ4++OAcYe0um0BpRIeame018nR6Xy1pGdmzFwLeFREP5zj3WmBe1fpc4Jlh5+4FzgZQ1uv8ZHohqZssWXwjIr7TIL4rgSsBFi1aNKLbfbPRap0xzMwayVPDICWIPEmi2t3AAknzgaeB04H3VO+QOtO3pT6ODwK3RURvSh5fAh6JiH/cw+vusVI5mNjlJikzs0ZyJYyRiIghSRcAN5HdVrskIh6SdF7afgVwBHC1pBJZQjonHX482Sx/D6bmKoBPRsSNrYjVc3qbmTXXsoQBkL7gbxxWdkXV8p3AghrH3U7tPpCW8JzeZmbNueEeT6BkZpaHEwZQLntObzOzZpwwqDy4V3QUZmbtzQkDz4dhZpaHEwbpLik3SZmZNeSEQTZFq/OFmVljThi4ScrMLA8nDNwkZWaWhxMG2dAgzhdmZo05YVCZQMkZw8ysEScM3CRlZpaHEwZpaBD/nzAza8hfk1QGH3QNw8ysEScMPEWrmVkeThhU+jCKjsLMrL05YZAGH3TGMDNryAmDbGgQN0mZmTXmhEHWJOWhQczMGnPCAN76qtkcMWda0WGYmbW1ls7pPVZ84fRjig7BzKztuYZhZma5OGGYmVkuThhmZpaLE4aZmeXihGFmZrk4YZiZWS5OGGZmlosThpmZ5aKIKDqGl4ykDcBTIzx8JvD8SxhOK4yFGGFsxDkWYoSxEedYiBHGRpxFxHhIRMzKs+O4ShgvhqRlEbGo6DgaGQsxwtiIcyzECGMjzrEQI4yNONs9RjdJmZlZLk4YZmaWixPGTlcWHUAOYyFGGBtxjoUYYWzEORZihLERZ1vH6D4MMzPLxTUMMzPLZa9PGJJOlrRS0ipJFxUcyxJJ6yWtqCrbX9IPJf0sve9Xte3iFPdKSW8dpRjnSbpF0iOSHpL04XaLU9IkSXdJuj/F+Ol2i3FYvJ2S7pP0/XaMU9JqSQ9KWi5pWTvGmK47Q9K3JD2a/n6+rp3ilHR4+n9YefVK+kg7xdhUROy1L6ATeBw4FJgA3A8sLDCeNwDHAiuqyj4LXJSWLwL+Pi0vTPFOBOanz9E5CjHOAY5Ny9OAx1IsbRMnIGCftNwN/BQ4rp1iHBbvx4B/A77fpn/mq4GZw8raKsZ07a8CH0zLE4AZ7Rhnun4n8CxwSLvGWDPuIi9e9At4HXBT1frFwMUFx9TDrgljJTAnLc8BVtaKFbgJeF0B8f478JZ2jROYAtwLvLYdYwTmAj8C3lSVMNoqzjoJo91inA48SeqXbdc4q673G8BP2jnGWq+9vUnqIGBN1fraVNZOZkfEOoD0fmAqLzx2ST3AMWS/4NsqztTMsxxYD/wwItouxuQLwJ8C5aqydoszgJsl3SPp3DaN8VBgA/Dl1Lx3laSpbRhnxenANWm5XWPczd6eMFSjbKzcNlZo7JL2Ab4NfCQiehvtWqOs5XFGRCkijib7Bb9Y0pENdi8kRknvANZHxD15D6lRNhp/5sdHxLHAKcD5kt7QYN+iYuwia879l4g4BthK1rxTT2H/fiRNAE4Frmu2a42yQr+f9vaEsRaYV7U+F3imoFjqeU7SHID0vj6VFxa7pG6yZPGNiPhOu8YJEBGbgKXAyW0Y4/HAqZJWA9cCb5L09XaLMyKeSe/rge8Ci9stxnTdtakmCfAtsgTSbnFClnjvjYjn0no7xljT3p4w7gYWSJqfsv7pwA0FxzTcDcD70/L7yfoMKuWnS5ooaT6wALir1cFIEvAl4JGI+Md2jFPSLEkz0vJk4M3Ao+0UI0BEXBwRcyOih+zv3n9FxHvbKU5JUyVNqyyTtb2vaKcYASLiWWCNpMNT0UnAw+0WZ3IGO5ujKrG0W4y1FdmB0g4v4G1kd/o8Dvx5wbFcA6wDBsl+XZwDHEDWKfqz9L5/1f5/nuJeCZwySjGeQFYtfgBYnl5va6c4gaOA+1KMK4C/TOVtE2ONmE9kZ6d328RJ1jdwf3o9VPk30k4xVl33aGBZ+nO/Htiv3eIkuwljI7BvVVlbxdjo5Se9zcwsl72Eq+hlAAAEgUlEQVS9ScrMzHJywjAzs1ycMMzMLBcnDDMzy8UJw8zMcnHCsHFL0lJJLZ8fWdKH0uio32j1tYqURoP946LjsOI4YZjVIKlrD3b/Y+BtEfH7rYqnTcwg+6y2l3LCsEJJ6km/zr+obO6Km9PT2bvUECTNTENoIOksSddL+p6kJyVdIOljadC5/5a0f9Ul3ivpDkkrJC1Ox09VNvfI3emY06rOe52k7wE314j1Y+k8KyR9JJVdQfZw2w2SPjps/05Jn1M2l8QDki5M5Sel6z6Y4piYyldL+jtJd0paJulYSTdJelzSeWmfEyXdJum7kh6WdIWkjrTtjHTOFZL+viqOFyT9rbL5Qf5b0uxUPkvSt9P/h7slHZ/K/yrFtVTSE5I+lE51CXCYsrkc/kHSnBTL8nTNXxvxXwQbG4p+ctCvvftFNpz7EHB0Wv8m8N60vBRYlJZnAqvT8lnAKrL5OGYBm4Hz0rZ/IhsQsXL8F9PyG0jDxgN/V3WNGWRP+k9N511L1ZO2VXH+KvBg2m8fsqeej0nbVjNs+O9U/kdkY251pfX9gUlkI5C+MpVdXRXvauCPqj7HA1WfcX0qPxHoI0tSncAPgd8GXg78PO3bBfwX8JvpmADemZY/C/xFWv434IS0fDDZcC8AfwXcQTYPw0yyJ5O72X3o/Y+z88nvTmBa0X+f/Grta0+q3Wat8mRELE/L95B9MTVzS0RsAbZI2gx8L5U/SDY0SMU1ABFxm6TpaYyp3yAb9O8TaZ9JZF+YkA2F/osa1zsB+G5EbAWQ9B3g18iGIKnnzcAVETGUYviFpNekz/tY2uerwPlkw5zDzrHMHiSbBKryGfsq42MBd0XEEymOa1Jsg8DSiNiQyr9BliSvBwaA76dj7yGbv6QS38JseDAAplfGjQL+IyL6gX5J64HZNT7f3cASZYNRXl/1Z2jjlBOGtYP+quUSMDktD7Gz2XRSg2PKVetldv17PXzsmyAbNvrdEbGyeoOk15INi11LraGmm1GN6zc7T/XnGP4ZK5+r3meqZzAiKseUqs7TQTYhz/ZdAswSyPA/k92+K1ISfgPwduBrkv4hIq5uEIeNce7DsHa2mqwpCLJml5H4PQBJJwCbI2Iz2cxlFyp9M0o6Jsd5bgN+U9KUNGrrbwE/bnLMzcB5lQ701LfyKNAj6RVpn/cBt+7hZ1qsbITlDrLPdzvZJFZvTH09nWQjojY7783ABZUVSUc32X8LWRNZZf9DyJrKvkg2gvGxe/g5bIxxDcPa2eeAb0p6H1mb/Ej8UtIdZFN4fiCV/TVZE9ADKWmsBt7R6CQRca+kr7BzeOmrIqJRcxTAVcAr03UGyfpTLpN0NnBdSiR3A1fs4We6k6wD+tVkiey7EVGWdDFwC1lt48aI+PcG5wD4EHC5pAfIvgtuA86rt3NEbJT0E0krgB+QjQT8J+mzvQCcuYefw8YYj1ZrNoZIOhH4REQ0THBmreAmKTMzy8U1DDMzy8U1DDMzy8UJw8zMcnHCMDOzXJwwzMwsFycMMzPLxQnDzMxy+f/6Zz0vuORF8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pca = PCA().fit(stack)\n",
    "plt.plot(np.cumsum(pca.explained_variance_ratio_))\n",
    "plt.xlabel('number of components')\n",
    "plt.ylabel('cumulative explained variance');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decompose original data using 20 PCs (~98.6% of the initial variance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(735, 20)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca = PCA(n_components=20).fit(stack)\n",
    "stackT = pca.transform(stack)\n",
    "\n",
    "stackT.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regenerate initial data using the decomposed components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(735, 160000)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec_stack = pca.inverse_transform(stackT)\n",
    "\n",
    "rec_stack.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write video comparing RGB time series of the original data (left) and recreated PCA (right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Lossy conversion from float64 to uint8. Range [0, 1]. Convert image to uint8 prior to saving to suppress this warning.\n"
     ]
    }
   ],
   "source": [
    "import imageio\n",
    "\n",
    "writer = imageio.get_writer('video.mp4', fps=2)\n",
    "\n",
    "for t in range(105):\n",
    "    orig_r = stack[0*105+t].reshape((400,400))\n",
    "    rec_r = rec_stack[0*105+t].reshape((400,400))\n",
    "    orig_g = stack[1*105+t].reshape((400,400))\n",
    "    rec_g = rec_stack[1*105+t].reshape((400,400))\n",
    "    orig_b = stack[2*105+t].reshape((400,400))\n",
    "    rec_b = rec_stack[2*105+t].reshape((400,400))\n",
    "    \n",
    "    writer.append_data(3*np.concatenate((np.dstack((orig_r,orig_g,orig_b)),np.dstack((rec_r,rec_g,rec_b))), axis=1))\n",
    "\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video src=\"video.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Video\n",
    "\n",
    "Video(\"video.mp4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefs = pca.transform(stack)\n",
    "comps = pca.components_\n",
    "\n",
    "def plot_img(band, t):     \n",
    "    f, axarr = plt.subplots(1,2)\n",
    "    orig = stack[band*105+t].reshape((400,400))\n",
    "    rec = rec_stack[band*105+t].reshape((400,400))\n",
    "    \n",
    "    axarr[0].imshow(orig, vmin=0, vmax=orig.max())\n",
    "    axarr[1].imshow(rec, vmin=0, vmax=orig.max())\n",
    "    \n",
    "    print(orig.max(), rec.max())\n",
    "    \n",
    "def plot_rgb(t):     \n",
    "    f, axarr = plt.subplots(1,2)\n",
    "    \n",
    "    orig_r = stack[0*105+t].reshape((400,400))\n",
    "    rec_r = rec_stack[0*105+t].reshape((400,400))\n",
    "    orig_g = stack[1*105+t].reshape((400,400))\n",
    "    rec_g = rec_stack[1*105+t].reshape((400,400))\n",
    "    orig_b = stack[2*105+t].reshape((400,400))\n",
    "    rec_b = rec_stack[2*105+t].reshape((400,400))\n",
    "    \n",
    "    axarr[0].imshow(np.dstack((orig_r,orig_g,orig_b))*3)\n",
    "    axarr[1].imshow(np.dstack((rec_r,rec_g,rec_b))*3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}