1171_paper_UOS_SAR.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.10"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_SHAc5qbiR8l"
      },
      "source": [
        "# Introduction\n",
        "\n",
        "Classify Urban Open Spaces from planetscope 3m imagery, with FCNN U-Net. Based on GEE FCNN notebook."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ir8E4Cghtnl_"
      },
      "source": [
        "##### Version\n",
        "\n",
        "v13\n",
        "For test site."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_MJ4kW1pEhwP"
      },
      "source": [
        "# Libraries & imports\n",
        "\n",
        "Authenticate and import packages."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6XrAU3TZ6MJc"
      },
      "source": [
        "# Cloud authentication.\n",
        "from google.colab import auth\n",
        "auth.authenticate_user()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jat01FEoUMqg",
        "scrolled": true
      },
      "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": "59PzK5sd50GN"
      },
      "source": [
        "## Specify your Cloud Storage Bucket"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OqNpnmu050GN"
      },
      "source": [
        "# INSERT YOUR BUCKET HERE:\n",
        "BUCKET = 'n-verde_bucket'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wmfKLl9XcnGJ"
      },
      "source": [
        "## Set global variables"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "psz7wJKalaoj"
      },
      "source": [
        "# Specify names locations for outputs.\n",
        "FOLDER = '1171-UOS-256_2500-SAR'\n",
        "TRAINING_BASE = 'training_patches'\n",
        "EVAL_BASE = 'eval_patches'\n",
        "\n",
        "# ---------------------- FEATURES --------------------------\n",
        "# Specify inputs (Planetscope & S1 bands) to the model and the response variable.\n",
        "opticalBands = ['b1', 'b2', 'b3', 'b4', 'NDVI']\n",
        "sarBands = []\n",
        "otherData = []  # predicted roads\n",
        "BANDS = opticalBands + sarBands + otherData\n",
        "RESPONSE = 'ugs'\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 = 2500   # 1500  # 2500  # 1200  # 16000\n",
        "EVAL_SIZE = 1071    # 643   # 1071  # 320    # 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 = 16 # 4\n",
        "# https://machinelearningmastery.com/how-to-stop-training-deep-neural-networks-at-the-right-time-using-early-stopping/\n",
        "EPOCHS = 500\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 = 2000\n",
        "\n",
        "# ---------------------- OPTIMIZER --------------------------\n",
        "\n",
        "# Optimizer with POLYNOMIAL learning rate decay\n",
        "# https://keras.io/api/optimizers/learning_rate_schedules/polynomial_decay/\n",
        "lr_schedule = tf.keras.optimizers.schedules.PolynomialDecay(\n",
        "    initial_learning_rate=0.01,   # 0.01\n",
        "    decay_steps=20,    # 200\n",
        "    end_learning_rate=0.0001,    # 0.00001\n",
        "    power=0.5,\n",
        "    cycle=False,\n",
        ")\n",
        "OPTIMIZER = tf.keras.optimizers.SGD(learning_rate=lr_schedule, momentum=0.9)\n",
        "\n",
        "# ---------------------- LOSS --------------------------\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",
        "## All\n",
        "\n",
        "Gather and setup the imagery to use for inputs (predictors) (Planetscope mosaic, indices, SAR etc).  Display it in the notebook for a sanity check.\n",
        "\n",
        "Prepare the response (what we want to predict). This is the Urban Open Spaces from OSM, rasterized. Display to check."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-IlgXu-vcUEY",
        "scrolled": false
      },
      "source": [
        "# Use folium to visualize the imagery.\n",
        "map = folium.Map(location=[37.981892554434936, 23.7269115882649])\n",
        "\n",
        "\n",
        "# ------------------------------------------\n",
        "# planetscope mosiac image of attica\n",
        "planetscope = ee.Image('users/n-verde/PhD_1171/attica_feather_mask_buff400')\n",
        "\n",
        "mapid = planetscope.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",
        "# ------------------------------------------\n",
        "# ugs ground truth\n",
        "gt = ee.Image('users/n-verde/PhD_1171/osm_open_areas_EDIT_v2_3m')\n",
        "gt = gt.rename('ugs')\n",
        "\n",
        "mapid = gt.getMapId({'bands': ['ugs'], '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",
        "# ------------------------------------------\n",
        "# S1 imagery\n",
        "s1_col = ee.ImageCollection(\"COPERNICUS/S1_GRD_FLOAT\")\n",
        "\n",
        "# get the s1 image closest to the date of the planetscope image\n",
        "# get only descending orbit\n",
        "s1_col = s1_col.filterDate('2020-08-27', '2020-09-02')\n",
        "s1_col = s1_col.filterMetadata('instrumentMode', 'equals', 'IW')\n",
        "s1_col = s1_col.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))\n",
        "s1_col = s1_col.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))\n",
        "s1_col = s1_col.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))\n",
        "s1_im = ee.Image(s1_col.mean()).mask(planetscope.select('b1')).select('VV', 'VH')\n",
        "\n",
        "mapid = s1_im.getMapId({'bands': ['VV'], 'min': 0, 'max': 3.5})\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='s1_im',\n",
        "  ).add_to(map)\n",
        "\n",
        "map.add_child(folium.LayerControl())\n",
        "\n",
        "map\n",
        "\n",
        "## SAR Processing"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ji-fPYGs3mMz"
      },
      "source": [
        "### Transformations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fQTv5ECV3tAZ"
      },
      "source": [
        "############################################################\n",
        "# This function calculates the VVH indicator.\n",
        "# Used for better capturing building height from Sentinel-1 at 500m spatial resolution\n",
        "# (Li et al., 2020)\n",
        "\n",
        "def addVVH(image):\n",
        "  vvh = None\n",
        "  vvh = image.expression('VV * (gamma ** VH)',{\n",
        "          'VV': image.select('VV'),\n",
        "          'VH': image.select('VH'),\n",
        "          'gamma': ee.Number(5)\n",
        "        })\n",
        "  return image.addBands(vvh.rename('VVH'))\n",
        "\n",
        "s1_im = addVVH(s1_im)\n",
        "\n",
        "############################################################\n",
        "# This function calculates building height.\n",
        "# Used for Sentinel-1 at 500m spatial resolution\n",
        "# (Li et al., 2020)\n",
        "\n",
        "def buildingHeightS1(image):\n",
        "  h = None\n",
        "  h = image.expression('e ** (a * (VVH ** b) + c)',{\n",
        "          'VVH': image.select('VVH'),\n",
        "          'e': ee.Number(2.71828),\n",
        "          'a': ee.Number(-23.61),\n",
        "          'b': ee.Number(-0.06),\n",
        "          'c': ee.Number(26.10),\n",
        "        })\n",
        "  return image.addBands(h.rename('buildingH_Natural'))\n",
        "\n",
        "s1_im = buildingHeightS1(s1_im)\n",
        "\n",
        "print(\"Calculated building height index\")\n",
        "\n",
        "# print(s1_im.getInfo())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fiKvWmEm2APM"
      },
      "source": [
        "### Turn to dB\n",
        "Turn SAR bands from natural values to dB."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PEDMrpJx2KnB"
      },
      "source": [
        "############################################################\n",
        "# This function converts to dB from Natural\n",
        "def toDB(img):\n",
        "    return ee.Image(img).log10().multiply(10.0)\n",
        "\n",
        "############################################################\n",
        "# This function cinverts from dB to Natural\n",
        "def toNatural(img):\n",
        "    return ee.Image(10.0).pow(img.divide(10.0));\n",
        "\n",
        "# this image is used for building heights\n",
        "buildingH_db = toDB(s1_im.select('buildingH_Natural'))\n",
        "s1_im_dB = s1_im.addBands(buildingH_db.rename('buildingH'))\n",
        "\n",
        "print(\"Converted SAR imagery from natural to dB\")\n",
        "# print(s1_im_dB.getInfo())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pA-ljY77bihF"
      },
      "source": [
        "### Normalization\n",
        "\n",
        "Find the min and max of the S1 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",
        "scrolled": true
      },
      "source": [
        "# define the study area (area of planetscope image)\n",
        "aoi = ee.FeatureCollection('users/n-verde/PhD_1171/Athens_buff400')\n",
        "\n",
        "print(\"finding min and max of SAR S1 bands, in order to normalize them...\")\n",
        "print(\"using percentiles instead of absolute min and max...\")\n",
        "print(\"----------\")\n",
        "\n",
        "def findMinSAR(imgBand, bandName, min_prec):\n",
        "    MIN = ee.Number(imgBand.reduceRegion(**{\n",
        "    'reducer': ee.Reducer.percentile([min_prec]),\n",
        "    'geometry': aoi,\n",
        "    'scale': 10,\n",
        "    'tileScale': 4,\n",
        "    # 'bestEffort': true,\n",
        "    'maxPixels': 1e9\n",
        "    }).get(bandName))\n",
        "\n",
        "    return MIN\n",
        "\n",
        "def findMaxSAR(imgBand, bandName, max_perc):\n",
        "    MAX = ee.Number(imgBand.reduceRegion(**{\n",
        "    'reducer': ee.Reducer.percentile([max_perc]),\n",
        "    'geometry': aoi,\n",
        "    'scale': 10,\n",
        "    'tileScale': 4,\n",
        "    # 'bestEffort': true,\n",
        "    'maxPixels': 1e9\n",
        "    }).get(bandName))\n",
        "\n",
        "    return MAX\n",
        "\n",
        "# VV ----------\n",
        "\n",
        "bandName = \"VV\"\n",
        "imgBand = s1_im_dB.select(bandName)\n",
        "\n",
        "# 15.85 to 84.15 is stretch to 1σ\n",
        "# 2.3 to 97.7 is stretch to 2σ\n",
        "# 0.15 to 99.85 is stretch to 3σ\n",
        "\n",
        "# MIN = ee.Number(findMinSAR(imgBand, bandName, 0.15))\n",
        "MIN = ee.Number(0.35441243231694636)\n",
        "\n",
        "# MAX = ee.Number(findMaxSAR(imgBand, bandName, 99.85))\n",
        "MAX = ee.Number(0.3544124323169465)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\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",
        "    MIN = ee.Number(0) # change min to 0\n",
        "\n",
        "min_VV = MIN\n",
        "max_VV = MAX\n",
        "\n",
        "print(\"min of VV:\", min_VV.getInfo())\n",
        "print(\"max of VV:\", max_VV.getInfo())\n",
        "\n",
        "# VH ----------\n",
        "\n",
        "bandName = \"VH\"\n",
        "imgBand = s1_im_dB.select(bandName)\n",
        "\n",
        "# MIN = findMinSAR(imgBand, bandName, 0.15)\n",
        "MIN = ee.Number(0.05221765024397127)\n",
        "\n",
        "# MAX = findMaxSAR(imgBand, bandName, 99.85)\n",
        "MAX = ee.Number(1.3602214673473254)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\n",
        "\n",
        "min_VH = MIN\n",
        "max_VH = MAX\n",
        "\n",
        "print(\"min of VH:\", min_VH.getInfo())\n",
        "print(\"max of VH:\", max_VH.getInfo())\n",
        "\n",
        "# buildingH ----------\n",
        "\n",
        "bandName = \"buildingH\"\n",
        "imgBand = s1_im_dB.select(bandName)\n",
        "\n",
        "# MIN = findMinSAR(imgBand, bandName, 0.15)\n",
        "MIN = ee.Number(-34.476049329083615)\n",
        "\n",
        "# MAX = findMaxSAR(imgBand, bandName, 99.85)\n",
        "MAX = ee.Number(31.480003022441355)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\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",
        "    MIN = ee.Number(0) # change min to 0\n",
        "\n",
        "min_buildingH = MIN\n",
        "max_buildingH = MAX\n",
        "\n",
        "print(\"min of buildingH:\", min_buildingH.getInfo())\n",
        "print(\"max of buildingH:\", max_buildingH.getInfo())\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 ----------\n",
        "\n",
        "# VV ----------\n",
        "VV_res = rescaleFixedMinMax(s1_im_dB, \"VV\", min_VV, max_VV)\n",
        "\n",
        "# VH ----------\n",
        "VH_res = rescaleFixedMinMax(s1_im_dB, \"VH\", min_VH, max_VH)\n",
        "\n",
        "# buildingH -----------\n",
        "buildingH_res = rescaleFixedMinMax(s1_im_dB,\"buildingH\", min_buildingH, max_buildingH)\n",
        "\n",
        "print(\"done rescaling SAR imagery to 0-1\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kj4qpQW0LH8D"
      },
      "source": [
        "## Optical Processing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z-AYUs37LPcl"
      },
      "source": [
        "### Transformations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FBE5HhIqRfdZ"
      },
      "source": [
        "#############################################################\n",
        "# This function calculates NDVI from a S2 image, and adds it to the image stack.\n",
        "def addNDVI(image):\n",
        "    return image.addBands(image.normalizedDifference(['b4', 'b3']).rename('NDVI')) # change new band name\n",
        "\n",
        "planetscope = addNDVI(planetscope)\n",
        "\n",
        "print(\"Calculated NDVI\")\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g_KQLZewLAd8"
      },
      "source": [
        "### Normalization"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AlLWk-tSLAv9"
      },
      "source": [
        "# define the study area (area of planetscope image)\n",
        "aoi = ee.FeatureCollection('users/n-verde/PhD_1171/Athens_buff400')\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 = planetscope.select(bandName)\n",
        "\n",
        "# MIN = findMin(imgBand, bandName)\n",
        "MIN = ee.Number(50)\n",
        "\n",
        "# MAX = findMax(imgBand, bandName)\n",
        "MAX = ee.Number(8233)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\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 = planetscope.select(bandName)\n",
        "\n",
        "# MIN = findMin(imgBand, bandName)\n",
        "MIN = ee.Number(251)\n",
        "\n",
        "# MAX = findMax(imgBand, bandName)\n",
        "MAX = ee.Number(9164)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\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 = planetscope.select(bandName)\n",
        "\n",
        "# MIN = findMin(imgBand, bandName)\n",
        "MIN = ee.Number(159)\n",
        "\n",
        "# MAX = findMax(imgBand, bandName)\n",
        "MAX = ee.Number(11765)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\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 = planetscope.select(bandName)\n",
        "\n",
        "# MIN = findMin(imgBand, bandName)\n",
        "MIN = ee.Number(1)\n",
        "\n",
        "# MAX = findMax(imgBand, bandName)\n",
        "MAX = ee.Number(12127)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\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",
        "# NDVI ----------\n",
        "\n",
        "bandName = \"NDVI\"\n",
        "imgBand = planetscope.select(bandName)\n",
        "\n",
        "# MIN = findMin(imgBand, bandName)\n",
        "MIN = ee.Number(-0.9998300186979432)\n",
        "\n",
        "# MAX = findMax(imgBand, bandName)\n",
        "MAX = ee.Number(0.8263786957433424)\n",
        "\n",
        "print('old min', MIN.getInfo())\n",
        "print('old max', MAX.getInfo())\n",
        "\n",
        "test = ee.Algorithms.If(MIN.lt(ee.Number(0)), ee.Number(1), ee.Number(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",
        "    MIN = ee.Number(0) # change min to 0\n",
        "\n",
        "min_NDVI = MIN\n",
        "max_NDVI = MAX\n",
        "\n",
        "print(\"min of NDVI:\", min_NDVI.getInfo())\n",
        "print(\"max of NDVI:\", max_NDVI.getInfo())\n",
        "\n",
        "print(\"----------\")\n",
        "print(\"done!\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yYafKXtxLdf6"
      },
      "source": [
        "Rescale to [0-1] (normalize)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9-pOkPvdLmcj"
      },
      "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(planetscope,\"b1\", min_b1, max_b1)\n",
        "\n",
        "# b2 ----------\n",
        "b2_res = rescaleFixedMinMax(planetscope,\"b2\", min_b2, max_b2)\n",
        "\n",
        "# b3 ----------\n",
        "b3_res = rescaleFixedMinMax(planetscope,\"b3\", min_b3, max_b3)\n",
        "\n",
        "# b4 ----------\n",
        "b4_res = rescaleFixedMinMax(planetscope,\"b4\", min_b4, max_b4)\n",
        "\n",
        "# NDVI ----------\n",
        "NDVI_res = rescaleFixedMinMax(planetscope,\"NDVI\", min_NDVI, max_NDVI)\n",
        "\n",
        "print(\"done rescaling optical imagery to 0-1\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GFfMcKiPV82n"
      },
      "source": [
        "## Stack"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ndsBPWFNV_hE"
      },
      "source": [
        "# merge the bands to a single image\n",
        "image_res = b1_res.addBands(b2_res).addBands(b3_res).addBands(b4_res).addBands(NDVI_res)\n",
        "image_res = image_res.addBands(VV_res).addBands(VH_res).addBands(buildingH_res)\n",
        "\n",
        "# display rescaled image\n",
        "mapid = image_res.getMapId({'bands': ['NDVI'], 'min': 0, 'max': 1})\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='rescaled',\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/BOXES_512x512_TRAIN_VAL')\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 image with response\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",
        "print('features for model: ', featureStack.bandNames().getInfo())\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": "RMXdPqLS50GX"
      },
      "source": [
        "# Split image to tiles\n",
        "Define a function for displaying Earth Engine image tiles on a folium map"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "huT1e5UT50GX"
      },
      "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": "vZRDIZ7h50GY"
      },
      "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 = planetscope.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": "dLiCDJoZ50GY"
      },
      "source": [
        "# show on Map\n",
        "\n",
        "no =29\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 = planetscope\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/UOS_grid_512_TRAINING')\n",
        "evalPolys = ee.FeatureCollection('users/n-verde/PhD_1171/UOS_grid_512_VALIDATION')\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": [
        "Take samples 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, 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 **1h** in GEE."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "cellView": "both",
        "id": "FyRpvwENxE-A"
      },
      "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 to Cloud Storage 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",
        "scrolled": true
      },
      "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",
        "    # https://stackoverflow.com/questions/53514495/what-does-batch-repeat-and-shuffle-do-with-tensorflow-dataset\n",
        "    # https://www.tensorflow.org/api_docs/python/tf/data/experimental/shuffle_and_repeat\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",
        "scrolled": true
      },
      "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(training.take(1)).next())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "abarD_ue50Gb"
      },
      "source": [
        "## Turn tfrecord samples to images\n",
        "If you want to see how the samples look like, turn them from tfrecords to images (Only works on Jupyter when run localy)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "id": "rSPoe7Hs50Gb"
      },
      "source": [
        "# import glob\n",
        "# import numpy as np\n",
        "\n",
        "# input_path = 'gs://' + BUCKET + '/' + FOLDER + '/' + TRAINING_BASE + '/'\n",
        "# output_path = 'gs://' + BUCKET + '/' + FOLDER + '/' + TRAINING_BASE + '/' + 'images'\n",
        "\n",
        "\n",
        "# files = glob.glob(input_path + '/*.tfrecord.gz')\n",
        "# if len(files) == 0:\n",
        "#     print('invalid input path : ' + input_path)\n",
        "#     quit(-1)\n",
        "# print('input files  : {}'.format(input_path + '/*.tfrecord.gz'))\n",
        "\n",
        "# if not os.path.exists(output_path):\n",
        "#     os.mkdir(output_path)\n",
        "#     print('output dir : {} (auto created)'.format(output_path))\n",
        "# else:\n",
        "#     print('output dir : {} '.format(output_path))\n",
        "\n",
        "# num = 0\n",
        "\n",
        "# for file in files:\n",
        "\n",
        "#     dataset = tf.data.TFRecordDataset(file, compression_type='GZIP')\n",
        "\n",
        "#     parsedDataset = dataset.map(lambda example: tf.io.parse_single_example(example, FEATURES_DICT))\n",
        "\n",
        "#     for image_features in parsedDataset:\n",
        "#         image_NDVI = image_features['NDVI'].numpy()\n",
        "#         image_ugs = image_features['ugs'].numpy()\n",
        "#         image_NIR = image_features['b4'].numpy()\n",
        "#         image_R = image_features['b3'].numpy()\n",
        "#         image_G = image_features['b2'].numpy()\n",
        "#         image_B = image_features['b1'].numpy()\n",
        "#         image_RGB = np.dstack((image_NIR, image_R, image_ugs))\n",
        "\n",
        "#         # write the image to file\n",
        "#         matplotlib.image.imsave(output_path + '/' + 'patch' + str(num)  + '.png', image_RGB)\n",
        "\n",
        "#         num = num + 1\n"
      ],
      "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 Urban Open Spaces (UOS) are representnted in a binary way (0=non-UOS, 1=UOS], a sigmoid activation function is suitable here."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-ZaAVy4nbtSL"
      },
      "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",
        "## Train\n",
        "You train a Keras model by calling `.fit()` on it.  Here we're going to train for 500 epochs (Takes around **4.5h** with a GPU on colab).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NzzaWxOhSxBy",
        "scrolled": true
      },
      "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()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0a08jjrVD6Xx"
      },
      "source": [
        "##Load a model and continue training"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j5AqngMAD_bs"
      },
      "source": [
        "EPOCHS_LEFT = 53\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='val_loss',\n",
        "    mode='min',\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": "2CRfTdeX50Gc"
      },
      "source": [
        "## Accuracy\n",
        "plot the train and validation loss as well as the train and validation accuracy"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KuIktBwu50Gd"
      },
      "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": "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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6nqoVrn7sLJp"
      },
      "source": [
        "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.\n",
        "\n",
        "Export and prediction is done by tile."
      ]
    },
    {
      "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 GEE FOLDER HERE.\n",
        "\n",
        "# Base file name to use for TFRecord files and assets.\n",
        "attica_image_base = '1171_UOS_'\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": "JHvpdlXRsjn3"
      },
      "source": [
        "## Export images"
      ]
    },
    {
      "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 Google Drive...')\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": "markdown",
      "metadata": {
        "id": "pMomKb2sPlHD"
      },
      "source": [
        "Run the export (all attica by tile, it takes **2-3h** on GEE)."
      ]
    },
    {
      "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",
        "\n",
        "  print('OK!')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QcL3dloG50Gd"
      },
      "source": [
        "## Load a saved model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U2XrwZHp66j4"
      },
      "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": "LQ8Mm3vI50Gd"
      },
      "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": "OoNafMURsc5A"
      },
      "source": [
        "## Prediction"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BeK6vPauPw4g"
      },
      "source": [
        "Run the prediction and upload to GEE (for whole attica, *prediction* takes **20m** and *upload* to GEE takes around **... h**).\n",
        "\n",
        "If the GEE image upload fails with the error \"Cannot read mixer\", try the specific tiles again, by typing their numbers in the \"tilesList\" variable and changing the for loop to v2."
      ]
    },
    {
      "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}\n",
        "\n",
        "print('DONE!')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KxACnxKFrQ_J",
        "scrolled": true
      },
      "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",
        "#####################################################\n",
        "## FOR ALL TILES --> (V1)\n",
        "\n",
        "# first tile\n",
        "tile_number = int(1)\n",
        "\n",
        "# loop\n",
        "for i in range(1,(len(splitted.getInfo())+1)): # V1\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",
        "\n",
        "# #####################################################\n",
        "# ## FOR SPECIFIC TILES --> (V2) (if failed to run)\n",
        "\n",
        "# tilesList = [14,29,34,36,45] # if you want to run for specific tiles\n",
        "\n",
        "# for i in tilesList: # v2 if you want to run for specific tiles\n",
        "\n",
        "#   tile_number = int(i)\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"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pQ2fju_1ssXe"
      },
      "source": [
        "END OF NOTEBOOK"
      ]
    }
  ]
}