notebooks/train_spectrogram_classifier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train Spectrogram Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import date\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.utils import shuffle\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "parent_dir = os.path.split(os.getcwd())[0]\n",
    "if parent_dir not in sys.path:\n",
    "    sys.path.insert(0, parent_dir)\n",
    "\n",
    "\n",
    "from scripts.viz_tools import stretch_histogram, normalize\n",
    "\n",
    "from scripts import dl_utils\n",
    "from scripts.dl_utils import predict_spectrogram, rect_from_point\n",
    "from scripts.nn_predict import make_predictions, visualize_predictions\n",
    "\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a Training Dataset\n",
    "Outputs will be: `x_train`, `y_train`, `x_test`, `y_test`, and optionally, `x_holdout`, `y_holdout`. Holdout data is only positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data_dir = '../data/training_data/pixel_arrays_3mo-mosaics_2x-int/'\n",
    "MOSAIC_PERIOD = 3\n",
    "SPECTROGRAM_INTERVAL = 2\n",
    "SPECTROGRAM_STEPS = 2 # WIP: Currently this is fixed - only pairs of mosaics are supported.\n",
    "\n",
    "data_files = [\n",
    "    'tpa_polygons_2019-01-01_2021-06-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'v_1.1.5_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'city_points_30_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'bali_bare_earth_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'bali_bootstrap_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'maluku_v0.0.7_pixel_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'papua_v0.0.7_pixel_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'java_v0.0.7_pixel_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'kalimantan_v0.0.7_pixel_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'sumatra_v0.0.7_pixel_negatives_2019-01-01_2021-01-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'v0.0.7_indonesia_confirmed_2019_polygons_2019-01-01_2021-06-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'sumatra_v0.0.7_pixel_negative_polygons_2019-01-01_2021-06-01_period_3_interval_2_method_min_pixel_arrays.pkl',\n",
    "    'java_v0.0.7_positive_polygons_2019_v2_2019-01-01_2021-06-01_period_3_interval_2_method_min_pixel_arrays.pkl'\n",
    "             ]\n",
    "label_files = [f.split('s.pkl')[0] + '_labels.pkl' for f in data_files]\n",
    "\n",
    "# Specify a path to a holdout set of positive-class-only pixel arrays. Or leave an empty string.\n",
    "holdout_file = ''\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_arrays = []\n",
    "labels = []\n",
    "for data, label in tqdm(zip(data_files, label_files), total=len(data_files)):\n",
    "    with open(os.path.join(train_data_dir, data), 'rb') as f:\n",
    "        pixel_arrays += pickle.load(f)\n",
    "    with open(os.path.join(train_data_dir, label), 'rb') as f:\n",
    "        labels += pickle.load(f)\n",
    "            \n",
    "pixel_arrays = np.array(pixel_arrays)\n",
    "labels = np.array(labels)\n",
    "positive_arrays = pixel_arrays[labels == 1]\n",
    "negative_arrays = pixel_arrays[labels == 0]\n",
    "\n",
    "print(f\"Loaded {len(positive_arrays):,} positive pixel arrays and {len(negative_arrays):,} negative pixel arrays\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if holdout_file:\n",
    "    with open(os.path.join(train_data_dir, holdout_file), 'rb') as f:\n",
    "        holdout_pixel_arrays = np.array(pickle.load(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter positive samples such that NDVI is within a range\n",
    "This is useful since the positive patches can include surrounding vegetation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_ndvi(pixel_arrays):\n",
    "    return (pixel_arrays[:,7] - pixel_arrays[:,3]) / (pixel_arrays[:,7] + pixel_arrays[:,3])\n",
    "\n",
    "def filter_ndvi(data, lower_bound=0, upper_bound=0.4):\n",
    "    ndvi = compute_ndvi(data)\n",
    "    index = np.logical_and(ndvi > lower_bound, ndvi < upper_bound)\n",
    "    filtered_data = data[index.all(axis=1)]\n",
    "    print(f\"{len(filtered_data) / len(data):.1%} of samples within NDVI range\")\n",
    "    return filtered_data\n",
    "\n",
    "def filter_bright(data, brightness_threshold=2500):\n",
    "    filtered_data = data[np.mean(data, axis=(1,2)) < brightness_threshold]\n",
    "    filtered_data.shape\n",
    "    print(f\"{len(filtered_data) / len(data) :.1%} of data below brightness limit\")\n",
    "    return filtered_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_positive_arrays = filter_ndvi(positive_arrays)\n",
    "filtered_positive_arrays = filter_bright(filtered_positive_arrays)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combine data and create train test split\n",
    "Also expand dimensions to account for batches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.concatenate((filtered_positive_arrays, negative_arrays))\n",
    "y = np.concatenate((np.ones(len(filtered_positive_arrays)), np.zeros(len(negative_arrays))))\n",
    "\n",
    "x, y = shuffle(x, y, random_state=42)\n",
    "x = normalize(x)\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.20, random_state=42)\n",
    "x_train = np.expand_dims(x_train, -1)\n",
    "x_test = np.expand_dims(x_test, -1)\n",
    "\n",
    "print(\"Num Train Samples:\\t\\t\", len(x_train))\n",
    "print(\"Num Test Samples:\\t\\t\", len(x_test))\n",
    "print(f\"Percent Negative Train:\\t {sum(y_train == 0.0) / len(y_train):.1%}\")\n",
    "print(f\"Percent Negative Test:\\t {sum(y_test == 0.0) / len(y_test):.1%}\")\n",
    "print(f\"Input data shape: {x_train.shape}\")\n",
    "\n",
    "# Note: I am accustomed to assigning two classes for binary classification. \n",
    "# This habit comes from an issue in theano a long time ago, but I'm too superstitious to change it.\n",
    "num_classes = 2\n",
    "y_train = keras.utils.to_categorical(y_train, num_classes)\n",
    "y_test = keras.utils.to_categorical(y_test, num_classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if holdout_file:\n",
    "    filtered_holdout_arrays = filter_ndvi(holdout_pixel_arrays)\n",
    "    filtered_holdout_arrays = filter_bright(filtered_holdout_arrays)\n",
    "    x_holdout = normalize(filtered_holdout_arrays)\n",
    "    x_holdout = np.expand_dims(x_holdout, -1)\n",
    "    y_holdout = keras.utils.to_categorical(np.ones(len(x_holdout)), num_classes)\n",
    "    print(\"Num Holdout Samples:\\t\\t\", len(x_holdout))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create and Train a Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_shape = np.shape(x_train[0])\n",
    "model = keras.Sequential(\n",
    "    [\n",
    "        keras.Input(shape=input_shape),\n",
    "        layers.Conv2D(16, kernel_size=(3,2), activation=\"relu\"),\n",
    "        layers.Conv2D(32, kernel_size=(3,1), activation=\"relu\"),\n",
    "        layers.Flatten(),\n",
    "        layers.Dense(32, activation=\"relu\"),\n",
    "        layers.Dense(32, activation=\"relu\"),\n",
    "        layers.Dropout(0.1),\n",
    "        layers.Dense(32, activation=\"relu\"),\n",
    "        layers.Dense(num_classes, activation=\"softmax\"),\n",
    "    ]\n",
    ")\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional Class Weighting\n",
    "Over experimental testing, I found that weighting classes seemed to degrade performance. This could use further investigation. Convert cell type from Raw NBConvert back to Code to run."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "from sklearn.utils import class_weight\n",
    "negative_weight, positive_weight = class_weight.compute_class_weight('balanced', \n",
    "                                                                     classes = np.unique(y_train),\n",
    "                                                                     y = y_train[:,1])\n",
    "print(f\"Negative Weight: {negative_weight:.2f}\")\n",
    "print(f\"Positive Weight: {positive_weight:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile model. Note that many of these metrics are extraneous. \n",
    "# Can be useful to track during training at times though\n",
    "model.compile(loss=\"binary_crossentropy\", \n",
    "              optimizer=\"adam\", \n",
    "              metrics=[keras.metrics.Recall(thresholds=(0.7), name='precision'), \n",
    "                       keras.metrics.Precision(thresholds=(0.7), name='recall'),\n",
    "                       keras.metrics.AUC(curve='PR', name='auc'),\n",
    "                       \"accuracy\"],\n",
    "              #loss_weights = sum(y_train) / len(y_train),\n",
    "              #weighted_metrics = ['accuracy']\n",
    "             )\n",
    "\n",
    "train_accuracy = []\n",
    "test_accuracy = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "batch_size = 128\n",
    "epochs = 15\n",
    "\n",
    "model.fit(x_train, \n",
    "          y_train, \n",
    "          batch_size=batch_size, \n",
    "          epochs=epochs, \n",
    "          validation_data=(x_test, y_test),\n",
    "          #validation_split=0.1,\n",
    "          #class_weight = {0: negative_weight, 1: positive_weight}\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_accuracy += model.history.history['accuracy']\n",
    "test_accuracy += model.history.history['val_accuracy']\n",
    "plt.figure(figsize=(8,5), dpi=100, facecolor=(1,1,1))\n",
    "plt.plot(train_accuracy, label='Train Acc')\n",
    "plt.plot(test_accuracy, c='r', label='Val Acc')\n",
    "percent_negative = (sum(y_train == 0.0) / len(y_train))[1]\n",
    "plt.plot([0, epochs-1], [percent_negative, percent_negative], '--', c='gray', label='Baseline')\n",
    "plt.grid()\n",
    "plt.xlabel('Epoch')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.legend()\n",
    "plt.title('Network Train and Val Accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.8\n",
    "print(\"Test Set Metrics:\")\n",
    "print(classification_report(y_test[:,1], model.predict(x_test)[:,1] > threshold, \n",
    "                            target_names=['No TPA', 'TPA']))\n",
    "\n",
    "if holdout_file:\n",
    "    print(\"\\nHoldout Positive Set Metrics:\")\n",
    "    print(classification_report(y_holdout[:,1], model.predict(x_holdout)[:,1] > threshold, \n",
    "                                target_names=['No TPA', 'TPA']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "version_number = '0.0.4'\n",
    "\n",
    "current_date = date.today()\n",
    "model_name = f\"spectrogram_v{version_number}_{current_date.isoformat()}\"\n",
    "assert not os.path.exists('../models/' + model_name + '.h5'), f\"Model of name {model_name} already exists\"\n",
    "\n",
    "with open('../models/' + model_name + '_config.txt', 'w') as f:\n",
    "    f.write(f'Image mosaic period: {MOSAIC_PERIOD}\\n')\n",
    "    f.write(f'Spectrogram interval, steps: {SPECTROGRAM_INTERVAL}, {SPECTROGRAM_STEPS}\\n')\n",
    "    f.write('Input Data:\\n')\n",
    "    [f.write(file + '\\n') for file in data_files]\n",
    "    f.write(f\"\\nBatch Size: {batch_size}\")\n",
    "    f.write(f\"\\nTraining Epochs: {len(train_accuracy)}\")\n",
    "    f.write('\\n\\nClassification Report\\n')\n",
    "    f.write(classification_report(y_test[:,1], model.predict(x_test)[:,1] > threshold, \n",
    "                            target_names=['No TPA', 'TPA']))\n",
    "model.save('../models/' + model_name + '.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Network Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment to restore a saved model.\n",
    "# model = keras.models.load_model('../models/spectrogram_v0.0.4_2021-05-19.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TPA Kebon Kongok on Lombok\n",
    "rect_width = 0.02\n",
    "coords = [116.0908,-8.6451]\n",
    "start_date = '2020-05-01'\n",
    "mosaic_end_date = '2021-05-01'\n",
    "\n",
    "mosaics, _ = dl_utils.download_mosaics(\n",
    "    rect_from_point(coords, rect_width), start_date, mosaic_end_date, MOSAIC_PERIOD)\n",
    "kebon = dl_utils.pair(mosaics, SPECTROGRAM_INTERVAL)\n",
    "print(f'{len(kebon)} pairs returned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds_kebon = [predict_spectrogram(pair, model) for pair in kebon]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TPA Alak, West Timor\n",
    "rect_width = 0.02\n",
    "coords = [123.546910, -10.200799]\n",
    "start_date = '2020-05-01'\n",
    "mosaic_end_date = '2021-05-01'\n",
    "mosaics, _ = dl_utils.download_mosaics(\n",
    "    rect_from_point(coords, rect_width), start_date, mosaic_end_date, MOSAIC_PERIOD)\n",
    "alak = dl_utils.pair(mosaics, SPECTROGRAM_INTERVAL)\n",
    "print(f'{len(alak)} pairs returned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds_alak = [predict_spectrogram(pair, model) for pair in alak]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TPA Lelowai, West Timor, large box\n",
    "rect_width = 0.04\n",
    "coords = [124.892853, -9.196869]\n",
    "start_date = '2020-05-01'\n",
    "mosaic_end_date = '2021-05-01'\n",
    "mosaics, _ = dl_utils.download_mosaics(\n",
    "    rect_from_point(coords, rect_width), start_date, mosaic_end_date, MOSAIC_PERIOD)\n",
    "lelowai = dl_utils.pair(mosaics, SPECTROGRAM_INTERVAL)\n",
    "print(f'{len(lelowai)} pairs returned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds_lelowai = [predict_spectrogram(pair, model) for pair in lelowai]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# single pair predictions\n",
    "visualize_predictions(kebon[2], [preds_kebon[2]], threshold=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# averaging over pairs' predictions\n",
    "visualize_predictions([p[1] for p in kebon], preds_kebon, threshold=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show timeseries predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_series(pairs, preds, threshold=0.8, savepath=''):\n",
    "    fig, axes = plt.subplots(1, len(preds), figsize=(30, 30))\n",
    "    for ax, (pair, pred) in zip(axes, zip(pairs, preds)):\n",
    "        rgb = normalize(pair[1][:,:,3:0:-1])\n",
    "        rgb[pred > threshold] = [0.9, 0, 0.1]\n",
    "        ax.imshow(np.clip(rgb, 0, 1))\n",
    "        ax.axis('off')\n",
    "    if savepath:\n",
    "        plt.savefig(savepath, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series(kebon, preds_kebon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series(alak, preds_alak)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series(lelowai, preds_lelowai)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline comparisons\n",
    "Specify an older non-spectrogram model for comparison. (It is run on the second image of each pair.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baseline_model_path = '../models/model_v1.1.9_baseline_5-15-2021.h5'\n",
    "baseline_model = keras.models.load_model(baseline_model_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds_kebon_baseline = make_predictions([p[1] for p in kebon], baseline_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# averaging over the predictions\n",
    "visualize_predictions([p[1] for p in kebon], preds_kebon_baseline, threshold=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series(kebon, preds_kebon_baseline)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('m1-plastics')",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "d425afa4a959a86aa036beaa1a58ff3469f38e31f3ec97f5785c695b9108eced"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}