lstm-cnn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNN-LSTM model for in-season crop type classification\n",
    "\n",
    "This notebook implements a dual-branch recurrent and convolutional neural network for crop type classification. Output classes are corn, soybean, or other. The input to the model includes optical (Harmonized Landsat and Sentinel-2) and microwave (Sentinel-1) observations at 30m spatial resolution on days of the year (DOYs) selected based on growth stages detected in each pixel (greenup, peak growth, and senescence). Models are trained and tested at these three stages of the growing season as well. The labels for training and evaluation are from the USDA NASS Cropland Data Layer (CDL), which provides crop-specific land cover classes at 30m resolution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import necessary libraries and packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard libraries\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import os\n",
    "from time import time\n",
    "\n",
    "# ML libraries\n",
    "import tensorflow as tf\n",
    "import keras\n",
    "from keras.layers.core import *\n",
    "from keras.models import Sequential, Model, load_model\n",
    "from keras.layers import Dense\n",
    "from keras.layers import concatenate\n",
    "from keras.layers import Input\n",
    "from keras.layers import LSTM\n",
    "from keras.layers import Dropout\n",
    "from keras.layers import Conv1D, Conv2D\n",
    "from keras.layers import MaxPool2D\n",
    "from keras.layers import Dropout\n",
    "from keras.utils import to_categorical\n",
    "from keras.optimizers import Adam\n",
    "from sklearn import metrics\n",
    "\n",
    "# Define some constants\n",
    "CORN = 0\n",
    "SOYBEAN = 1\n",
    "OTHER = 2\n",
    "N_CLASSES = 3\n",
    "N_BANDS = 6\n",
    "N_TIMESTEPS = 3\n",
    "BATCH_SIZE = 4096\n",
    "N_EPOCHS = 25\n",
    "HEIGHT = 3660\n",
    "WIDTH = 3660\n",
    "\n",
    "# Fix random seed for reproducibility\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the dataset\n",
    "The data are stored as numpy arrays with dimension height x width x bands x timesteps. All of the reflectance values are in the range [0,1]. We also add two spectral indices (NDWI and LSWI) and SAR bands (VV, VH, incidence angle/IA)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_lswi_channel(X):\n",
    "    _X = np.ndarray([HEIGHT, WIDTH, X.shape[2]+1, N_TIMESTEPS])\n",
    "    # copy the values from the original array\n",
    "    for i in range(X.shape[2]):\n",
    "        _X[:,:,i,:] = X[:,:,i,:]\n",
    "    # calculate values for LSWI channel\n",
    "    for i in range(N_TIMESTEPS):\n",
    "        lswi = (X[:,:,NIR,i]-X[:,:,SWIR1,i])/(X[:,:,NIR,i]+X[:,:,SWIR1,i])\n",
    "        _X[:,:,-1,i] = lswi\n",
    "    # make sure we didn't introduce any NaNs\n",
    "    _X[np.where(np.isnan(_X))] = 0\n",
    "    return _X\n",
    "\n",
    "def add_ndwi_channel(X):\n",
    "    _X = np.ndarray([HEIGHT, WIDTH, X.shape[2]+1, N_TIMESTEPS])\n",
    "    # copy the values from the original array\n",
    "    for i in range(X.shape[2]):\n",
    "        _X[:,:,i,:] = X[:,:,i,:]\n",
    "    # calculate values for NDWI channel\n",
    "    for i in range(N_TIMESTEPS):\n",
    "        ndwi = (X[:,:,GREEN,i]-X[:,:,SWIR1,i])/(X[:,:,GREEN,i]+X[:,:,SWIR1,i])\n",
    "        _X[:,:,-1,i] = ndwi\n",
    "    # make sure we didn't introduce any NaNs\n",
    "    _X[np.where(np.isnan(_X))] = 0\n",
    "    return _X\n",
    "\n",
    "def add_sar_channel(X, band, path):\n",
    "    _X = np.ndarray([HEIGHT, WIDTH, X.shape[2]+1, N_TIMESTEPS])\n",
    "    # copy the values from the original array\n",
    "    for i in range(X.shape[2]):\n",
    "        _X[:,:,i,:] = X[:,:,i,:]\n",
    "    # load the corresponding SAR band\n",
    "    if band=='vv' or band=='VV':\n",
    "        sarpath = path.replace('pheno_timeseries', 'vv_timeseries')\n",
    "        #sarpath = path.replace('fixed_timeseries', 'vv_timeseries')\n",
    "    elif band=='vh' or band=='VH':\n",
    "        sarpath = path.replace('pheno_timeseries', 'vh_timeseries')\n",
    "        #sarpath = path.replace('fixed_timeseries', 'vh_timeseries')\n",
    "    elif band=='ia' or band=='IA':\n",
    "        sarpath = path.replace('pheno_timeseries', 'ia_timeseries')\n",
    "        #sarpath = path.replace('fixed_timeseries', 'ia_timeseries')\n",
    "    sar = np.load(sarpath).astype(np.float32)\n",
    "    for i in range(N_TIMESTEPS):\n",
    "        _X[:,:,-1,i] = sar[...,i]\n",
    "    # make sure we didn't introduce any NaNs\n",
    "    _X[np.where(np.isnan(_X))] = 0\n",
    "    return _X\n",
    "\n",
    "def load_data(x_path, y_path, flatten=True, convert_nans=True):\n",
    "    # Load the time series image data\n",
    "    X = np.load(x_path).astype(np.float32)\n",
    "    # Load the associated labels\n",
    "    Y = np.load(y_path).astype(np.int8)\n",
    "    \n",
    "    # Convert all the NaNs to zeros\n",
    "    if convert_nans:\n",
    "        X[np.where(np.isnan(X))] = 0\n",
    "        \n",
    "    X[np.where(X==0)] = 0.00000001\n",
    "    # Add band indices\n",
    "    X = add_lswi_channel(X)\n",
    "    X = add_ndwi_channel(X)\n",
    "    X = add_sar_channel(X, 'vv', x_path)\n",
    "    X = add_sar_channel(X, 'vh', x_path)\n",
    "    X = add_sar_channel(X, 'ia', x_path)\n",
    "    if flatten:\n",
    "        # Reduce the h x w x b x t dataset to h*w x b x t\n",
    "        X = np.reshape(X, (X.shape[0]*X.shape[1], X.shape[2], X.shape[3]))\n",
    "        Y = np.reshape(Y, (Y.shape[0]*Y.shape[1]))\n",
    "    assert X.shape[0] == Y.shape[0] \n",
    "    return X, Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training data used was collected over MGRS tile 16TBL in northern Illinois in 2018 and 2017. For validation, we will use the pixels in the southwest quadrant of the 2017-2018 16TBL data. We load 2017 and 2018 separately first and then split them up geographically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = load_data(x_path='data/inputs_3doy/pheno_timeseries_16TBL_2018.npy', y_path='data/cdl_labels_16TBL_2018.npy', flatten=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_val, y_val = load_data(x_path='data/inputs_3doy/pheno_timeseries_16TBL_2017.npy', y_path='data/cdl_labels_16TBL_2017.npy', flatten=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dois_2018 = np.load('data/phenology/16TBL_2018_phenology.npy').astype(np.int16)\n",
    "dois_2017 = np.load('data/phenology/16TBL_2017_phenology.npy').astype(np.int16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable delta_s is the number of days elapsed between the senescence and greenup DOYs in each pixel. Similarly, delta_p is the number of days elapsed between the peak and greenup DOYs in each pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_2018_s = dois_2018[...,2]-dois_2018[...,0]\n",
    "delta_2017_s = dois_2017[...,2]-dois_2017[...,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_2018_p = dois_2018[...,1]-dois_2018[...,0]\n",
    "delta_2017_p = dois_2017[...,1]-dois_2017[...,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The homogeneity mask is used to select training pixels only within regions of the CDL that are spatially consistent and thus more likely to be correctly classified in the CDL reference labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_2018 = np.load('data/cdl_labels_16TBL_2018_homogmask_3x3.npy')\n",
    "mask_2017 = np.load('data/cdl_labels_16TBL_2017_homogmask_3x3.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the eastern half and northwest quadrant of the 2017 and 2018 data out for training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use eastern half and NW quadrant of 2017 and 2018 for training\n",
    "east_2018 = X_train[:,int(WIDTH/2):]\n",
    "east_2018_y = y_train[:,int(WIDTH/2):]\n",
    "east_2018_mask = mask_2018[:,int(WIDTH/2):]\n",
    "east_2018_delta_p = delta_2018_p[:,int(WIDTH/2):]\n",
    "east_2018_delta_s = delta_2018_s[:,int(WIDTH/2):]\n",
    "\n",
    "east_2017 = X_val[:,int(WIDTH/2):]\n",
    "east_2017_y = y_val[:,int(WIDTH/2):]\n",
    "east_2017_mask = mask_2017[:,int(WIDTH/2):]\n",
    "east_2017_delta_s = delta_2017_s[:,int(WIDTH/2):]\n",
    "east_2017_delta_p = delta_2017_p[:,int(WIDTH/2):]\n",
    "\n",
    "nw_2018 = X_train[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2018_y = y_train[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2018_mask = mask_2018[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2018_delta_s = delta_2018_s[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2018_delta_p = delta_2018_p[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "\n",
    "nw_2017 = X_val[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2017_y = y_val[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2017_mask = mask_2017[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2017_delta_s = delta_2017_s[:int(HEIGHT/2),:int(WIDTH/2)]\n",
    "nw_2017_delta_p = delta_2017_p[:int(HEIGHT/2),:int(WIDTH/2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the southwestern quadrant of the 2017 and 2018 data out for validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use SW quadrant of 2017 and 2018 for validation\n",
    "sw_2018 = X_train[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2018_y = y_train[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2018_mask = mask_2018[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2018_delta_s = delta_2018_s[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2018_delta_p = delta_2018_p[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "\n",
    "sw_2017 = X_val[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2017_y = y_val[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2017_mask = mask_2017[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2017_delta_s = delta_2017_s[int(HEIGHT/2):,:int(WIDTH/2)]\n",
    "sw_2017_delta_p = delta_2017_p[int(HEIGHT/2):,:int(WIDTH/2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual input to the model is a pair containing the spectrum at 3 DOYs in one pixel and a 5x5 patch surrounding that pixel. This function creates these input pairs from the larger rasters. There is an option to select a random subset or to get the pairs for all pixels in the raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make pixel and patch pairs from the raster\n",
    "def make_inputs(raster, labels, deltas_s, deltas_p, subset=None, kernel_size=5, mask=None):\n",
    "    # Iterate through pixels with buffer for kernel size\n",
    "    w = int(kernel_size/2)\n",
    "    margin = w+kernel_size%2\n",
    "    if subset:\n",
    "        # Random rows and cols gives subset**2 examples\n",
    "        rand_i = np.random.randint(0, raster.shape[0]-kernel_size, size=subset)\n",
    "        rand_j = np.random.randint(0, raster.shape[1]-kernel_size, size=subset)\n",
    "        # Initialize arrays to hold the examples\n",
    "        pixels = []\n",
    "        patches = []\n",
    "        y = []\n",
    "        delta_s = []\n",
    "        delta_p = []\n",
    "        for idx_i, i in enumerate(rand_i):\n",
    "            for idx_j, j in enumerate(rand_j):\n",
    "                # only take examples from homogeneous regions of CDL\n",
    "                if np.any(mask[i-w:i+w+1, j-w:j+w+1] != 1):\n",
    "                    continue\n",
    "                # Store the pixel representation\n",
    "                pixel = raster[i,j]\n",
    "                # Store the patch representation\n",
    "                patch = raster[i-w:i+w+1, j-w:j+w+1]\n",
    "                if patch.shape[0] != 5 or patch.shape[1] != 5:\n",
    "                    continue\n",
    "                patch = np.reshape(patch, [patch.shape[0], patch.shape[1], patch.shape[2]*patch.shape[3]],order='F')\n",
    "                patches.append(patch)\n",
    "                pixels.append(pixel)\n",
    "                # Store the label\n",
    "                y.append(labels[i,j])\n",
    "                # Store the senescence-greenup deltas\n",
    "                delta_s.append(deltas_s[i,j])\n",
    "                delta_p.append(deltas_p[i,j])\n",
    "    else:\n",
    "        # Initialize arrays to hold the examples\n",
    "        pixels = np.ndarray([HEIGHT*WIDTH, N_BANDS, N_TIMESTEPS])\n",
    "        patches = np.ndarray([HEIGHT*WIDTH, kernel_size, kernel_size, N_BANDS*N_TIMESTEPS])\n",
    "        y = np.ndarray([HEIGHT*WIDTH]).astype(np.int8)\n",
    "        delta_s = np.ndarray([HEIGHT*WIDTH]).astype(np.int32)\n",
    "        delta_p = np.ndarray([HEIGHT*WIDTH]).astype(np.int32)\n",
    "        for i in range(margin, raster.shape[0]-margin):\n",
    "            for j in range(margin, raster.shape[1]-margin):\n",
    "                # Store the pixel representation\n",
    "                pixel = raster[i,j]\n",
    "                pixels[i*WIDTH+j] = pixel\n",
    "                # Store the patch representation\n",
    "                patch = raster[i-w:i+w+1, j-w:j+w+1]\n",
    "                patch = np.reshape(patch, [patch.shape[0], patch.shape[1], patch.shape[2]*patch.shape[3]], order='F')\n",
    "                patches[i*WIDTH+j] = patch\n",
    "                # Store the label\n",
    "                y[i*WIDTH+j] = labels[i,j]   \n",
    "                # Store the senescence-greenup deltas\n",
    "                delta_s[i*WIDTH+j] = deltas_s[i,j]\n",
    "                delta_p[i*WIDTH+j] = deltas_p[i,j]\n",
    "                \n",
    "    return np.array(pixels), np.array(patches), np.array(y), np.array(delta_s), np.array(delta_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the input pairs for each of the training regions. Note some need to be run multiple times if a pixel in the random subset falls on the edge of the raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "east_2018_px, east_2018_pa, east_2018_labels, east_2018_delta_s, east_2018_delta_p = make_inputs(east_2018, east_2018_y, east_2018_delta_s, east_2018_delta_p, subset=400, mask=east_2018_mask) # 2 times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "east_2017_px, east_2017_pa, east_2017_labels, east_2017_delta_s, east_2017_delta_p = make_inputs(east_2017, east_2017_y, east_2017_delta_s, east_2017_delta_p, subset=400, mask=east_2017_mask) # 4 times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nw_2018_px, nw_2018_pa, nw_2018_labels, nw_2018_delta_s, nw_2018_delta_p = make_inputs(nw_2018, nw_2018_y, nw_2018_delta_s, nw_2018_delta_p, subset=100, mask=nw_2018_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nw_2017_px, nw_2017_pa, nw_2017_labels, nw_2017_delta_s, nw_2017_delta_p = make_inputs(nw_2017, nw_2017_y, nw_2017_delta_s, nw_2017_delta_p, subset=100, mask=nw_2017_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concatenate the eastern and northwestern regions to form the training examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_px = np.concatenate([east_2018_px, east_2017_px, nw_2018_px, nw_2017_px], axis=0)\n",
    "X_train_pa = np.concatenate([east_2018_pa, east_2017_pa, nw_2018_pa, nw_2017_pa], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = np.concatenate([east_2018_labels, east_2017_labels, nw_2018_labels, nw_2017_labels], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_train_s = np.concatenate([east_2018_delta_s, east_2017_delta_s, nw_2018_delta_s, nw_2017_delta_s], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_train_p = np.concatenate([east_2018_delta_s, east_2017_delta_s, nw_2018_delta_s, nw_2017_delta_p], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsample the dataset to equalize the number of examples we have for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subsample(X_px, X_pa, y, delta_s, delta_p):\n",
    "    # Get the minimum number of samples in a class\n",
    "    n_corn = len(np.where(y==CORN)[0])\n",
    "    n_soy = len(np.where(y==SOYBEAN)[0])\n",
    "    n_other = len(np.where(y==OTHER)[0])\n",
    "    lim = np.min([n_corn, n_soy, n_other])\n",
    "    # Sub-sample the classes to the same number\n",
    "    X_px = np.concatenate([X_px[np.where(y==CORN)][:lim], X_px[np.where(y==SOYBEAN)][:lim], X_px[np.where(y==OTHER)][:lim]])\n",
    "    X_pa = np.concatenate([X_pa[np.where(y==CORN)][:lim], X_pa[np.where(y==SOYBEAN)][:lim], X_pa[np.where(y==OTHER)][:lim]])\n",
    "    delta_s = np.concatenate([delta_s[np.where(y==CORN)][:lim], delta_s[np.where(y==SOYBEAN)][:lim], delta_s[np.where(y==OTHER)][:lim]])\n",
    "    delta_p = np.concatenate([delta_p[np.where(y==CORN)][:lim], delta_p[np.where(y==SOYBEAN)][:lim], delta_p[np.where(y==OTHER)][:lim]])\n",
    "    y = np.concatenate([y[np.where(y==CORN)][:lim], y[np.where(y==SOYBEAN)][:lim], y[np.where(y==OTHER)][:lim]])\n",
    "    return X_px, X_pa, y, delta_s, delta_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_px, X_train_pa, y_train, delta_train_s, delta_train_p = subsample(X_train_px, X_train_pa, y_train, delta_train_s, delta_train_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the input pairs for each of the validation regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_2018_px, sw_2018_pa, sw_2018_labels, sw_2018_delta_s, sw_2018_delta_p = make_inputs(sw_2018, sw_2018_y, sw_2018_delta_s, sw_2018_delta_p, subset=100, mask=sw_2018_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_2017_px, sw_2017_pa, sw_2017_labels, sw_2017_delta_s, sw_2017_delta_p = make_inputs(sw_2017, sw_2017_y, sw_2017_delta_s, sw_2017_delta_p, subset=100, mask=sw_2017_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concatenate the southwestern regions to form the validation examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_val_px = np.concatenate([sw_2018_px, sw_2017_px], axis=0)\n",
    "X_val_pa = np.concatenate([sw_2018_pa, sw_2017_pa], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_val = np.concatenate([sw_2018_labels, sw_2017_labels], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_val_s = np.concatenate([sw_2018_delta_s, sw_2017_delta_s], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_val_p = np.concatenate([sw_2018_delta_p, sw_2017_delta_p], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsample the dataset to equalize the number of examples we have for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_val_px, X_val_pa, y_val, delta_val_s, delta_val_p = subsample(X_val_px, X_val_pa, y_val, delta_val_s, delta_val_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Swap temporal and spectral axes\n",
    "Currently our data have dimension # pixels x # bands (features) x # time step. The LSTM expects them in the order # pixels x # time steps x # bands (features) so we need to re-order the last two dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_px = np.swapaxes(X_train_px, 1, 2)\n",
    "X_val_px = np.swapaxes(X_val_px, 1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make one-hot labels\n",
    "The labels for our three crop type classes---corn, soybean, and other---are currently integers. We need to make these one-hot vectors (e.g., corn: [1, 0, 0])."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_val = to_categorical(y_val, num_classes=N_CLASSES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = to_categorical(y_train, num_classes=N_CLASSES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Label smoothing has been shown to improve model generalization and learning speed: http://papers.nips.cc/paper/8717-when-does-label-smoothing-help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Smooth the labels\n",
    "def smooth_labels(labels, factor=0.1):\n",
    "    # smooth the labels\n",
    "    labels = labels * (1 - factor)\n",
    "    labels = labels + (factor / labels.shape[1])\n",
    "    return labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = smooth_labels(y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_val = smooth_labels(y_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardize the data\n",
    "Standardize the data by subtracting the mean and dividing by the standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Different scaler for each timestep\n",
    "gscaler = StandardScaler()\n",
    "pscaler = StandardScaler()\n",
    "sscaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_px[:,0] = gscaler.fit_transform(X_train_px[:,0])\n",
    "X_train_px[:,1] = pscaler.fit_transform(X_train_px[:,1])\n",
    "X_train_px[:,2] = sscaler.fit_transform(X_train_px[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_val_px[:,0] = gscaler.transform(X_val_px[:,0])\n",
    "X_val_px[:,1] = pscaler.transform(X_val_px[:,1])\n",
    "X_val_px[:,2] = sscaler.transform(X_val_px[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_pa = np.mean(X_train_pa, axis=0)\n",
    "std_pa = np.std(X_train_pa, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standardize(X, mu, std):\n",
    "    return (X - mu)/std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_pa = standardize(X_train_pa, mu_pa, std_pa)\n",
    "X_val_pa = standardize(X_val_pa, mu_pa, std_pa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize an example of the standardised data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 40\n",
    "y_train[t]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X_train_px[t,0],c='g',label='Greenup')\n",
    "plt.plot(X_train_px[t,1],c='y',label='Peak')\n",
    "plt.plot(X_train_px[t,2],c='r',label='Senescence')\n",
    "plt.axis('off')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=N_TIMESTEPS, ncols=int(X_train_pa.shape[-1]/N_TIMESTEPS))\n",
    "for b in range(int(X_train_pa.shape[-1]/N_TIMESTEPS)):\n",
    "    for d in range(N_TIMESTEPS):\n",
    "        axes[d,b].imshow(X_train_pa[t,:,:,d*int(X_train_pa.shape[-1]/N_TIMESTEPS) + b], cmap='gray')\n",
    "        axes[d,b].axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model\n",
    "We have separated out different components of the full model to evaluate the contribution of each component. Within each model or sub-model, there is a model for early-season, mid-season, and late-season classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PlotLearning(keras.callbacks.Callback):\n",
    "    def on_train_begin(self, logs={}):\n",
    "        self.i = 0\n",
    "        self.x = []\n",
    "        self.losses = []\n",
    "        self.val_losses = []\n",
    "        self.acc = []\n",
    "        self.val_acc = []\n",
    "        self.fig = plt.figure()\n",
    "        \n",
    "        self.logs = []\n",
    "\n",
    "    def on_epoch_end(self, epoch, logs={}):\n",
    "        \n",
    "        self.logs.append(logs)\n",
    "        self.x.append(self.i)\n",
    "        self.losses.append(logs.get('loss'))\n",
    "        self.val_losses.append(logs.get('val_loss'))\n",
    "        self.acc.append(logs.get('accuracy'))\n",
    "        self.val_acc.append(logs.get('val_accuracy'))\n",
    "        self.i += 1\n",
    "        f, (ax1, ax2) = plt.subplots(1, 2, sharex=True)\n",
    "        f.tight_layout()\n",
    "        clear_output(wait=True)\n",
    "        \n",
    "        ax1.set_yscale('log')\n",
    "        ax1.plot(self.x, self.losses, label=\"training loss\")\n",
    "        ax1.plot(self.x, self.val_losses, label=\"validation loss\")\n",
    "        ax1.legend()\n",
    "        \n",
    "        ax2.plot(self.x, self.acc, label=\"training accuracy\")\n",
    "        ax2.plot(self.x, self.val_acc, label=\"validation accuracy\")\n",
    "        ax2.legend()\n",
    "        \n",
    "        plt.show();\n",
    "        \n",
    "plot = PlotLearning()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Early-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## CNN branch only ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 1*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(conv1_2d)\n",
    "mp1 = MaxPool2D()(conv2_2d)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(flat_cnn)\n",
    "\n",
    "# output layer -> corn, soybeans, other crop, other non-crop\n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[cnn_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_pa[:,:,:,:1*N_BANDS]], y=y_train, validation_data=([X_val_pa[:,:,:,:1*N_BANDS]], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## CNN branch only ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 2*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(conv1_2d)\n",
    "mp1 = MaxPool2D()(conv2_2d)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(flat_cnn)\n",
    "\n",
    "# output layer -> corn, soybeans, other crop, other non-crop\n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[cnn_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_pa[:,:,:,:2*N_BANDS]], y=y_train, validation_data=([X_val_pa[:,:,:,:2*N_BANDS]], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## CNN branch only ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, N_TIMESTEPS*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(conv1_2d)\n",
    "mp1 = MaxPool2D()(conv2_2d)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(flat_cnn)\n",
    "\n",
    "# output layer -> corn, soybeans, other crop, other non-crop\n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[cnn_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_pa], y=y_train, validation_data=([X_val_pa], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSTM only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LSTM-only model does not have early-season classification because it does not make sense to use an LSTM with only one timestep (we only use the greenup DOY for early-season classification)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## LSTM only ########\n",
    "from keras.layers import BatchNormalization\n",
    "\n",
    "input_shape = (2, N_BANDS)\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(conv1)\n",
    "drop1 = Dropout(0.3)(lstm1)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "drop2 = Dropout(0.3)(lstm2)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "drop3 = Dropout(0.3)(lstm3)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "drop4 = Dropout(0.3)(lstm4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "#fully-connected layer\n",
    "fc = Dense(64, activation='relu')(flat_lstm)\n",
    "\n",
    "# output layer -> corn, soybeans, other crop, other non-crop\n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px[:,:2,:]], y=y_train, validation_data=([X_val_px[:,:2,:]], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## LSTM only ########\n",
    "from keras.layers import BatchNormalization\n",
    "\n",
    "input_shape = (N_TIMESTEPS, N_BANDS)\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(conv1)\n",
    "# bn2 = BatchNormalization()(lstm1)\n",
    "drop1 = Dropout(0.3)(lstm1)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "drop2 = Dropout(0.3)(lstm2)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "drop3 = Dropout(0.3)(lstm3)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "drop4 = Dropout(0.3)(lstm4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "#fully-connected layer\n",
    "fc = Dense(64, activation='relu')(flat_lstm)\n",
    "\n",
    "# output layer -> corn, soybeans, other crop, other non-crop\n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px], y=y_train, validation_data=([X_val_px], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN-LSTM-delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, there is no early-season classification because of the LSTM component."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season LSTM-CNN_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## LSTM branch ########\n",
    "from keras.layers import BatchNormalization\n",
    "\n",
    "input_shape = (2, N_BANDS)\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "bn1 = BatchNormalization()(conv1)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn1)\n",
    "bn2 = BatchNormalization()(lstm1)\n",
    "drop1 = Dropout(0.3)(bn2)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "bn3 = BatchNormalization()(lstm2)\n",
    "drop2 = Dropout(0.3)(bn3)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "bn4 = BatchNormalization()(lstm3)\n",
    "drop3 = Dropout(0.3)(bn4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "bn5 = BatchNormalization()(lstm4)\n",
    "drop4 = Dropout(0.3)(bn5)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "# Build the LSTM branch\n",
    "lstm = Model(inputs=lstm_inputs, outputs=flat_lstm)\n",
    "\n",
    "######## CNN branch ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 2*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "bn6 = BatchNormalization()(conv1_2d)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn6)\n",
    "bn7 = BatchNormalization()(conv2_2d)\n",
    "mp1 = MaxPool2D()(bn7)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "# Build the CNN branch\n",
    "cnn = Model(inputs=cnn_inputs, outputs=flat_cnn)\n",
    "\n",
    "########## Scalar timeline branch ###########\n",
    "scalar_inputs = Input(shape=(1,))\n",
    "\n",
    "combined = concatenate([lstm.output, cnn.output, scalar_inputs])\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(combined)\n",
    "\n",
    "# output layer -> corn, soybeans, other \n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm.input, cnn.input, scalar_inputs], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px[:,:2], X_train_pa[...,:2*N_BANDS], delta_train_p], y=y_train, validation_data=([X_val_px[:,:2], X_val_pa[...,:2*N_BANDS], delta_val_p], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season LSTM-CNN (no delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## LSTM branch ########\n",
    "from keras.layers import BatchNormalization\n",
    "\n",
    "input_shape = (2, N_BANDS)\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "bn1 = BatchNormalization()(conv1)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn1)\n",
    "bn2 = BatchNormalization()(lstm1)\n",
    "drop1 = Dropout(0.3)(bn2)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "bn3 = BatchNormalization()(lstm2)\n",
    "drop2 = Dropout(0.3)(bn3)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "bn4 = BatchNormalization()(lstm3)\n",
    "drop3 = Dropout(0.3)(bn4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "bn5 = BatchNormalization()(lstm4)\n",
    "drop4 = Dropout(0.3)(bn5)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "# Build the LSTM branch\n",
    "lstm = Model(inputs=lstm_inputs, outputs=flat_lstm)\n",
    "\n",
    "######## CNN branch ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 2*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "bn6 = BatchNormalization()(conv1_2d)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn6)\n",
    "bn7 = BatchNormalization()(conv2_2d)\n",
    "mp1 = MaxPool2D()(bn7)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "# Build the CNN branch\n",
    "cnn = Model(inputs=cnn_inputs, outputs=flat_cnn)\n",
    "\n",
    "combined = concatenate([lstm.output, cnn.output])\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(combined)\n",
    "\n",
    "# output layer -> corn, soybeans, other \n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm.input, cnn.input], outputs=preds)\n",
    "\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px[:,:2], X_train_pa[...,:2*N_BANDS]], y=y_train, validation_data=([X_val_px[:,:2], X_val_pa[...,:2*N_BANDS]], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season LSTM-CNN_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "######## LSTM branch ########\n",
    "from keras.layers import BatchNormalization\n",
    "np.random.seed(42)\n",
    "\n",
    "input_shape = (3, X_train_px.shape[-1])\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "bn1 = BatchNormalization()(conv1)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn1)\n",
    "bn2 = BatchNormalization()(lstm1)\n",
    "drop1 = Dropout(0.3)(bn2)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "bn3 = BatchNormalization()(lstm2)\n",
    "drop2 = Dropout(0.3)(bn3)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "bn4 = BatchNormalization()(lstm3)\n",
    "drop3 = Dropout(0.3)(bn4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "bn5 = BatchNormalization()(lstm4)\n",
    "drop4 = Dropout(0.3)(bn5)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "# Build the LSTM branch\n",
    "lstm = Model(inputs=lstm_inputs, outputs=flat_lstm)\n",
    "\n",
    "######## CNN branch ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 3*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "bn6 = BatchNormalization()(conv1_2d)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn6)\n",
    "bn7 = BatchNormalization()(conv2_2d)\n",
    "mp1 = MaxPool2D()(bn7)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "# Build the CNN branch\n",
    "cnn = Model(inputs=cnn_inputs, outputs=flat_cnn)\n",
    "\n",
    "########## Scalar timeline branch ###########\n",
    "scalar_inputs = Input(shape=(1,))\n",
    "\n",
    "combined = concatenate([lstm.output, cnn.output, scalar_inputs])\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(combined)\n",
    "\n",
    "# output layer -> corn, soybeans, other \n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm.input, cnn.input, scalar_inputs], outputs=preds)\n",
    "np.random.seed(42)\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px, X_train_pa, delta_train_s], y=y_train, validation_data=([X_val_px, X_val_pa, delta_val_s], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Late-season LSTM-CNN (no delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######## LSTM branch ########\n",
    "from keras.layers import BatchNormalization\n",
    "\n",
    "input_shape = (3, N_BANDS)\n",
    "lstm_inputs = Input(shape=input_shape)\n",
    "\n",
    "conv1 = Conv1D(filters=8, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_inputs)\n",
    "bn1 = BatchNormalization()(conv1)\n",
    "\n",
    "lstm1 = LSTM(16, input_shape=conv1.shape, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn1)\n",
    "bn2 = BatchNormalization()(lstm1)\n",
    "drop1 = Dropout(0.3)(bn2)\n",
    "\n",
    "lstm2 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(drop1)\n",
    "bn3 = BatchNormalization()(lstm2)\n",
    "drop2 = Dropout(0.3)(bn3)\n",
    "\n",
    "# Concatenate the output of LSTM1 and LSTM2\n",
    "lstm_dense_1 = concatenate([drop1, drop2])\n",
    "\n",
    "lstm3 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_1)\n",
    "bn4 = BatchNormalization()(lstm3)\n",
    "drop3 = Dropout(0.3)(bn4)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, and LSTM3\n",
    "lstm_dense_2 = concatenate([drop1, drop2, drop3])\n",
    "\n",
    "lstm4 = LSTM(16, return_sequences=True, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(lstm_dense_2)\n",
    "bn5 = BatchNormalization()(lstm4)\n",
    "drop4 = Dropout(0.3)(bn5)\n",
    "\n",
    "# Concatenate the output of LSTM1, LSTM2, LSTM3, and LSTM4\n",
    "lstm_dense_3 = concatenate([drop1, drop2, drop3, drop4])\n",
    "\n",
    "flat_lstm = Flatten()(lstm_dense_3)\n",
    "\n",
    "# Build the LSTM branch\n",
    "lstm = Model(inputs=lstm_inputs, outputs=flat_lstm)\n",
    "\n",
    "######## CNN branch ########\n",
    "\n",
    "cnn_inputs = Input(shape=(5, 5, 3*N_BANDS))\n",
    "\n",
    "conv1_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(cnn_inputs)\n",
    "bn6 = BatchNormalization()(conv1_2d)\n",
    "\n",
    "conv2_2d = Conv2D(filters=64, kernel_size=3, padding='same', activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(bn6)\n",
    "bn7 = BatchNormalization()(conv2_2d)\n",
    "mp1 = MaxPool2D()(bn7)\n",
    "\n",
    "flat_cnn = Flatten()(mp1)\n",
    "\n",
    "# Build the CNN branch\n",
    "cnn = Model(inputs=cnn_inputs, outputs=flat_cnn)\n",
    "\n",
    "combined = concatenate([lstm.output, cnn.output])\n",
    "\n",
    "fc = Dense(64, activation='relu', kernel_initializer=keras.initializers.glorot_uniform(seed=4))(combined)\n",
    "\n",
    "# output layer -> corn, soybeans, other \n",
    "preds = Dense(N_CLASSES, activation='softmax')(fc)\n",
    "\n",
    "# Build the model\n",
    "model = Model(inputs=[lstm.input, cnn.input], outputs=preds)\n",
    "np.random.seed(42)\n",
    "# Compile the model with loss function and optimizer\n",
    "adam = Adam()\n",
    "model.compile(loss='categorical_crossentropy', optimizer=adam, metrics=['accuracy'])\n",
    "\n",
    "model.build(input_shape) \n",
    "\n",
    "# Train the model\n",
    "model.fit(x=[X_train_px, X_train_pa], y=y_train, validation_data=([X_val_px, X_val_pa], y_val), epochs=N_EPOCHS, batch_size=BATCH_SIZE, callbacks=[plot])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save('results/models/cnn_lstm_train=16tbl201718_val=16tbl201718_epochs=%d_batchsize=%d' % (N_EPOCHS, BATCH_SIZE))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Test the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = load_model('results/models/cnn_lstm_train=16tbl201718_val=16tbl201718_epochs=%d_batchsize=%d' % (N_EPOCHS, BATCH_SIZE))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test, y_labels = load_data(x_path='data/inputs_3doy/pheno_timeseries_16TBL_2019.npy', y_path='data/cdl_labels_16TBL_2019.npy', flatten=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dois_2019 = np.load('data/phenology/16TBL_2019_phenology.npy').astype(np.int16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_test_s = dois_2019[...,2]-dois_2019[...,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_test_p = dois_2019[...,1]-dois_2019[...,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_px, X_test_pa, y_test, delta_test_s, delta_test_p = make_inputs(X_test, y_labels, delta_test_s, delta_test_p, subset=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_px = np.swapaxes(X_test_px, 1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test = to_categorical(y_test, num_classes=N_CLASSES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardize the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_px[:,0] = gscaler.transform(X_test_px[:,0])\n",
    "X_test_px[:,1] = pscaler.transform(X_test_px[:,1])\n",
    "X_test_px[:,2] = sscaler.transform(X_test_px[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_pa = standardize(X_test_pa, mu_pa, std_pa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Early-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_pa[...,:1*N_BANDS]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_pa[...,:2*N_BANDS]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_pa])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSTM only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px[:,:2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN-LSTM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px[:,:2], X_test_pa[...,:2*N_BANDS]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px, X_test_pa])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN-LSTM-delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mid-season LSTM-CNN_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px[:,:2], X_test_pa[...,:2*N_BANDS], delta_test_p])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Late-season LSTM-CNN_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = model.predict([X_test_px, X_test_pa, delta_test_s])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute accuracy metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_true_class = np.argmax(y_test, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_true_class = np.reshape(y_true_class, [HEIGHT, WIDTH])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_class = np.argmax(pred, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_class = np.reshape(y_pred_class, [HEIGHT, WIDTH])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Test accuracy: %f\" % metrics.accuracy_score(y_true_class[3:-3,3:-3].flatten(), y_pred_class[3:-3,3:-3].flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confusion matrix (un-normalized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(metrics.confusion_matrix(y_true_class[3:-3,3:-3].flatten(), y_pred_class[3:-3,3:-3].flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalized confusion matrix by row (recall/producers accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def confusion_matrix_recall(y_true, y_pred):\n",
    "    conf = metrics.confusion_matrix(y_true.flatten(), y_pred.flatten())\n",
    "    _conf = np.zeros(conf.shape)\n",
    "    for row in range(conf.shape[0]): # each row is for one class\n",
    "        total = np.sum(conf[row])\n",
    "        _conf[row] = conf[row]/float(total)\n",
    "    return _conf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(confusion_matrix_recall(y_true_class[3:-3,3:-3].flatten(), y_pred_class[3:-3,3:-3].flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalized confusion matrix by column (precision/producer's accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def confusion_matrix_precision(y_true, y_pred):\n",
    "    conf = metrics.confusion_matrix(y_true.flatten(), y_pred.flatten())\n",
    "    _conf = np.zeros(conf.shape)\n",
    "    for col in range(conf.shape[1]): # each row is for one class\n",
    "        total = np.sum(conf[:,col])\n",
    "        _conf[:,col] = conf[:,col]/float(total)\n",
    "    return _conf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(confusion_matrix_precision(y_true_class[3:-3,3:-3].flatten(), y_pred_class[3:-3,3:-3].flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the results compared to CDL labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def class_to_rgb(Y):\n",
    "    Y_rgb = np.ndarray([Y.shape[0], Y.shape[1], 3])\n",
    "    Y_rgb[np.where(Y==CORN)] = [1, 0.82, 0] # yellow\n",
    "    Y_rgb[np.where(Y==OTHER)] = [0.8, 0.8, 0.8] # gray\n",
    "    Y_rgb[np.where(Y==SOYBEAN)] = [0.149, 0.439, 0] # green\n",
    "    return Y_rgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(class_to_rgb(y_true_class))\n",
    "plt.title('CDL 2019')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(class_to_rgb(y_true_class[500:1000,500:1000]))\n",
    "plt.title('CDL 2019')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(class_to_rgb(y_pred_class[500:1000,500:1000]))\n",
    "plt.title('Predicted 2019')\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(class_to_rgb(y_pred_class))\n",
    "plt.title('Predicted 2019')\n",
    "plt.axis('off')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}