1171_paper_LAS.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_SHAc5qbiR8l"
      },
      "source": [
        "# Introduction\n",
        "\n",
        "Classify roads from planetscope 3m imagery, with FCNN U-Net. Based on GEE FCNN notebook.\n",
        "This notebook shows:\n",
        "\n",
        "1.   Exporting training/testing patches from Earth Engine, suitable for training an FCNN model.\n",
        "2.   Preprocessing.\n",
        "3.   Training and validating an FCNN model.\n",
        "4.   Making predictions with the trained model and importing them to Earth Engine."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ir8E4Cghtnl_"
      },
      "source": [
        "## Version\n",
        "\n",
        "v5\n",
        "for all Attica (splitted in tiles)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yuO7BD7ee4CS"
      },
      "source": [
        "Create more RAM for colab by crushing it"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_MJ4kW1pEhwP"
      },
      "source": [
        "# Libraries & imports\n",
        "\n",
        "Authenticate and import as necessary."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "neIa46CpciXq"
      },
      "source": [
        "# Cloud authentication.\n",
        "from google.colab import auth\n",
        "auth.authenticate_user()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jat01FEoUMqg"
      },
      "source": [
        "# Import, authenticate and initialize the Earth Engine library.\n",
        "import ee\n",
        "ee.Authenticate()\n",
        "ee.Initialize()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8RnZzcYhcpsQ"
      },
      "source": [
        "import tensorflow as tf\n",
        "print(\"Temsorflow version:\")\n",
        "print(tf.__version__)\n",
        "\n",
        "import folium\n",
        "print(\"Folium version:\")\n",
        "print(folium.__version__)\n",
        "\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "import os\n",
        "import numpy as np"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iT8ycmzClYwf"
      },
      "source": [
        "# Variables\n",
        "\n",
        "Declare the variables that will be in use throughout the notebook."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qKs6HuxOzjMl"
      },
      "source": [
        "## Specify your Cloud Storage Bucket\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "obDDH1eDzsch"
      },
      "source": [
        "# INSERT YOUR BUCKET HERE:\n",
        "BUCKET = 'n-verde_bucket'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wmfKLl9XcnGJ"
      },
      "source": [
        "## Set other global variables"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "psz7wJKalaoj"
      },
      "source": [
        "# Specify names locations for outputs in Cloud Storage.\n",
        "FOLDER =  '1171-LAS-256_16000'\n",
        "TRAINING_BASE = 'training_patches'\n",
        "EVAL_BASE = 'eval_patches'\n",
        "\n",
        "# Specify inputs (Planetscope bands) to the model and the response variable.\n",
        "opticalBands = ['b1', 'b2', 'b3', 'b4']\n",
        "BANDS = opticalBands\n",
        "RESPONSE = 'roads'\n",
        "FEATURES = BANDS + [RESPONSE]\n",
        "\n",
        "# Specify the number of tiles the imagery will be split in\n",
        "# predictions will run per tile\n",
        "# will be split in SxS tiles\n",
        "S = 8\n",
        "\n",
        "# Specify the size and shape of patches expected by the model.\n",
        "KERNEL_SIZE =  256\n",
        "KERNEL_SHAPE = [KERNEL_SIZE, KERNEL_SIZE]\n",
        "COLUMNS = [\n",
        "  tf.io.FixedLenFeature(shape=KERNEL_SHAPE, dtype=tf.float32) for k in FEATURES\n",
        "]\n",
        "FEATURES_DICT = dict(zip(FEATURES, COLUMNS))\n",
        "\n",
        "# ---------------------- SAMPLE SIZE --------------------------\n",
        "# Sizes of the training and evaluation datasets.\n",
        "TRAIN_SIZE = 1500    # 16000\n",
        "EVAL_SIZE = 643      # 4000\n",
        "\n",
        "# Specify model training parameters.\n",
        "# https://www.kite.com/python/docs/keras.backend.moving_averages.distribution_strategy_context.distribute_lib.dataset_ops.BatchDataset.shuffle\n",
        "# BATCH_SIZE is dependent on your GPU memory. (e.g. on my PC it can't be larger than 4, on Colab it can be 32)\n",
        "BATCH_SIZE = 32\n",
        "# https://machinelearningmastery.com/how-to-stop-training-deep-neural-networks-at-the-right-time-using-early-stopping/\n",
        "EPOCHS = 50\n",
        "# For perfect shuffling, set the buffer size equal to the full size of the dataset.\n",
        "# https://www.tensorflow.org/api_docs/python/tf/data/experimental/shuffle_and_repeat\n",
        "BUFFER_SIZE = TRAIN_SIZE\n",
        "\n",
        "# ---------------------- OPTIMIZER --------------------------\n",
        "\n",
        "OPTIMIZER = tf.keras.optimizers.Adam()\n",
        "\n",
        "# ---------------------- LOSS --------------------------\n",
        "\n",
        "# LOSS = tf.keras.losses.get('binary_crossentropy')\n",
        "\n",
        "# DICE LOSS\n",
        "def dice_coef(y_true, y_pred, smooth=1):\n",
        "    y_true_f = tf.keras.backend.flatten(y_true)\n",
        "    y_pred_f = tf.keras.backend.flatten(y_pred)\n",
        "    intersection = tf.keras.backend.sum(y_true_f * y_pred_f)\n",
        "    return (2. * intersection + smooth) / (tf.keras.backend.sum(y_true_f) + tf.keras.backend.sum(y_pred_f) + smooth)\n",
        "def dice_coef_loss(y_true, y_pred):\n",
        "    return -dice_coef(y_true, y_pred)\n",
        "LOSS = dice_coef_loss\n",
        "\n",
        "# ---------------------- METRICS --------------------------\n",
        "# https://keras.io/api/metrics/accuracy_metrics/\n",
        "# METRICS = [metrics.get('binary_accuracy')]\n",
        "\n",
        "def dice_coef(y_true, y_pred, smooth=1):\n",
        "    y_true_f = tf.keras.backend.flatten(y_true)\n",
        "    y_pred_f = tf.keras.backend.flatten(y_pred)\n",
        "    intersection = tf.keras.backend.sum(y_true_f * y_pred_f)\n",
        "    return (2. * intersection + smooth) / (tf.keras.backend.sum(y_true_f) + tf.keras.backend.sum(y_pred_f) + smooth)\n",
        "\n",
        "METRICS = [tf.keras.metrics.get('binary_accuracy'),\n",
        "           dice_coef,\n",
        "           tf.keras.metrics.MeanIoU(num_classes=2),\n",
        "           tf.keras.metrics.Recall()]\n",
        "\n",
        "print('OK!')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hgoDc7Hilfc4"
      },
      "source": [
        "# Imagery\n",
        "\n",
        "Gather and setup the imagery to use for inputs (predictors) (Planetscope mosaic).  Display it in the notebook for a sanity check.\n",
        "\n",
        "Prepare the response (what we want to predict). This is the roads rasterized. Display to check."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-IlgXu-vcUEY"
      },
      "source": [
        "# Use folium to visualize the imagery.\n",
        "map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "\n",
        "# planetscope mosiac image of attica\n",
        "image = ee.Image('users/n-verde/PhD_1171/attica_feather_mask')\n",
        "\n",
        "mapid = image.getMapId({'bands': ['b4', 'b3', 'b2'], 'min': 730, 'max': 3500})\n",
        "folium.TileLayer(\n",
        "    tiles=mapid['tile_fetcher'].url_format,\n",
        "    attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "    overlay=True,\n",
        "    name='planetscope',\n",
        "  ).add_to(map)\n",
        "\n",
        "# roads ground truth\n",
        "gt = ee.Image('users/n-verde/PhD_1171/200_400m_OSM_LAS_edit_RASTER_UINT16')\n",
        "gt = gt.rename('roads')\n",
        "\n",
        "mapid = gt.getMapId({'bands': ['roads'], 'min': 0, 'max': 1})\n",
        "folium.TileLayer(\n",
        "    tiles=mapid['tile_fetcher'].url_format,\n",
        "    attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "    overlay=True,\n",
        "    name='ground_truth',\n",
        "  ).add_to(map)\n",
        "\n",
        "map.add_child(folium.LayerControl())\n",
        "\n",
        "map"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3WFrZOjtLuwl"
      },
      "source": [
        "## Preprocessing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pA-ljY77bihF"
      },
      "source": [
        "### Normalization\n",
        "\n",
        "Find the min and max of the planetscope image for rescaling to [0-1]. When normalizing data to [0-1] the training of the NN is made easier.\n",
        "If you are running this section for the first time, uncomment the reducers and comment-out the `min=ee.Number(...)` and `max=ee.Number(...)`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fLJVrCkpbhdW"
      },
      "source": [
        "# define the study area (area of planetscope image)\n",
        "aoi = ee.FeatureCollection('users/n-verde/PhD_1171/Attica_urbanExtent')\n",
        "\n",
        "print(\"finding min and max of planetscope bands, in order to normalize them...\")\n",
        "print(\"----------\")\n",
        "\n",
        "def findMin(imgBand, bandName):\n",
        "  min = ee.Number(imgBand.reduceRegion(**{\n",
        "  'reducer': ee.Reducer.min(),\n",
        "  'geometry': aoi,\n",
        "  'scale': 9,\n",
        "  'tileScale': 4,\n",
        "  # 'bestEffort': true,\n",
        "  'maxPixels': 1e9\n",
        "  }).get(bandName))\n",
        "\n",
        "  return min\n",
        "\n",
        "def findMax(imgBand, bandName):\n",
        "  max = ee.Number(imgBand.reduceRegion(**{\n",
        "  'reducer': ee.Reducer.max(),\n",
        "  'geometry': aoi,\n",
        "  'scale': 9,\n",
        "  'tileScale': 4,\n",
        "  # 'bestEffort': true,\n",
        "  'maxPixels': 1e9\n",
        "  }).get(bandName))\n",
        "\n",
        "  return max\n",
        "\n",
        "# b1 ----------\n",
        "\n",
        "bandName = \"b1\"\n",
        "imgBand = image.select(bandName)\n",
        "\n",
        "# min = findMin(imgBand, bandName)\n",
        "min = ee.Number(181)\n",
        "\n",
        "# max = findMax(imgBand, bandName)\n",
        "max = ee.Number(7988)\n",
        "\n",
        "test = ee.Algorithms.If(min.lt(ee.Number(0)), 1, 0)\n",
        "\n",
        "if (test.getInfo()==1):\n",
        "  imgBand = imgBand.add(min.abs()) # turn negative values to possitive\n",
        "  max = max.add(min.abs()) # change new max\n",
        "\n",
        "min_b1 = min\n",
        "max_b1 = max\n",
        "\n",
        "print(\"min of b1:\", min_b1.getInfo())\n",
        "print(\"max of b1:\", max_b1.getInfo())\n",
        "\n",
        "# b2 ----------\n",
        "\n",
        "bandName = \"b2\"\n",
        "imgBand = image.select(bandName)\n",
        "\n",
        "# min = findMin(imgBand, bandName)\n",
        "min = ee.Number(309)\n",
        "\n",
        "# max = findMax(imgBand, bandName)\n",
        "max = ee.Number(9432)\n",
        "\n",
        "test = ee.Algorithms.If(min.lt(ee.Number(0)), 1, 0)\n",
        "\n",
        "if (test.getInfo()==1):\n",
        "  imgBand = imgBand.add(min.abs()) # turn negative values to possitive\n",
        "  max = max.add(min.abs()) # change new max\n",
        "\n",
        "min_b2 = min\n",
        "max_b2 = max\n",
        "\n",
        "print(\"min of b2:\", min_b2.getInfo())\n",
        "print(\"max of b2:\", max_b2.getInfo())\n",
        "\n",
        "# b3 ----------\n",
        "\n",
        "bandName = \"b3\"\n",
        "imgBand = image.select(bandName)\n",
        "\n",
        "# min = findMin(imgBand, bandName)\n",
        "min = ee.Number(238)\n",
        "\n",
        "# max = findMax(imgBand, bandName)\n",
        "max = ee.Number(11849)\n",
        "\n",
        "test = ee.Algorithms.If(min.lt(ee.Number(0)), 1, 0)\n",
        "\n",
        "if (test.getInfo()==1):\n",
        "  imgBand = imgBand.add(min.abs()) # turn negative values to possitive\n",
        "  max = max.add(min.abs()) # change new max\n",
        "\n",
        "min_b3 = min\n",
        "max_b3 = max\n",
        "\n",
        "print(\"min of b3:\", min_b3.getInfo())\n",
        "print(\"max of b3:\", max_b3.getInfo())\n",
        "\n",
        "# b4 ----------\n",
        "\n",
        "bandName = \"b4\"\n",
        "imgBand = image.select(bandName)\n",
        "\n",
        "# min = findMin(imgBand, bandName)\n",
        "min = ee.Number(1)\n",
        "\n",
        "# max = findMax(imgBand, bandName)\n",
        "max = ee.Number(11453)\n",
        "\n",
        "test = ee.Algorithms.If(min.lt(ee.Number(0)), 1, 0)\n",
        "\n",
        "if (test.getInfo()==1):\n",
        "  imgBand = imgBand.add(min.abs()) # turn negative values to possitive\n",
        "  max = max.add(min.abs()) # change new max\n",
        "\n",
        "min_b4 = min\n",
        "max_b4 = max\n",
        "\n",
        "print(\"min of b4:\", min_b4.getInfo())\n",
        "print(\"max of b4:\", max_b4.getInfo())\n",
        "\n",
        "\n",
        "print(\"----------\")\n",
        "print(\"done!\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M-_nPEgrbuwe"
      },
      "source": [
        "Rescale to [0-1] (normalize)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l5_JKDrIbxyk"
      },
      "source": [
        "def rescaleFixedMinMax(img, bandName, min, max):\n",
        "  Min = ee.Number(min)\n",
        "  Max = ee.Number(max)\n",
        "\n",
        "  imgBand = img.select(bandName)\n",
        "\n",
        "  imgBand = imgBand.float()\n",
        "\n",
        "  imgBandNorm = imgBand.divide(Max) # normalize the data to 0 - 1\n",
        "\n",
        "  return imgBandNorm\n",
        "\n",
        "# Apply the rescale function to all bands of the planetscope image\n",
        "# b1 ----------\n",
        "b1_res = rescaleFixedMinMax(image,\"b1\", min_b1, max_b1)\n",
        "\n",
        "# b2 ----------\n",
        "b2_res = rescaleFixedMinMax(image,\"b2\", min_b2, max_b2)\n",
        "\n",
        "# b3 ----------\n",
        "b3_res = rescaleFixedMinMax(image,\"b3\", min_b3, max_b3)\n",
        "\n",
        "# b4 ----------\n",
        "b4_res = rescaleFixedMinMax(image,\"b4\", min_b4, max_b4)\n",
        "\n",
        "# merge the bands to a single image\n",
        "image_res = b1_res.addBands(b2_res).addBands(b3_res).addBands(b4_res)\n",
        "\n",
        "# # display rescaled image\n",
        "# mapid = image_res.getMapId({'bands': ['b4', 'b3', 'b2'], 'min': 0, 'max': 0.4})\n",
        "# map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "# folium.TileLayer(\n",
        "#     tiles=mapid['tile_fetcher'].url_format,\n",
        "#     attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "#     overlay=True,\n",
        "#     name='planetscope_masked',\n",
        "#   ).add_to(map)\n",
        "# map.add_child(folium.LayerControl())\n",
        "# map"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tza7Ooasc-cz"
      },
      "source": [
        "### Masking\n",
        "\n",
        "Mask the rescaled planetscope image to the samples extent (boxes image)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sCL2d23oc_U5"
      },
      "source": [
        "boxes = ee.Image('users/n-verde/PhD_1171/200_400m_RASTER_boxes')\n",
        "masked_image = image_res.mask(boxes)\n",
        "\n",
        "# display the masked image\n",
        "mapid = masked_image.getMapId({'bands': ['b4', 'b3', 'b2'], 'min': 0, 'max': 0.4})\n",
        "map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "folium.TileLayer(\n",
        "    tiles=mapid['tile_fetcher'].url_format,\n",
        "    attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "    overlay=True,\n",
        "    name='planetscope_masked',\n",
        "  ).add_to(map)\n",
        "map.add_child(folium.LayerControl())\n",
        "map"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CTS7_ZzPDhhg"
      },
      "source": [
        "### Stack\n",
        "\n",
        "Stack the 2D images (rescaled planetscope image and road samples) to create a single image from which samples can be taken. Convert the image into an array image in which each pixel stores 256x256 patches of pixels for each band.  This is a key step that bears emphasis: to export training patches, convert a multi-band image to [an array image](https://developers.google.com/earth-engine/arrays_array_images#array-images) using [`neighborhoodToArray()`](https://developers.google.com/earth-engine/api_docs#eeimageneighborhoodtoarray), then sample the image at points."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eGHYsdAOipa4"
      },
      "source": [
        "featureStack = ee.Image.cat([\n",
        "  masked_image.select(BANDS),\n",
        "  gt.select(RESPONSE)\n",
        "]).float()\n",
        "\n",
        "list = ee.List.repeat(1, KERNEL_SIZE)\n",
        "lists = ee.List.repeat(list, KERNEL_SIZE)\n",
        "kernel = ee.Kernel.fixed(KERNEL_SIZE, KERNEL_SIZE, lists)\n",
        "\n",
        "arrays = featureStack.neighborhoodToArray(kernel)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zgCkAHZ5qM7T"
      },
      "source": [
        "### Split image to tiles"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S0IC24vfyyUu"
      },
      "source": [
        "Define a function for displaying Earth Engine image tiles on a folium map\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bDJ5-4pGyyqb"
      },
      "source": [
        "# Define a method for displaying Earth Engine image tiles on a folium map.\n",
        "def add_ee_layer(self, ee_object, vis_params, name):\n",
        "\n",
        "    try:\n",
        "        # display ee.Image()\n",
        "        if isinstance(ee_object, ee.image.Image):\n",
        "            map_id_dict = ee.Image(ee_object).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "            ).add_to(self)\n",
        "        # display ee.ImageCollection()\n",
        "        elif isinstance(ee_object, ee.imagecollection.ImageCollection):\n",
        "            ee_object_new = ee_object.mosaic()\n",
        "            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "            ).add_to(self)\n",
        "        # display ee.Geometry()\n",
        "        elif isinstance(ee_object, ee.geometry.Geometry):\n",
        "            folium.GeoJson(\n",
        "            data = ee_object.getInfo(),\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "        ).add_to(self)\n",
        "        # display ee.FeatureCollection()\n",
        "        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):\n",
        "            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n",
        "            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n",
        "            folium.raster_layers.TileLayer(\n",
        "            tiles = map_id_dict['tile_fetcher'].url_format,\n",
        "            attr = 'Google Earth Engine',\n",
        "            name = name,\n",
        "            overlay = True,\n",
        "            control = True\n",
        "        ).add_to(self)\n",
        "\n",
        "    except:\n",
        "        print(\"Could not display {}\".format(name))\n",
        "\n",
        "# Add EE drawing method to folium.\n",
        "folium.Map.add_ee_layer = add_ee_layer"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dlH3TvdJqSA3"
      },
      "source": [
        "# This function splits a geometry into equal-sized subrectangles.\n",
        "# splits a geometry to parts^2 subregions\n",
        "# Assumes that the geometry has only 4 vertices.\n",
        "# returns a list of polygons (you can map over it)\n",
        "def split(geom, nSplits):\n",
        "\n",
        "  one = ee.Number(1)\n",
        "\n",
        "  # Return (nSplit+1)^2 coordinate pairs, given 4 vertices.\n",
        "  def toPts(rect, nSplits):\n",
        "    k = ee.List.sequence(0, None, one.divide(nSplits), one.add(nSplits));\n",
        "\n",
        "    def f1(x):\n",
        "\n",
        "      def f2(y):\n",
        "        xp = one.subtract(x)\n",
        "        yp = one.subtract(y)\n",
        "        coeffs = ee.Array([[xp.multiply(yp), yp.multiply(x), xp.multiply(y), ee.Number(x).multiply(y)]])\n",
        "        return coeffs.matrixMultiply(rect).project([1])\n",
        "\n",
        "      return k.map(f2)\n",
        "\n",
        "    return k.map(f1).flatten()\n",
        "\n",
        "  # Return nSplit^2 polygons, given the list of vertices built by toPts().\n",
        "  def toRects(pts, nSplits):\n",
        "    offsets = ee.List([0, 1, one.add(nSplits).add(one), one.add(nSplits)])\n",
        "    k1 = ee.List.sequence(0, None, one.add(nSplits), nSplits)\n",
        "\n",
        "    def f3(i):\n",
        "      k2 = ee.List.sequence(i, None, 1, nSplits)\n",
        "\n",
        "      def f4(j):\n",
        "\n",
        "        def f5(offset):\n",
        "          return ee.Array(pts.get(ee.Number(j).add(offset))).toList()\n",
        "\n",
        "        return ee.Geometry.Polygon(offsets.map(f5))\n",
        "\n",
        "      return k2.map(f4)\n",
        "\n",
        "    return k1.map(f3).flatten()\n",
        "\n",
        "\n",
        "  # Get the 4 vertices.  Assumes that the scene geometry has only 4 vertices.\n",
        "  rect = ee.List(geom.coordinates().get(0))\n",
        "  rect = ee.Array([rect.get(0), rect.get(1), rect.get(3), rect.get(2)])\n",
        "\n",
        "  pts = toPts(rect, nSplits)\n",
        "  rects = toRects(pts, nSplits)\n",
        "\n",
        "  print('Parts that geometry was splitted to ', rects.getInfo())\n",
        "\n",
        "  return ee.List(rects)\n",
        "\n",
        "# get planetscope image bounds in a box\n",
        "i = image.geometry().bounds().transform('EPSG:4326',100);\n",
        "\n",
        "# split to SxS tiles\n",
        "splitted = split(i,S);\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DoDGO3lxz7o2"
      },
      "source": [
        "# show on Map\n",
        "\n",
        "no = 36\n",
        "spl = ee.Geometry(ee.List(splitted.get(no-1)))\n",
        "desc = 'tile' + str(int(no-1)) + 'from the split'\n",
        "\n",
        "map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "\n",
        "# planetscope mosiac image of attica\n",
        "image = ee.Image('users/n-verde/PhD_1171/attica_feather_mask')\n",
        "\n",
        "mapid = image.getMapId({'bands': ['b4', 'b3', 'b2'], 'min': 730, 'max': 3500})\n",
        "folium.TileLayer(\n",
        "    tiles=mapid['tile_fetcher'].url_format,\n",
        "    attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "    overlay=True,\n",
        "    name='planetscope',\n",
        "  ).add_to(map)\n",
        "\n",
        "# planetscope bounds\n",
        "map.add_ee_layer(i,{}, 'planetscope bounds')\n",
        "\n",
        "# tile\n",
        "map.add_ee_layer(spl, {}, desc)\n",
        "\n",
        "map.add_child(folium.LayerControl())\n",
        "map"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uNC20UFGMDST"
      },
      "source": [
        "# Sampling"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F4djSxBRG2el"
      },
      "source": [
        "Use some pre-made geometries to sample the stack in strategic locations.  Specifically, these are hand-made polygons in which to take the 256x256 samples.  Display the sampling polygons on a map, red for training polygons, blue for evaluation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ure_WaD0itQY"
      },
      "source": [
        "trainingPolys = ee.FeatureCollection('users/n-verde/PhD_1171/200_400m_TRAIN70')\n",
        "evalPolys = ee.FeatureCollection('users/n-verde/PhD_1171/200_400m_VAL30')\n",
        "\n",
        "polyImage = ee.Image(0).byte().paint(trainingPolys, 1).paint(evalPolys, 2)\n",
        "polyImage = polyImage.updateMask(polyImage)\n",
        "\n",
        "# display the sample boxes\n",
        "mapid = polyImage.getMapId({'min': 1, 'max': 2, 'palette': ['red', 'blue']})\n",
        "map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "folium.TileLayer(\n",
        "    tiles=mapid['tile_fetcher'].url_format,\n",
        "    attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
        "    overlay=True,\n",
        "    name='training & validation polygons',\n",
        "  ).add_to(map)\n",
        "map.add_child(folium.LayerControl())\n",
        "map"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZV890gPHeZqz"
      },
      "source": [
        "The mapped data look reasonable so take a sample from each polygon and merge the results into a single export.  The key step is sampling the array image at points, to get all the pixels in a 256x256 neighborhood at each point.  It's worth noting that to build the training and testing data for the FCNN, you export a single TFRecord file that contains patches of pixel values in each record.  You do NOT need to export each training/testing patch to a different image.  Since each record potentially contains a lot of data (especially with big patches or many input bands), some manual sharding of the computation is necessary to avoid the `computed value too large` error.  Specifically, the following code takes multiple (smaller) samples within each geometry, merging the results to get a single export.\n",
        "For whole Attica takes around **40m-1h** in GEE."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FyRpvwENxE-A",
        "cellView": "both"
      },
      "source": [
        "# Convert the feature collections to lists for iteration.\n",
        "trainingPolysList = trainingPolys.toList(trainingPolys.size())\n",
        "evalPolysList = evalPolys.toList(evalPolys.size())\n",
        "\n",
        "# These numbers determined experimentally.\n",
        "n = 10 # Number of shards in each polygon.\n",
        "N =  100 # Total sample size in each polygon.\n",
        "\n",
        "# Export all the training data (in many pieces), with one task\n",
        "# per geometry.\n",
        "for g in range(trainingPolys.size().getInfo()):\n",
        "  geomSample = ee.FeatureCollection([])\n",
        "  for i in range(n):\n",
        "    sample = arrays.sample(\n",
        "      region = ee.Feature(trainingPolysList.get(g)).geometry(),\n",
        "      scale = 1,\n",
        "      numPixels = N / n, # Size of the shard.\n",
        "      seed = i,\n",
        "      tileScale = 16\n",
        "    )\n",
        "    geomSample = geomSample.merge(sample)\n",
        "\n",
        "  desc = TRAINING_BASE + '_g' + str(g)\n",
        "  task = ee.batch.Export.table.toCloudStorage(\n",
        "    collection = geomSample,\n",
        "    description = desc,\n",
        "    bucket = BUCKET,\n",
        "    fileNamePrefix = FOLDER + '/' + desc,\n",
        "    fileFormat = 'TFRecord',\n",
        "    selectors = BANDS + [RESPONSE]\n",
        "  )\n",
        "  task.start()\n",
        "\n",
        "# Export all the evaluation data.\n",
        "for g in range(evalPolys.size().getInfo()):\n",
        "  geomSample = ee.FeatureCollection([])\n",
        "  for i in range(n):\n",
        "    sample = arrays.sample(\n",
        "      region = ee.Feature(evalPolysList.get(g)).geometry(),\n",
        "      scale = 1,\n",
        "      numPixels = N / n,\n",
        "      seed = i,\n",
        "      tileScale = 16\n",
        "    )\n",
        "    geomSample = geomSample.merge(sample)\n",
        "\n",
        "  desc = EVAL_BASE + '_g' + str(g)\n",
        "  task = ee.batch.Export.table.toCloudStorage(\n",
        "    collection = geomSample,\n",
        "    description = desc,\n",
        "    bucket = BUCKET,\n",
        "    fileNamePrefix = FOLDER + '/' + desc,\n",
        "    fileFormat = 'TFRecord',\n",
        "    selectors = BANDS + [RESPONSE]\n",
        "  )\n",
        "  task.start()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rWXrvBE4607G"
      },
      "source": [
        "## Training data\n",
        "\n",
        "Load the data exported from Earth Engine into a `tf.data.Dataset`.  The following are helper functions for that."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WWZ0UXCVMyJP"
      },
      "source": [
        "def parse_tfrecord(example_proto):\n",
        "  \"\"\"The parsing function.\n",
        "  Read a serialized example into the structure defined by FEATURES_DICT.\n",
        "  Args:\n",
        "    example_proto: a serialized Example.\n",
        "  Returns:\n",
        "    A dictionary of tensors, keyed by feature name.\n",
        "  \"\"\"\n",
        "  return tf.io.parse_single_example(example_proto, FEATURES_DICT)\n",
        "\n",
        "\n",
        "def to_tuple(inputs):\n",
        "  \"\"\"Function to convert a dictionary of tensors to a tuple of (inputs, outputs).\n",
        "  Turn the tensors returned by parse_tfrecord into a stack in HWC shape.\n",
        "  Args:\n",
        "    inputs: A dictionary of tensors, keyed by feature name.\n",
        "  Returns:\n",
        "    A tuple of (inputs, outputs).\n",
        "  \"\"\"\n",
        "  inputsList = [inputs.get(key) for key in FEATURES]\n",
        "  stacked = tf.stack(inputsList, axis=0)\n",
        "  # Convert from CHW to HWC\n",
        "  stacked = tf.transpose(stacked, [1, 2, 0])\n",
        "  return stacked[:,:,:len(BANDS)], stacked[:,:,len(BANDS):]\n",
        "\n",
        "\n",
        "def get_dataset(pattern):\n",
        "  \"\"\"Function to read, parse and format to tuple a set of input tfrecord files.\n",
        "  Get all the files matching the pattern, parse and convert to tuple.\n",
        "  Args:\n",
        "    pattern: A file pattern to match in a Cloud Storage bucket.\n",
        "  Returns:\n",
        "    A tf.data.Dataset\n",
        "  \"\"\"\n",
        "  glob = tf.io.gfile.glob(pattern)\n",
        "  dataset = tf.data.TFRecordDataset(glob, compression_type='GZIP')\n",
        "  dataset = dataset.map(parse_tfrecord, num_parallel_calls=5)\n",
        "  dataset = dataset.map(to_tuple, num_parallel_calls=5)\n",
        "  return dataset"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Xg1fa18336D2"
      },
      "source": [
        "Use the helpers to read in the training dataset.  Print the first record to check."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rm0qRF0fAYcC"
      },
      "source": [
        "def get_training_dataset():\n",
        "\t\"\"\"Get the preprocessed training dataset\n",
        "  Returns:\n",
        "    A tf.data.Dataset of training data.\n",
        "  \"\"\"\n",
        "\tglob = 'gs://' + BUCKET + '/' + FOLDER + '/' + TRAINING_BASE + '*'\n",
        "\tdataset = get_dataset(glob)\n",
        "\tdataset = dataset.shuffle(BUFFER_SIZE).batch(BATCH_SIZE).repeat()\n",
        "\treturn dataset\n",
        "\n",
        "training = get_training_dataset()\n",
        "\n",
        "print(iter(training.take(1)).next())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j-cQO5RL6vob"
      },
      "source": [
        "## Evaluation data\n",
        "\n",
        "Now do the same thing to get an evaluation dataset.  Note that unlike the training dataset, the evaluation dataset has a batch size of 1, is not repeated and is not shuffled."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ieKTCGiJ6xzo"
      },
      "source": [
        "def get_eval_dataset():\n",
        "\t\"\"\"Get the preprocessed evaluation dataset\n",
        "  Returns:\n",
        "    A tf.data.Dataset of evaluation data.\n",
        "  \"\"\"\n",
        "\tglob = 'gs://' + BUCKET + '/' + FOLDER + '/' + EVAL_BASE + '*'\n",
        "\tdataset = get_dataset(glob)\n",
        "\tdataset = dataset.batch(1).repeat()\n",
        "\treturn dataset\n",
        "\n",
        "evaluation = get_eval_dataset()\n",
        "\n",
        "print(iter(evaluation.take(1)).next())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9JIE7Yl87lgU"
      },
      "source": [
        "# Model\n",
        "\n",
        "Here we use the Keras implementation of the U-Net model.  The U-Net model takes 256x256 pixel patches as input and outputs per-pixel class labels output.  We can implement the model essentially unmodified.  Since roads are representnted in a binary way (0=non-road, 1=road], a ReLU activation function is suitable here."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wsnnnz56yS3l"
      },
      "source": [
        "def conv_block(input_tensor, num_filters):\n",
        "\tencoder = tf.keras.layers.Conv2D(num_filters, (3, 3), padding='same')(input_tensor)\n",
        "\tencoder = tf.keras.layers.BatchNormalization()(encoder)\n",
        "\tencoder = tf.keras.layers.Activation('relu')(encoder)\n",
        "\tencoder = tf.keras.layers.Conv2D(num_filters, (3, 3), padding='same')(encoder)\n",
        "\tencoder = tf.keras.layers.BatchNormalization()(encoder)\n",
        "\tencoder = tf.keras.layers.Activation('relu')(encoder)\n",
        "\treturn encoder\n",
        "\n",
        "def encoder_block(input_tensor, num_filters):\n",
        "\tencoder = conv_block(input_tensor, num_filters)\n",
        "\tencoder_pool = tf.keras.layers.MaxPooling2D((2, 2), strides=(2, 2))(encoder)\n",
        "\treturn encoder_pool, encoder\n",
        "\n",
        "def decoder_block(input_tensor, concat_tensor, num_filters):\n",
        "\tdecoder = tf.keras.layers.Conv2DTranspose(num_filters, (2, 2), strides=(2, 2), padding='same')(input_tensor)\n",
        "\tdecoder = tf.keras.layers.concatenate([concat_tensor, decoder], axis=-1)\n",
        "\tdecoder = tf.keras.layers.BatchNormalization()(decoder)\n",
        "\tdecoder = tf.keras.layers.Activation('relu')(decoder)\n",
        "\tdecoder = tf.keras.layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)\n",
        "\tdecoder = tf.keras.layers.BatchNormalization()(decoder)\n",
        "\tdecoder = tf.keras.layers.Activation('relu')(decoder)\n",
        "\tdecoder = tf.keras.layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)\n",
        "\tdecoder = tf.keras.layers.BatchNormalization()(decoder)\n",
        "\tdecoder = tf.keras.layers.Activation('relu')(decoder)\n",
        "\treturn decoder\n",
        "\n",
        "def get_model():\n",
        "\tinputs = tf.keras.layers.Input(shape=[None, None, len(BANDS)]) # 256\n",
        "\tencoder0_pool, encoder0 = encoder_block(inputs, 32) # 128\n",
        "\tencoder1_pool, encoder1 = encoder_block(encoder0_pool, 64) # 64\n",
        "\tencoder2_pool, encoder2 = encoder_block(encoder1_pool, 128) # 32\n",
        "\tencoder3_pool, encoder3 = encoder_block(encoder2_pool, 256) # 16\n",
        "\tencoder4_pool, encoder4 = encoder_block(encoder3_pool, 512) # 8\n",
        "\tcenter = conv_block(encoder4_pool, 1024) # center\n",
        "\tdecoder4 = decoder_block(center, encoder4, 512) # 16\n",
        "\tdecoder3 = decoder_block(decoder4, encoder3, 256) # 32\n",
        "\tdecoder2 = decoder_block(decoder3, encoder2, 128) # 64\n",
        "\tdecoder1 = decoder_block(decoder2, encoder1, 64) # 128\n",
        "\tdecoder0 = decoder_block(decoder1, encoder0, 32) # 256\n",
        "\toutputs = tf.keras.layers.Conv2D(1, (1, 1), activation='sigmoid')(decoder0)\n",
        "    # https://machinelearningmastery.com/choose-an-activation-function-for-deep-learning/\n",
        "    # https://www.quora.com/Does-it-make-sense-to-use-Relu-activation-on-the-output-neuron-for-binary-classification-If-not-why?share=1\n",
        "\n",
        "\tmodel = tf.keras.Model(inputs=[inputs], outputs=[outputs])\n",
        "\n",
        "\tmodel.compile(\n",
        "\t\toptimizer=OPTIMIZER,\n",
        "\t\tloss=LOSS,\n",
        "\t\tmetrics=METRICS)\n",
        "\n",
        "\treturn model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uu_E7OTDBCoS"
      },
      "source": [
        "# Training\n",
        "\n",
        "You train a Keras model by calling `.fit()` on it.  Here we're going to train for 20 epochs (Takes around **2-3h** with a GPU).  For production use, you probably want to optimize this parameter, for example through [hyperparamter tuning](https://cloud.google.com/ml-engine/docs/tensorflow/using-hyperparameter-tuning).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NzzaWxOhSxBy"
      },
      "source": [
        "m = get_model()\n",
        "\n",
        "# path for best model\n",
        "bpath = 'gs://' + BUCKET + '/' + FOLDER + '/' + 'best_model-epoch{epoch:02d}-loss{loss:.3f}-dice_coef{dice_coef:.3f}'\n",
        "\n",
        "# best model is saved in every checkpoint\n",
        "model_checkpoint_callback = tf.keras.callbacks.ModelCheckpoint(\n",
        "    filepath=bpath,\n",
        "    save_weights_only=True,\n",
        "    monitor='val_loss',\n",
        "    mode='min',\n",
        "    save_best_only=True)\n",
        "\n",
        "\n",
        "# history is a dictionary holding loss & accuracy at each epoch\n",
        "history = m.fit(\n",
        "    x=training,\n",
        "    epochs=EPOCHS,\n",
        "    steps_per_epoch=int(TRAIN_SIZE / BATCH_SIZE),\n",
        "    validation_data=evaluation,\n",
        "    validation_steps=EVAL_SIZE,\n",
        "    # for verbose see https://github.com/tensorflow/tensorflow/issues/37876\n",
        "    verbose=1,\n",
        "    callbacks=[model_checkpoint_callback]\n",
        ")\n",
        "\n",
        "# For training loss, keras does a running average over the batches. For validation loss, a conventional average over\n",
        "# all the batches in validation data is performed. The training accuracy is the average of the accuracy values for each\n",
        "# batch of training data during training. (https://github.com/keras-team/keras/issues/10426)\n",
        "\n",
        "# save history to file\n",
        "hpath = 'gs://' + BUCKET + '/' + FOLDER + '/' + 'MODEL_HISTORY.npy'\n",
        "np.save(hpath,history.history)\n",
        "\n",
        "# get the latest (best) model saved\n",
        "path = 'gs://' + BUCKET + '/' + FOLDER\n",
        "\n",
        "# path of best model\n",
        "latest = tf.train.latest_checkpoint(path)\n",
        "m.load_weights(latest)\n",
        "\n",
        "m.summary()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "718KsSBmybar"
      },
      "source": [
        "## Load a model and continue training"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kbdvyq-hyhV2"
      },
      "source": [
        "EPOCHS_LEFT = 29\n",
        "\n",
        "# load the checkpoint from disk\n",
        "print(\"loading model...\")\n",
        "\n",
        "m = get_model()\n",
        "\n",
        "# get the latest (best) model saved\n",
        "path = 'gs://' + BUCKET + '/' + FOLDER\n",
        "\n",
        "# path of best model\n",
        "latest = tf.train.latest_checkpoint(path)\n",
        "\n",
        "m.load_weights(latest)\n",
        "\n",
        "# loss,acc = m.evaluate(evaluation, verbose=2)\n",
        "# print(\"Restored model, accuracy: {:5.2f}%\".format(100*acc))\n",
        "\n",
        "# path for best model\n",
        "bpath = 'gs://' + BUCKET + '/' + FOLDER + '/' + 'best_model-epoch{epoch:02d}-loss{loss:.3f}-dice_coef{dice_coef:.3f}'\n",
        "\n",
        "# best model is saved in every checkpoint\n",
        "model_checkpoint_callback = tf.keras.callbacks.ModelCheckpoint(\n",
        "    filepath=bpath,\n",
        "    save_weights_only=True,\n",
        "    monitor='dice_coef',\n",
        "    mode='max',\n",
        "    save_best_only=True)\n",
        "\n",
        "# resume the training where we left off\n",
        "history2 = m.fit(\n",
        "  x=training,\n",
        "  epochs=EPOCHS_LEFT,\n",
        "  steps_per_epoch=int(TRAIN_SIZE / BATCH_SIZE),\n",
        "  validation_data=evaluation,\n",
        "  validation_steps=EVAL_SIZE,\n",
        "  # for verbose see https://github.com/tensorflow/tensorflow/issues/37876\n",
        "  verbose=1,\n",
        "  callbacks=[model_checkpoint_callback]\n",
        ")\n",
        "\n",
        "\n",
        "# # save history to file\n",
        "# hpath = 'gs://' + BUCKET + '/' + FOLDER + '/' + 'MODEL_HISTORY.npy'\n",
        "# np.save(hpath,history2.history)\n",
        "\n",
        "# get the latest (best) model saved\n",
        "path = 'gs://' + BUCKET + '/' + FOLDER\n",
        "\n",
        "# path of best model\n",
        "latest = tf.train.latest_checkpoint(path)\n",
        "m.load_weights(latest)\n",
        "\n",
        "m.summary()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gmV3zeT1ysis"
      },
      "source": [
        "## Accuracy\n",
        "plot the train and validation loss as well as the train and validation accuracy"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zQbpn4tSywWq"
      },
      "source": [
        "# # load history from file\n",
        "# hpath = 'gs://' + BUCKET + '/' + FOLDER + '/' + 'MODEL_HISTORY.npy'\n",
        "# history=np.load(hpath,allow_pickle='TRUE').item()\n",
        "\n",
        "# Plotting both losses simultaneously\n",
        "l = plt.plot(history.history['loss'], color='gray', label='dice loss')\n",
        "\n",
        "l = plt.title('model loss')\n",
        "l = plt.xlabel('epoch')\n",
        "l = plt.legend(loc='upper left')\n",
        "\n",
        "# save to file\n",
        "plt.savefig(os.path.join('gs://' + BUCKET + '/' + FOLDER + '/' +  'loss.jpg'), dpi=300)\n",
        "\n",
        "l = plt.show()\n",
        "\n",
        "# Plotting both accuracy simultaneously\n",
        "a = plt.plot(history.history['recall'], color='blue', label='recall')\n",
        "a = plt.plot(history.history['binary_accuracy'], color='orange', label='binary_accuracy')\n",
        "a = plt.plot(history.history['mean_io_u'], color='red', label='IoU')\n",
        "a = plt.plot(history.history['dice_coef'], color='gray', label='dice coeficient')\n",
        "\n",
        "a = plt.title('model accuracy')\n",
        "a = plt.ylabel('accuracy')\n",
        "a = plt.xlabel('epoch')\n",
        "a = plt.legend(loc='lower right')\n",
        "\n",
        "# save to file\n",
        "plt.savefig(os.path.join('gs://' + BUCKET + '/' + FOLDER + '/' + 'accuracy.jpg'), dpi=300)\n",
        "\n",
        "a = plt.show()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HnlYZTZRy71K"
      },
      "source": [
        "The following code loads the best epoch for the model trained above, which you\n",
        "can use for predictions right away. Loading model takes around **15-40m**.\n",
        "\n",
        "Model at **Epoch 461**\n",
        "\n",
        "Accuracy:\n",
        "`loss: -0.5551 - binary_accuracy: 0.9444 - dice_coef: 0.5551 - mean_io_u: 0.4755 - recall: 0.5634`\n",
        "`Total params: 31,127,361 - Trainable params: 31,111,361 - Non-trainable params: 16,000`\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-RJpNfEUS1qp"
      },
      "source": [
        "m = get_model()\n",
        "\n",
        "# get the latest (best) model saved\n",
        "path = 'gs://' + BUCKET + '/' + FOLDER\n",
        "\n",
        "# path of best model\n",
        "latest = tf.train.latest_checkpoint(path)\n",
        "m.load_weights(latest)\n",
        "\n",
        "m.summary()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J1ySNup0xCqN"
      },
      "source": [
        "# Prediction\n",
        "\n",
        "The prediction pipeline is:\n",
        "\n",
        "1.  Export imagery on which to do predictions from Earth Engine in TFRecord format to a Cloud Storge bucket.\n",
        "2.  Use the trained model to make the predictions.\n",
        "3.  Write the predictions to a TFRecord file in a Cloud Storage.\n",
        "4.  Upload the predictions TFRecord file to Earth Engine.\n",
        "\n",
        "The following functions handle this process.  It's useful to separate the export from the predictions so that you can experiment with different models without running the export every time."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M3WDAa-RUpXP"
      },
      "source": [
        "def doExport(out_image_base, kernel_buffer, region, tileNum):\n",
        "  \"\"\"Run the image export task.  Block until complete.\n",
        "  \"\"\"\n",
        "  task = ee.batch.Export.image.toCloudStorage(\n",
        "    image = image_res.select(BANDS),\n",
        "    description = out_image_base + str(tileNum) + '_',\n",
        "    bucket = BUCKET,\n",
        "    fileNamePrefix = FOLDER + '/' + out_image_base + str(tileNum) + '_',\n",
        "    region = region.getInfo()['coordinates'],\n",
        "    scale = 1,\n",
        "    fileFormat = 'TFRecord',\n",
        "    maxPixels = 1e10,\n",
        "    formatOptions = {\n",
        "      'patchDimensions': KERNEL_SHAPE,\n",
        "      'kernelSize': kernel_buffer,\n",
        "      'compressed': True,\n",
        "      'maxFileSize': 104857600\n",
        "    }\n",
        "  )\n",
        "  task.start()\n",
        "\n",
        "  # # Block until the task completes.\n",
        "  # print('Running image export to Cloud Storage...')\n",
        "  # import time\n",
        "  # while task.active():\n",
        "  #   time.sleep(30)\n",
        "\n",
        "  # # Error condition\n",
        "  # if task.status()['state'] != 'COMPLETED':\n",
        "  #   print('Error with image export.')\n",
        "  # else:\n",
        "  #   print('Image export completed.')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zb_9_FflygVw"
      },
      "source": [
        "def doPrediction(out_image_base, user_folder, kernel_buffer, region):\n",
        "  \"\"\"Perform inference on exported imagery, upload to Earth Engine.\n",
        "  \"\"\"\n",
        "\n",
        "  print('Looking for TFRecord files...')\n",
        "\n",
        "  # Get a list of all the files in the output bucket.\n",
        "  filesList = !gsutil ls 'gs://'{BUCKET}'/'{FOLDER}\n",
        "\n",
        "  # Get only the files generated by the image export for the specific tile.\n",
        "  exportFilesList = [s for s in filesList if out_image_base in s]\n",
        "\n",
        "  # Get the list of image files and the JSON mixer file.\n",
        "  imageFilesList = []\n",
        "  jsonFile = None\n",
        "  for f in exportFilesList:\n",
        "    if f.endswith('.tfrecord.gz'):\n",
        "      imageFilesList.append(f)\n",
        "    elif f.endswith('.json'):\n",
        "      jsonFile = f\n",
        "\n",
        "  # Make sure the files are in the right order.\n",
        "  imageFilesList.sort()\n",
        "\n",
        "  from pprint import pprint\n",
        "  pprint(imageFilesList)\n",
        "  print(jsonFile)\n",
        "\n",
        "  import json\n",
        "  # Load the contents of the mixer file to a JSON object.\n",
        "  jsonText = !gsutil cat {jsonFile}\n",
        "  # Get a single string w/ newlines from the IPython.utils.text.SList\n",
        "  mixer = json.loads(jsonText.nlstr)\n",
        "  pprint(mixer)\n",
        "  patches = mixer['totalPatches']\n",
        "\n",
        "  # Get set up for prediction.\n",
        "  x_buffer = int(kernel_buffer[0] / 2)\n",
        "  y_buffer = int(kernel_buffer[1] / 2)\n",
        "\n",
        "  buffered_shape = [\n",
        "      KERNEL_SHAPE[0] + kernel_buffer[0],\n",
        "      KERNEL_SHAPE[1] + kernel_buffer[1]]\n",
        "\n",
        "  imageColumns = [\n",
        "    tf.io.FixedLenFeature(shape=buffered_shape, dtype=tf.float32)\n",
        "      for k in BANDS\n",
        "  ]\n",
        "\n",
        "  imageFeaturesDict = dict(zip(BANDS, imageColumns))\n",
        "\n",
        "  def parse_image(example_proto):\n",
        "    return tf.io.parse_single_example(example_proto, imageFeaturesDict)\n",
        "\n",
        "  def toTupleImage(inputs):\n",
        "    inputsList = [inputs.get(key) for key in BANDS]\n",
        "    stacked = tf.stack(inputsList, axis=0)\n",
        "    stacked = tf.transpose(stacked, [1, 2, 0])\n",
        "    return stacked\n",
        "\n",
        "   # Create a dataset from the TFRecord file(s) in Cloud Storage.\n",
        "  imageDataset = tf.data.TFRecordDataset(imageFilesList, compression_type='GZIP')\n",
        "  imageDataset = imageDataset.map(parse_image, num_parallel_calls=5)\n",
        "  imageDataset = imageDataset.map(toTupleImage).batch(1)\n",
        "\n",
        "  # Perform inference.\n",
        "  print('Running predictions...')\n",
        "  predictions = m.predict(imageDataset, steps=patches, verbose=1)\n",
        "  # print(predictions[0])\n",
        "\n",
        "  print('Writing predictions...')\n",
        "  out_image_file = 'gs://' + BUCKET + '/' + FOLDER + '/' + out_image_base + '.TFRecord'\n",
        "  writer = tf.io.TFRecordWriter(out_image_file)\n",
        "  patches = 0\n",
        "  for predictionPatch in predictions:\n",
        "    print('Writing patch ' + str(patches) + '...')\n",
        "    predictionPatch = predictionPatch[\n",
        "        x_buffer:x_buffer+KERNEL_SIZE, y_buffer:y_buffer+KERNEL_SIZE]\n",
        "\n",
        "    # Create an example.\n",
        "    example = tf.train.Example(\n",
        "      features=tf.train.Features(\n",
        "        feature={\n",
        "          'r': tf.train.Feature(\n",
        "              float_list=tf.train.FloatList(\n",
        "                  value=predictionPatch.flatten()))\n",
        "        }\n",
        "      )\n",
        "    )\n",
        "    # Write the example.\n",
        "    writer.write(example.SerializeToString())\n",
        "    patches += 1\n",
        "\n",
        "  writer.close()\n",
        "\n",
        "  # Start the upload.\n",
        "  out_image_asset = user_folder + '/' + out_image_base\n",
        "  !earthengine upload image --asset_id={out_image_asset} {out_image_file} {jsonFile}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LZqlymOehnQO"
      },
      "source": [
        "Now there's all the code needed to run the prediction pipeline, all that remains is to specify the output region in which to do the prediction, the names of the output files, where to put them, and the shape of the outputs.  In terms of the shape, the model is trained on 256x256 patches, but can work (in theory) on any patch that's big enough with even dimensions ([reference](https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Long_Fully_Convolutional_Networks_2015_CVPR_paper.pdf)).  Because of tile boundary artifacts, give the model slightly larger patches for prediction, then clip out the middle 256x256 patch.  This is controlled with a kernel buffer, half the size of which will extend beyond the kernel buffer.  For example, specifying a 128x128 kernel will append 64 pixels on each side of the patch, to ensure that the pixels in the output are taken from inputs completely covered by the kernel."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FPANwc7B1-TS"
      },
      "source": [
        "# Output assets folder: YOUR FOLDER\n",
        "user_folder = 'users/n-verde' # INSERT YOUR FOLDER HERE.\n",
        "\n",
        "# Base file name to use for TFRecord files and assets.\n",
        "attica_image_base = '1171_LAS_'\n",
        "# Half this will extend on the sides of each patch.\n",
        "attica_kernel_buffer = [128, 128]\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pMomKb2sPlHD"
      },
      "source": [
        "Run the export (takes around **8m** for each tile, so approx. **8.5h** for 64 tiles)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lLNEOLkXWvSi"
      },
      "source": [
        "# Run the export for tiles of Attica\n",
        "\n",
        "# Attica in tiles\n",
        "\n",
        "tilesList = [36] # if you want to run for specific tiles\n",
        "\n",
        "# first tile\n",
        "tile_number = int(1)\n",
        "\n",
        "# loop\n",
        "# for i in range(1,(len(splitted.getInfo())+1)):\n",
        "for i in tilesList: # if you want to run for specific tiles\n",
        "\n",
        "  # # option a. without buffer\n",
        "  # tile = splitted.getInfo()[i-1]\n",
        "  # region = tile[\"coordinates\"]\n",
        "  # # print(region)\n",
        "\n",
        "  # option b. with buffer\n",
        "  # buffer the tile 250m to avoid data gaps between tiles\n",
        "  # the distance of 250m was found empirically by running the code and seeing the gaps\n",
        "  prj = ee.Projection('EPSG:3035');  # European projection\n",
        "  buffer = ee.Geometry(ee.List(splitted.get(i-1))).buffer(250,None,prj).bounds(0.1)\n",
        "  buffered_tile = buffer.getInfo()\n",
        "  buffered_region = buffered_tile[\"coordinates\"]\n",
        "  region = buffered_region\n",
        "  # print(region)\n",
        "\n",
        "  tile_region = ee.Geometry.Polygon(region, None, False)\n",
        "\n",
        "  # convert \"1\" to \"01\"\n",
        "  tile_number_string = str(tile_number).zfill(2)\n",
        "\n",
        "  doExport(attica_image_base, attica_kernel_buffer, tile_region, tile_number_string)\n",
        "\n",
        "  print('done exporting for tile: ', tile_number_string, region)\n",
        "  tile_number = tile_number + 1\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BeK6vPauPw4g"
      },
      "source": [
        "Run the prediction and upload to GEE (for whole attica (64 tiles), *prediction* takes around **10s-2m** per tile (with GPU) so **1h** in total (with GPU) or **10m** per tile (without GPU) so **10.5h** total (without GPU), and *upload* to GEE takes around **2h-3h** for each tile. For whole Attica in tiles takes around **10h** in total (GEE can run 20 ingestions simultaneously)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KxACnxKFrQ_J"
      },
      "source": [
        "# Run the prediction for tiles of Attica\n",
        "\n",
        "# just to invoke earth engine so it doesn't disconnect\n",
        "ee.Authenticate()\n",
        "ee.Initialize()\n",
        "\n",
        "tilesList = [36] # if you want to run for specific tiles\n",
        "\n",
        "# first tile\n",
        "tile_number = int(1)\n",
        "\n",
        "# loop\n",
        "# for i in range(1,(len(splitted.getInfo())+1)):\n",
        "for i in tilesList: # if you want to run for specific tiles\n",
        "\n",
        "  # # option a. without buffer\n",
        "  # tile = splitted.getInfo()[i-1]\n",
        "  # region = tile[\"coordinates\"]\n",
        "  # # print(region)\n",
        "\n",
        "  # option b. with buffer\n",
        "  # buffer the tile 250m to avoid data gaps between tiles\n",
        "  # the distance of 250m was found empirically by running the code and seeing the gaps\n",
        "  prj = ee.Projection('EPSG:3035');  # European projection\n",
        "  buffer = ee.Geometry(ee.List(splitted.get(i-1))).buffer(250,None,prj).bounds(0.1)\n",
        "  buffered_tile = buffer.getInfo()\n",
        "  buffered_region = buffered_tile[\"coordinates\"]\n",
        "  region = buffered_region\n",
        "  # print(region)\n",
        "\n",
        "  tile_region = ee.Geometry.Polygon(region, None, False)\n",
        "\n",
        "  # convert \"1\" to \"01\"\n",
        "  tile_number_string = str(tile_number).zfill(2)\n",
        "\n",
        "  attica_image_base2 = attica_image_base + str(tile_number_string)\n",
        "\n",
        "  doPrediction(attica_image_base2, user_folder, attica_kernel_buffer, tile_region)\n",
        "\n",
        "  print('done exporting for tile: ', tile_number_string, region)\n",
        "  tile_number = tile_number + 1\n"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}