Demonstration_of_DSen2_Superresolution_.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Demonstration of DSen2 Superresolution .ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H2t2au9hMITi",
        "colab_type": "text"
      },
      "source": [
        "#Demonstration of DSen2 Superresolution \n",
        "*Johannes Mast*\n",
        "\n",
        "\n",
        "*30.03.2020*\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lLupyujYarFr",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## Introduction\n",
        "This script is an implementation of *Super-resolution of Sentinel-2 images: Learning a globally applicable deep neural network* [[Link]](https://arxiv.org/abs/1803.04271) by Charis Lanaras, José Bioucas-Dias, Silvano Galliani, Emmanuel Baltsavias and Konrad Schindler. The code was made available online on [github](https://github.com/lanha/DSen2). \n",
        "\n",
        "This implementation is designed to do three things:\n",
        "\n",
        "1. **Demonstrate** the functionality of the code. The code published in the abovementioned Github is designed to be run from shell and -while documented- is fairly opaque. In this script, we want to execute the individual parts step by step to gain a better understanding of the algorithm. \n",
        "\n",
        "2. **Introduce** the basic building blocks for flexible **semi-automisation** of the procedure. We develop functionality for semi-automated data acquisition and batch-processing.\n",
        "\n",
        "3. **Modify** the original code, in which some variables are hardcoded, to be more flexible in execution. This will also help achieving our other goals.\n",
        "\n",
        "The implementation is designed to run in google colab, connected to the users google drive. Google colab allows for the connection to a moderately powerful cloud runtime - This allows the user more flexibility in using a comparatively demanding algorithm.\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CL0f6UgrNZi5",
        "colab_type": "text"
      },
      "source": [
        "## Overview\n",
        "\n",
        "\n",
        "Sentinel 2 imagery comes in three different resolutions: 10m, 20m, and 60m. To obtain a data cube in which all bands are at the highest resolution, the 20m and 60m bands need to be superresolved. In this script, we will see how deep neural networks can be applied to this purpose. We begin with an overview of the model, the training procedure, and the possible outputs.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0pOnSPN6aiQj",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "###The Model\n",
        "\n",
        "The model's network architecture consists mainly of convolutional layers,\n",
        "ReLU non-linearities and skip connections packed into a series of residual blocks (ResBlocks) and notably a long skip connection from the input directly to the output.\n",
        "\n",
        "\n",
        "![alt text](https://ars.els-cdn.com/content/image/1-s2.0-S0924271618302636-gr5.sml)\n",
        "\n",
        "The model exists in a deep version **DSen2** and a very deep **VDSen2** version. Quoting the publication:\n",
        "\n",
        ">*VDSen2 has a lot\n",
        "higher capacity, and was designed with maximum accuracy\n",
        "in mind. It is closer in terms of size and training time to\n",
        "modern high-end CNNs for other image analysis tasks (Simonyan and Zisserman, 2015; He et al., 2016; Huang et al.,\n",
        "2017), but is approximately two times slower and five times\n",
        "slower in both training and prediction respectively, compared to its shallower counterpart (DSen2) [...] On the one hand, the very deep variant is\n",
        "consistently a bit better, while training and applying it is not\n",
        "more difficult, if adequate resources (i.e., high-end GPUs)\n",
        "are available.*\n",
        "\n",
        "The model comes with two different sets of weights, one which is trained on the superresolution from 60m to 10m (*030*) and one which is trained on the superresolution from 20m to 10m (*032*) .\n",
        "\n",
        "The weights for the VDSen2 (very deep) models should be able to be downloaded from  [s2_033_lr_1e-04.hdf5](http://n.ethz.ch/~lanarasc/DSen2/s2_033_lr_1e-04.hdf5) (*033*) and [s2_034_lr_1e-04.hdf5](http://n.ethz.ch/~lanarasc/DSen2/s2_034_lr_1e-04.hdf5) (*034*) but are offline as of the creation of this notebook.\n",
        "\n",
        "In this experiment, we will not use the VDSen2 model, as the much lighter DSen2 produces only slightly worse results.\n",
        "\n",
        "For a detailed description of the models we refer to the publication. We will also look at the model code more closely in the Prediction section in Chapter 3.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V8S3BcJTaZ99",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "### Training\n",
        "\n",
        "The available model comes with already trained weights. Training is beneficial but not necessary. To quote the publication:\n",
        " \n",
        "\n",
        "> *These shall enable out-of-the-box super-resolution of Sentinel-2 images world-wide, with minimal knowledge\n",
        "of neural network tools. Of course, if a study is focussed\n",
        "only in a specific geographic location, biome or land-cover\n",
        "type, even better result can be expected by training the network only with images showing those specific conditions.\n",
        "The literature suggests that in that case, it may be best\n",
        "to start from our globally trained network and fine-tune it\n",
        "through further training iterations on task-specific imagery*\n",
        "\n",
        "For the training, no labeled data is necessary. This is possible because -the developers assume- the superresolution is scale invariant. Therefore, for the network can learn the superresolution from downsampled versions of image bands to their initial resolution. No labeled training data is needed, and potentially any Sentinel-2 tile can be used for training.\n",
        "\n",
        "For a detailed description of the training procedure which produced the provided weights, we again refer to the publication. A list of the Sentinel-2 images which were used for training can be found [as a supplement to the code on the github](https://github.com/lanha/DSen2/blob/master/S2_tiles_training.txt). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dTfYBLiZac3r",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "### Bands \n",
        "\n",
        "The principles of the algorithm can be applied to any data, however, the publication and the DSen2 model deal specifically with the superresolution of Sentinel-2 bands. Let us take a quick look at these bands.\n",
        "\n",
        "**10m Resolution:**\n",
        "\n",
        "These are visible and VNIR bands. They are already at the highest possible resolution of 10m and will not be predicted, although they can be copied to the output file. \n",
        "\n",
        "![alt text](https://earth.esa.int/image/image_gallery?uuid=c5fa6c3e-2978-4fb8-ac95-3be9c5171be2&groupId=247904&t=1345630320883)\n",
        "\n",
        "\n",
        "**20m Resolution:**\n",
        "\n",
        "These are VNIR bands 5,6,7, 8a, 11 and 12. They are prime candidates for the superresolution and can all be predicted at 10m.\n",
        "![alt text](https://earth.esa.int/image/image_gallery?uuid=15dad96b-be6a-4b04-931d-d8c4db39e9e2&groupId=247904&t=1345630328076)\n",
        "\n",
        "\n",
        "**60m Resolution:**\n",
        "\n",
        "These are bands 1, 9 and 10 which are mainly used to measure water vapor, aerosols and clouds. Bands 1 and 9 can optionally be predicted to 10m. Band 10 is too noisy and cannot be superresolved.\n",
        "\n",
        "\n",
        "![alt text](https://earth.esa.int/image/image_gallery?uuid=f6117fbe-1513-4a84-acc4-845e14e5c876&groupId=247904&t=1345630315020)\n",
        "\n",
        "\n",
        "\n",
        "Note that due to python indexing, the band indices within the code are usually shifted by one position compared to the documentation except for gdal and rasterio, which start indexing at one. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KVLJW0G-XPrE",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "### Structure of this script\n",
        "\n",
        "The process can be clearly structured into a number of tasks:\n",
        "\n",
        "* **Setup**:  Firstly, we have to setup the environment. This must be done anytime the VM is started.\n",
        "*   **Data Acquisition**: Secondly, we acquire imagery for our desired study area by using [sentinelsat](https://pypi.org/project/sentinelsat/) to query and download Sentinel-2 images from the [copernicus open access hub](https://pypi.org/project/sentinelsat/).\n",
        "*  **Training**: The functionality for training the network is given by DSen2, it requires, however, a comparatively large amount of data and time. It is therefore **optional**.\n",
        "*  **Superresolution**: We subsequently use the implementation of DSen2 to superresolve the desired bands of the images we have downloaded. The results will be saved on the drive.\n",
        "* **Visualisation**: Of course we want take a quick look at the result.\n",
        "* **Streamlining**: To finally apply the algorithm to many files, we combine its essential parts back into a reduced function. We can apply this function in a loop to process a list of files.\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8TNth4wUX3kk",
        "colab_type": "text"
      },
      "source": [
        "## The Script"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iEWcjo4lZDSK",
        "colab_type": "text"
      },
      "source": [
        "### Chapter 1: Setup\n",
        "\n",
        "Before we begin, we must set up the environment. The VMs come with some preinstalled packages, but other elements must be reestablished anytime we connect to a new VM."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Sb8ZT5YuUsod",
        "colab_type": "text"
      },
      "source": [
        "#### Connecting to Google Drive\n",
        "\n",
        "Every time the notebook runtime is started, it must be reconnected to the Google Drive.\n",
        "Therefore, this must be done at the beginning of every script. An authorisation code is required as an input by the user."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CoXmNYqZKngt",
        "colab_type": "code",
        "outputId": "0a4b4896-ac7d-400b-971a-6e62701017b3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 110
        }
      },
      "source": [
        "%%capture\n",
        "from google.colab import drive\n",
        "from importlib import reload\n",
        "\n",
        "drive.mount('/content/gdrive',force_remount=True)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&response_type=code&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly\n",
            "\n",
            "Enter your authorization code:\n",
            "··········\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Yz9YxFQTU7x_",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "We want to create all required directories, as well as import the DSen2 and the input data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TGsPi0KcPFLJ",
        "colab_type": "text"
      },
      "source": [
        "#### Project Directory\n",
        "\n",
        "We create a project directory *DSen2*, which will contain everything which is related to this project. Alternatively we can navigate to an already existing directory on the drive, which has been created beforehand (for instance, the directory containing the present notebook)\n",
        "\n",
        "We further create a subdirectory *Data* in which our data will be saved."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dgBw8x4eL-UR",
        "colab_type": "code",
        "outputId": "7412ea81-320c-483b-e16c-51b8cc29149d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 127
        }
      },
      "source": [
        "#Change to the google drive main\n",
        "%cd \"gdrive/My Drive\"\n",
        "\n",
        "#Create a new directory\n",
        "!mkdir DSen2Main\n",
        "\n",
        "#Enter the new directory\n",
        "%cd DSen2Main\n",
        "\n",
        "#make a dir for the data\n",
        "!mkdir Data\n",
        "!mkdir Data/Training\n",
        "#make an output dir\n",
        "!mkdir Outputs"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/content/gdrive/My Drive\n",
            "mkdir: cannot create directory ‘DSen2Main’: File exists\n",
            "/content/gdrive/My Drive/DSen2Main\n",
            "mkdir: cannot create directory ‘Data’: File exists\n",
            "mkdir: cannot create directory ‘Data/Training’: File exists\n",
            "mkdir: cannot create directory ‘Outputs’: File exists\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1qkNt103PP4q",
        "colab_type": "text"
      },
      "source": [
        "#### Importing the DSen2\n",
        "\n",
        "The DSen2 is available on github. It forms the backbone of the algorithm, and comes with pretrained weights. We clone it into our project directory."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J5yGH2_IKVX7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git clone https://github.com/lanha/DSen2"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HWKaHyN_ReSn",
        "colab_type": "text"
      },
      "source": [
        "#### Installing required packages\n",
        "Packages can be installed using the *!pip install* command. This needs to be done everytime a notebook is loaded which requires certain dependencies. Fortunately, many packages are already installed on the VMs. We install several common packages for handling geospatial data, which should only take a minute.\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xB3RNE8GRdwa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%capture\n",
        "!pip install rasterio\n",
        "!pip install geopandas\n",
        "!pip install descartes\n",
        "!pip install gdal\n",
        "!pip install sentinelsat\n",
        "!pip install geojson\n",
        "!pip install shapely\n",
        "!pip install Pillow"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cwmgHkj24-wH",
        "colab_type": "text"
      },
      "source": [
        "#### Selecting  Tensorflow Version\n",
        "\n",
        "Colab VMs offer different versions of tensorflow - here, we select tensorflow 1.15.0, the default as of 27.03.2020 is tensorflow 2."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Matuuq0e5Ren",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%tensorflow_version 1.x\n",
        "import tensorflow as tf\n",
        "print(\"Tensorflow version \" + tf.__version__)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SPNwfG7poXwE",
        "colab_type": "text"
      },
      "source": [
        "### Chapter 2: Data Acquisition \n",
        "\n",
        "We use the [sentinelsat](https://pypi.org/project/sentinelsat/) package to connect to the api and then query and download images. For the spatial query we use a polygon string of our study area. We will be looking for three images in this example, but the number could be increased or reduced without issue."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_IBlNgBGFtwa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "n_dl_images=3"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tmTyKk4JnmbS",
        "colab_type": "text"
      },
      "source": [
        "#### Defining the AOI"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "redouJUhQoOA",
        "colab_type": "text"
      },
      "source": [
        "For querying spatially, we use an AOI. The author of this script is German, and thinks his country is very representative of the world, so we choose some German coordinates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D-UJQkKi_rxm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "lon_min=9\n",
        "lon_max=10\n",
        "lat_min= 48\n",
        "lat_max= 52"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IHznury5fmkk",
        "colab_type": "text"
      },
      "source": [
        "#### Connecting to the API\n",
        "To query the images, we connect to the API via the sentinelsat package. For this, we need our username and a password."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SDwJnEDXkll1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sentinelsat.sentinel import SentinelAPI, read_geojson, geojson_to_wkt\n",
        "from datetime import date\n",
        "import getpass\n",
        "# connect to the API\n",
        "api = SentinelAPI(getpass.getpass(prompt=\"Please enter Copernicus Open Access Hub username\"),\n",
        "                  getpass.getpass(prompt=\"Please enter Copernicus Open Access Hub password\"),\n",
        "                  'https://scihub.copernicus.eu/dhus')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zGzirr_nAGZ-",
        "colab_type": "text"
      },
      "source": [
        "We create an polygon from our AOI coordinates and visualize it on a map to make sure we got the coordinates right."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VInlwjTHWHYc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "aoi = 'POLYGON((%s %s,%s %s,%s %s,%s %s,%s %s))'  %(lon_min,lat_min,lon_min,lat_max,lon_max,lat_max,lon_max,lat_min,lon_min,lat_min)\n",
        "import plotly.graph_objects as go\n",
        "fig = go.Figure(go.Scattermapbox(\n",
        "    fill = \"toself\",text=\"Query Area\",\n",
        "    lon = [lon_min, lon_max, lon_max, lon_min],\n",
        "    lat = [lat_max, lat_max, lat_min, lat_min],\n",
        "    marker = { 'size': 10, 'color': \"red\" }))\n",
        "fig.update_layout(\n",
        "     title=\"AOI for the Query\",font=dict(family=\"Arial\",size=18,),\n",
        "     mapbox = {\n",
        "        'style': \"stamen-terrain\",\n",
        "        'center': {'lon': ((lon_min+lon_max)/2), 'lat': ((lat_min+lat_max)/2) },\n",
        "        'zoom': 4},\n",
        "    showlegend = False)\n",
        "fig.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tv8fabESQ72F",
        "colab_type": "text"
      },
      "source": [
        "Looks alright!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PxpRLXDefwJ5",
        "colab_type": "text"
      },
      "source": [
        "#### Querying the images\n",
        "\n",
        "Now we are ready to query the API.\n",
        "We use our polygon to query spatially and specify a date range to query temporally. We also only want scenes with low cloud cover. For our example run, we have that liberty."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jkhGsKV0rV5l",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# search by polygon, time, and SciHub query keywords\n",
        "products = api.query(aoi,\n",
        "                     date=('20190819', date(2019, 12, 29)),\n",
        "                     platformname='Sentinel-2',\n",
        "                     cloudcoverpercentage=(0, 30))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GwWGiV4IcjYF",
        "colab_type": "text"
      },
      "source": [
        "We convert the query result to a dataframe, sort it by cloud cover and take the best 3 results."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-gnuGHATvBec",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "# convert to Pandas DataFrame\n",
        "products_df = api.to_dataframe(products)\n",
        "\n",
        "# sort and limit to first n_dl_images sorted products\n",
        "products_df_sorted = products_df.sort_values(['cloudcoverpercentage', 'ingestiondate'], ascending=[True, True])\n",
        "products_df_sorted = products_df_sorted.head(n_dl_images)\n",
        "products_df_sorted"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VpP2G7PsdfDJ",
        "colab_type": "text"
      },
      "source": [
        "#### Download\n",
        "\n",
        "Then we start the download into the data directory. This should be relatively quick, only taking a few minutes. Note: Uncomment `%%capture` to hide progress"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9cOZZFQBpcXk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# download sorted and reduced products\n",
        "#%%capture \n",
        "api.download_all(products_df_sorted.index,\"Data/\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VSh9l2r_PCgD",
        "colab_type": "text"
      },
      "source": [
        "Check the downloaded files - We should have one for each of the selected query results (3 in our example)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bf5gGFR9ohUU",
        "colab_type": "code",
        "outputId": "da9bea7e-c759-4e98-9449-61d37918d552",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 72
        }
      },
      "source": [
        "import glob\n",
        "downloaded_files=(glob.glob(\"Data/*.zip\"))\n",
        "downloaded_files"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['Data/S2B_MSIL1C_20190823T103029_N0208_R108_T32UNB_20190823T124349.zip',\n",
              " 'Data/S2B_MSIL1C_20190823T103029_N0208_R108_T32UMB_20190823T124349.zip',\n",
              " 'Data/S2B_MSIL1C_20190823T103029_N0208_R108_T32UNA_20190823T124349.zip']"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j3TiRBHzPE4v",
        "colab_type": "text"
      },
      "source": [
        "Extracting the files not necessary, as the Gdal and therefore DSen2 can also accept zipped files.  We want to do it here regardless, because for the training procedure we need the extracted SAFE files. \n",
        "\n",
        "Since it also takes a while, and can occasionally break things, we do want to do it if we are not doing training. In that case, we can comment the following codeblock.\n",
        "\n",
        "If we extract, we then list the MSI xml files which we use instead of the zip files- If it worked, we see one file for each previous zip file, again, 3 in this example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ug10BMXy4b3c",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import zipfile\n",
        "for i in range(0, len(downloaded_files)):\n",
        "  with zipfile.ZipFile(downloaded_files[i], 'r') as zip_ref:\n",
        "      zip_ref.extractall(\"Data\")\n",
        "downloaded_files=(glob.glob(\"Data/**/*MSI*.xml\",recursive=True))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4TL1u6vsOuFI",
        "colab_type": "text"
      },
      "source": [
        "###Chapter 3: Training\n",
        "\n",
        "To perform the superresolution, our network needs some weights. Here, we have a number of options:\n",
        "\n",
        "* We can use the **pretrained weights** [provided by the developers](https://github.com/lanha/DSen2/tree/f6b8790c28a136b0ee57f8d4aa801348efcb4b74/models) , which have been extensively trained on a large number of Sentinel-2 images from many different geographical areas. These are perfectly fine to use! And so this entire chapter can be skipped.\n",
        "\n",
        "\n",
        "But maybe we are dealing with a particularly unusual region, for which the network was not trained. Then, it might make sense to train the network to also work well on this region. \n",
        "To fine-tune our model to our AOI, we can train it on the images we have just downloaded! \n",
        "\n",
        "* We can **train** from scratch, in which case the weights are randomly initialized using the [HeUniform](https://arxiv.org/abs/1502.01852) method.\n",
        "\n",
        "* It is more advisable, however, to **transfer** from  the pretrained weights, as that will give us a very solid base to start from, and will therefore greatly reduce the required training time.\n",
        "\n",
        "In this script, we will, as an example, train the 030 model which comes pretrained for the superresolution of 60m bands to 10m resolution, with two of our previously downloaded images.\n",
        "\n",
        "\n",
        "\n",
        ">  At this point, a note: As of the time of the writing of this script, the code from the developers github implements rather strictly the procedure of the published study.\n",
        "This creates reproducibility for the study, but also means that many variables, parameters, and especially path names are hardcoded. This conflicts with the goals of our script to flexibly execute the procedure flexibly and in a reasonable time. We therefore have to make some changes to the code. And for doing that, the simple way is to execute it in code blocks. So this chapter will be fairly bulky."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PigxTCbsHDZi",
        "colab_type": "text"
      },
      "source": [
        "We need to accomplish a number of tasks:\n",
        "\n",
        "*   **Creating the patches**: From our previously downloaded files, we create a number of smaller image-patches.\n",
        "*   **Creating the Train/Test split**: We create an index file, which indexes the previously created patches. Using this file, we can split the data into a training and a testing partition.\n",
        "*   **Training**: We train the 030 Network using the patches.\n",
        "*   **Replacing the original file**: We replace the original weight file. We do this to avoid any unnecessary changes in file paths for the remainder of this script.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MtQHDHWG1V3w",
        "colab_type": "text"
      },
      "source": [
        "#### Creating the Patch Dataset\n",
        "\n",
        "First, we sample a number of patches from our image, which we can then use for training.\n",
        "To quote the publication:\n",
        ">Sentinel-2 images are too big to fit them into GPU memory for training and testing, and in fact it is unlikely than\n",
        "long-range context over distances of a kilometer or more\n",
        "plays any significant role for super-resolution at the 10 m\n",
        "level. With this in mind, we train the network on small\n",
        "patches of w×h = (32×32) for T2×, respectively (96×96)\n",
        "pixels for S6×. We note that this corresponds to a receptive field of several hundred metres on the ground, sufficient to capture the local low-level texture and potentially\n",
        "also small semantic structures such as individual buildings\n",
        "or small waterbodies, but not large-scale topographic features. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VkaWvm4cH8YX",
        "colab_type": "text"
      },
      "source": [
        "First, we set a couple of parameters.\n",
        "\n",
        "Essential are:\n",
        "\n",
        "*   **data_file**: The Sentinel-2 file from which patches are created. This can be either the original ZIP file, or the S2A[...].xml file in a SAFE directory extracted from that ZIP. When running from command line, this is the only necessary argument. However, it only accepts one file at a time. We change it so that now it accepts a list of files (see **data_files** below)\n",
        "\n",
        "No output file is specified, the code automatically saves into a test folder. A path can be specified with **save_prefix** below.\n",
        "\n",
        "Optional are:\n",
        "\n",
        "* **roi_x_y**: A string which sets the region of interest (x_1,y_1,x_2,y_2) to extract as pixels locations on the 10m bands.\n",
        "* **test_data**: If `true` stores test patches in a separate dir? Not necessary for us.\n",
        "* **save_prefix**: We can use this to add a prefix to the output file. This could also be a directory, as in `\"Test/Outputs/\"`. If we execute from command line, the default is to save into `\"../data/\"`, but we change this to` Data/Training/`.\n",
        "* **write_images**: If `true`, write PNG images for the original and the superresolved bands, together with a composite rgb image (first three 10m bands), all with a quick and dirty clipping to 99%% of the original bands dynamic range and a quantization of the values to 256 levels. Since our space on the drive is limited, let's not do that for now.\n",
        "* **run_60**: If `true`, creates patches also from the 60m channels. This is what we want to do, so we set this to True.\n",
        "* **true_data**: If `true`, creates patches for S2 without GT. This option is not really useful here, please check the testing folder for predicting S2 images.\n",
        "\n",
        "We added/changed the following parameters:\n",
        "\n",
        "*   **data_files**: The Sentinel-2 files from which patches are created. We get the path to the downloaded and extracted SAFE files, but only take the first two.\n",
        "\n",
        "* **NR_CROP**: In the provided implementation, the number of patches is given by the default of NR_CROP in `save_random_patches`, which is 8000. To control the number of patches, we add this variable, which is passed on towards `save_random_patches`.\n",
        "The result will be a file containing the patches which is output to the specified folder. We use 1500 here - a comparatively low number, but we want to save time.\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9rRpE3Suso4L",
        "colab_type": "code",
        "outputId": "0f50469d-f294-4d11-d966-5274fa82829a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 54
        }
      },
      "source": [
        "training_files=glob.glob(\"Data/*SAFE\",recursive=True)[0:2]\n",
        "training_files"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['Data/S2B_MSIL1C_20190823T103029_N0208_R108_T32UNB_20190823T124349.SAFE',\n",
              " 'Data/S2B_MSIL1C_20190823T103029_N0208_R108_T32UMB_20190823T124349.SAFE']"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EkGGOF604wlO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "data_files=training_files\n",
        "test_data=False\n",
        "roi_x_y=\"\"\n",
        "save_prefix=\"Data/Training/\"\n",
        "write_images=False\n",
        "run_60=True\n",
        "true_data=False\n",
        "NR_CROP=1500"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "plENYtciZDJu",
        "colab_type": "text"
      },
      "source": [
        "The bulk of the work is done by the function below. We have modified it slightly compared to the source code, to include the `NR_CROP` argument."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8SGx2zJJ3cDA",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "#@title\n",
        "\n",
        "from __future__ import division\n",
        "import argparse\n",
        "import numpy as np\n",
        "from osgeo import gdal\n",
        "import sys\n",
        "from collections import defaultdict\n",
        "import re\n",
        "import os\n",
        "import imageio\n",
        "import json\n",
        "sys.path.append('../')\n",
        "from DSen2.utils.patches import downPixelAggr, save_test_patches, save_random_patches, save_random_patches60, save_test_patches60\n",
        "\n",
        "data_filename = '/MTD_MSIL1C.xml'\n",
        "\n",
        "# sleep(randint(0, 20))\n",
        "\n",
        "def readS2fromFile(data_file,\n",
        "                   test_data=False,\n",
        "                   roi_x_y=None,\n",
        "                   save_prefix=\"../data/\",\n",
        "                   write_images=False,\n",
        "                   run_60=False,\n",
        "                   true_data=False,\n",
        "                   NR_CROP=8000):\n",
        "\n",
        "    if run_60:\n",
        "        select_bands = \"B1,B2,B3,B4,B5,B6,B7,B8,B8A,B9,B11,B12\"\n",
        "    else:\n",
        "        select_bands = \"B2,B3,B4,B5,B6,B7,B8,B8A,B11,B12\"\n",
        "\n",
        "    raster = gdal.Open(data_file + data_filename)\n",
        "\n",
        "    datasets = raster.GetSubDatasets()\n",
        "    tenMsets = []\n",
        "    twentyMsets = []\n",
        "    sixtyMsets = []\n",
        "    unknownMsets = []\n",
        "    for (dsname, dsdesc) in datasets:\n",
        "        if '10m resolution' in dsdesc:\n",
        "            tenMsets += [ (dsname, dsdesc) ]\n",
        "        elif '20m resolution' in dsdesc:\n",
        "            twentyMsets += [ (dsname, dsdesc) ]\n",
        "        elif '60m resolution' in dsdesc:\n",
        "            sixtyMsets += [ (dsname, dsdesc) ]\n",
        "        else:\n",
        "            unknownMsets += [ (dsname, dsdesc) ]\n",
        "\n",
        "    if roi_x_y:\n",
        "        roi_x1, roi_y1, roi_x2, roi_y2 = [float(x) for x in re.split(',', args.roi_x_y)]\n",
        "\n",
        "    # case where we have several UTM in the data set\n",
        "    # => select the one with maximal coverage of the study zone\n",
        "    utm_idx = 0\n",
        "    utm = \"\"\n",
        "    all_utms = defaultdict(int)\n",
        "    xmin, ymin, xmax, ymax = 0, 0, 0, 0\n",
        "    largest_area = -1\n",
        "    # process even if there is only one 10m set, in order to get roi -> pixels\n",
        "    for (tmidx, (dsname, dsdesc)) in enumerate(tenMsets + unknownMsets):\n",
        "        ds = gdal.Open(dsname)\n",
        "        if roi_x_y:\n",
        "            tmxmin = max(min(roi_x1, roi_x2, ds.RasterXSize - 1), 0)\n",
        "            tmxmax = min(max(roi_x1, roi_x2, 0), ds.RasterXSize - 1)\n",
        "            tmymin = max(min(roi_y1, roi_y2, ds.RasterYSize - 1), 0)\n",
        "            tmymax = min(max(roi_y1, roi_y2, 0), ds.RasterYSize - 1)\n",
        "            # enlarge to the nearest 60 pixel boundary for the super-resolution\n",
        "            tmxmin = int(tmxmin / 36) * 36\n",
        "            tmxmax = int((tmxmax + 1) / 36) * 36 - 1\n",
        "            tmymin = int(tmymin / 36) * 36\n",
        "            tmymax = int((tmymax + 1) / 36) * 36 - 1\n",
        "        else:\n",
        "            tmxmin = 0\n",
        "            tmxmax = ds.RasterXSize - 1\n",
        "            tmymin = 0\n",
        "            tmymax = ds.RasterYSize - 1\n",
        "\n",
        "        area = (tmxmax - tmxmin + 1) * (tmymax - tmymin + 1)\n",
        "        current_utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "        if area > all_utms[current_utm]:\n",
        "            all_utms[current_utm] = area\n",
        "        if area > largest_area:\n",
        "            xmin, ymin, xmax, ymax = tmxmin, tmymin, tmxmax, tmymax\n",
        "            largest_area = area\n",
        "            utm_idx = tmidx\n",
        "            utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "\n",
        "    # convert comma separated band list into a list\n",
        "    select_bands = [x for x in re.split(',',select_bands) ]\n",
        "\n",
        "    print(\"Selected UTM Zone:\".format(utm))\n",
        "    print(\"Selected pixel region: xmin=%d, ymin=%d, xmax=%d, ymax=%d:\" % (xmin, ymin, xmax, ymax))\n",
        "    print(\"Selected pixel region: tmxmin=%d, tmymin=%d, tmxmax=%d, tmymax=%d:\" % (tmxmin, tmymin, tmxmax, tmymax))\n",
        "    print(\"Image size: width=%d x height=%d\" % (xmax - xmin + 1, ymax - ymin + 1))\n",
        "\n",
        "    if xmax < xmin or ymax < ymin:\n",
        "        print(\"Invalid region of interest / UTM Zone combination\")\n",
        "        sys.exit(0)\n",
        "\n",
        "    selected_10m_data_set = None\n",
        "    if not tenMsets:\n",
        "        selected_10m_data_set = unknownMsets[0]\n",
        "    else:\n",
        "        selected_10m_data_set = tenMsets[utm_idx]\n",
        "    selected_20m_data_set = None\n",
        "    for (dsname, dsdesc) in enumerate(twentyMsets):\n",
        "        if utm in dsdesc:\n",
        "            selected_20m_data_set = (dsname, dsdesc)\n",
        "    # if not found, assume the listing is in the same order\n",
        "    # => OK if only one set\n",
        "    if not selected_20m_data_set: selected_20m_data_set = twentyMsets[utm_idx]\n",
        "    selected_60m_data_set = None\n",
        "    for (dsname, dsdesc) in enumerate(sixtyMsets):\n",
        "        if utm in dsdesc:\n",
        "            selected_60m_data_set = (dsname, dsdesc)\n",
        "    if not selected_60m_data_set: selected_60m_data_set = sixtyMsets[utm_idx]\n",
        "\n",
        "    ds10 = gdal.Open(selected_10m_data_set[0])\n",
        "    ds20 = gdal.Open(selected_20m_data_set[0])\n",
        "    ds60 = gdal.Open(selected_60m_data_set[0])\n",
        "\n",
        "    def validate_description(description):\n",
        "        m = re.match(\"(.*?), central wavelength (\\d+) nm\", description)\n",
        "        if m:\n",
        "            return m.group(1) + \" (\" + m.group(2) + \" nm)\"\n",
        "        # Some HDR restrictions... ENVI band names should not include commas\n",
        "\n",
        "        pos = description.find(',')\n",
        "        return description[:pos] + description[(pos + 1):]\n",
        "\n",
        "    def get_band_short_name(description):\n",
        "        if ',' in description:\n",
        "            return description[:description.find(',')]\n",
        "        if ' ' in description:\n",
        "            return description[:description.find(' ')]\n",
        "        return description[:3]\n",
        "\n",
        "    validated_10m_bands = []\n",
        "    validated_10m_indices = []\n",
        "    validated_20m_bands = []\n",
        "    validated_20m_indices = []\n",
        "    validated_60m_bands = []\n",
        "    validated_60m_indices = []\n",
        "    validated_descriptions = defaultdict(str)\n",
        "\n",
        "    sys.stdout.write(\"Selected 10m bands:\")\n",
        "    for b in range(0, ds10.RasterCount):\n",
        "        desc = validate_description(ds10.GetRasterBand(b + 1).GetDescription())\n",
        "        shortname = get_band_short_name(desc)\n",
        "        if shortname in select_bands:\n",
        "            sys.stdout.write(\" \" + shortname)\n",
        "            select_bands.remove(shortname)\n",
        "            validated_10m_bands += [shortname]\n",
        "            validated_10m_indices += [b]\n",
        "            validated_descriptions[shortname] = desc\n",
        "    sys.stdout.write(\"\\nSelected 20m bands:\")\n",
        "    for b in range(0, ds20.RasterCount):\n",
        "        desc = validate_description(ds20.GetRasterBand(b + 1).GetDescription())\n",
        "        shortname = get_band_short_name(desc)\n",
        "        if shortname in select_bands:\n",
        "            sys.stdout.write(\" \" + shortname)\n",
        "            select_bands.remove(shortname)\n",
        "            validated_20m_bands += [shortname]\n",
        "            validated_20m_indices += [b]\n",
        "            validated_descriptions[shortname] = desc\n",
        "    sys.stdout.write(\"\\nSelected 60m bands:\")\n",
        "    for b in range(0, ds60.RasterCount):\n",
        "        desc = validate_description(ds60.GetRasterBand(b + 1).GetDescription())\n",
        "        shortname = get_band_short_name(desc)\n",
        "        if shortname in select_bands:\n",
        "            sys.stdout.write(\" \" + shortname)\n",
        "            select_bands.remove(shortname)\n",
        "            validated_60m_bands += [shortname]\n",
        "            validated_60m_indices += [b]\n",
        "            validated_descriptions[shortname] = desc\n",
        "    sys.stdout.write(\"\\n\")\n",
        "\n",
        "    if validated_10m_indices:\n",
        "        print(\"Loading selected data from: %s\" % selected_10m_data_set[1])\n",
        "        data10 = np.rollaxis(\n",
        "            ds10.ReadAsArray(xoff=xmin, yoff=ymin, xsize=xmax - xmin + 1, ysize=ymax - ymin + 1, buf_xsize=xmax - xmin + 1,\n",
        "                             buf_ysize=ymax - ymin + 1), 0, 3)[:, :, validated_10m_indices]\n",
        "\n",
        "    if validated_20m_indices:\n",
        "        print(\"Loading selected data from: %s\" % selected_20m_data_set[1])\n",
        "        data20 = np.rollaxis(\n",
        "            ds20.ReadAsArray(xoff=xmin // 2, yoff=ymin // 2, xsize=(xmax - xmin + 1) // 2, ysize=(ymax - ymin + 1) // 2,\n",
        "                             buf_xsize=(xmax - xmin + 1) // 2, buf_ysize=(ymax - ymin + 1) // 2), 0, 3)[:, :,\n",
        "                 validated_20m_indices]\n",
        "\n",
        "    if validated_60m_indices:\n",
        "        print(\"Loading selected data from: %s\" % selected_60m_data_set[1])\n",
        "        data60 = np.rollaxis(\n",
        "            ds60.ReadAsArray(xoff=xmin // 6, yoff=ymin // 6, xsize=(xmax - xmin + 1) // 6, ysize=(ymax - ymin + 1) // 6,\n",
        "                             buf_xsize=(xmax - xmin + 1) // 6, buf_ysize=(ymax - ymin + 1) // 6), 0, 3)[:, :,\n",
        "                 validated_60m_indices]\n",
        "\n",
        "    # The percentile_data argument is used to plot superresolved and original data\n",
        "    # with a comparable black/white scale\n",
        "    def save_band(data, name, percentile_data=None):\n",
        "        if percentile_data is None:\n",
        "            percentile_data = data\n",
        "        mi, ma = np.percentile(percentile_data, (1, 99))\n",
        "        band_data = np.clip(data, mi, ma)\n",
        "        band_data = (band_data - mi) / (ma - mi)\n",
        "        imageio.imsave(save_prefix + name + \".png\", band_data)  # img_as_uint(band_data))\n",
        "\n",
        "    chan3 = data10[:, :, 0]\n",
        "    vis = (chan3 < 1).astype(np.int)\n",
        "    if np.sum(vis) > 0:\n",
        "        print('The selected image has some blank pixels')\n",
        "        # sys.exit()\n",
        "\n",
        "    scale20 = 2\n",
        "    scale60 = 6\n",
        "\n",
        "    data10_gt = data10\n",
        "    data20_gt = data20\n",
        "\n",
        "    if not true_data:\n",
        "        if run_60:\n",
        "            data60_gt = data60\n",
        "            data10_lr = downPixelAggr(data10_gt, SCALE=scale60)\n",
        "            data20_lr = downPixelAggr(data20_gt, SCALE=scale60)\n",
        "            data60_lr = downPixelAggr(data60_gt, SCALE=scale60)\n",
        "        else:\n",
        "            data10_lr = downPixelAggr(data10_gt, SCALE=scale20)\n",
        "            data20_lr = downPixelAggr(data20_gt, SCALE=scale20)\n",
        "            if scale20 > 2:\n",
        "                data20_lr = downPixelAggr(data20_gt, SCALE=scale20//2)\n",
        "\n",
        "    if data_file.endswith('/'):\n",
        "        tmp = os.path.split(data_file)[0]\n",
        "        data_file = os.path.split(tmp)[1]\n",
        "    else:\n",
        "        data_file = os.path.split(data_file)[1]\n",
        "    print(data_file)\n",
        "\n",
        "    if test_data:\n",
        "        if run_60:\n",
        "            out_per_image0 = save_prefix + 'test60/'\n",
        "            out_per_image = save_prefix + 'test60/' + data_file + '/'\n",
        "        else:\n",
        "            out_per_image0 = save_prefix + 'test/'\n",
        "            out_per_image = save_prefix + 'test/' + data_file + '/'\n",
        "        if not os.path.isdir(out_per_image0):\n",
        "            os.mkdir(out_per_image0)\n",
        "        if not os.path.isdir(out_per_image):\n",
        "            os.mkdir(out_per_image)\n",
        "\n",
        "        print('Writing files for testing to:{}'.format(out_per_image))\n",
        "        if run_60:\n",
        "            save_test_patches60(data10_lr, data20_lr, data60_lr, out_per_image)\n",
        "            with open(out_per_image + 'roi.json', 'w') as f:\n",
        "                json.dump([tmxmin // scale60, tmymin // scale60, (tmxmax + 1) // scale60, (tmymax + 1) // scale60], f)\n",
        "        else:\n",
        "            save_test_patches(data10_lr, data20_lr, out_per_image)\n",
        "            with open(out_per_image + 'roi.json', 'w') as f:\n",
        "                json.dump([tmxmin // scale20, tmymin // scale20, (tmxmax+1) // scale20, (tmymax+1) // scale20], f)\n",
        "\n",
        "        if not os.path.isdir(out_per_image + 'no_tiling/'):\n",
        "            os.mkdir(out_per_image + 'no_tiling/')\n",
        "\n",
        "        print(\"Now saving the whole image without tiling...\")\n",
        "        if run_60:\n",
        "            np.save(out_per_image + 'no_tiling/' + 'data60_gt', data60_gt.astype(np.float32))\n",
        "            np.save(out_per_image + 'no_tiling/' + 'data60', data60_lr.astype(np.float32))\n",
        "        else:\n",
        "            np.save(out_per_image + 'no_tiling/' + 'data20_gt', data20_gt.astype(np.float32))\n",
        "            save_band(data10_lr[:, :, 0:3], '/test/' + data_file + '/RGB')\n",
        "        np.save(out_per_image + 'no_tiling/' + 'data10', data10_lr.astype(np.float32))\n",
        "        np.save(out_per_image + 'no_tiling/' + 'data20', data20_lr.astype(np.float32))\n",
        "\n",
        "    elif write_images:\n",
        "        print('Creating RGB images...')\n",
        "        save_band(data10_lr[:, :, 0:3], '/raw/rgbs/' + data_file + 'RGB')\n",
        "        save_band(data20_lr[:, :, 0:3], '/raw/rgbs/' + data_file + 'RGB20')\n",
        "\n",
        "    elif true_data:\n",
        "        out_per_image0 = save_prefix + 'true/'\n",
        "        out_per_image = save_prefix + 'true/' + data_file + '/'\n",
        "        if not os.path.isdir(out_per_image0):\n",
        "            os.mkdir(out_per_image0)\n",
        "        if not os.path.isdir(out_per_image):\n",
        "            os.mkdir(out_per_image)\n",
        "\n",
        "        print('Writing files for testing to:{}'.format(out_per_image))\n",
        "        save_test_patches60(data10_gt, data20_gt, data60_gt, out_per_image, patchSize=384, border=12)\n",
        "\n",
        "        with open(out_per_image + 'roi.json', 'w') as f:\n",
        "            json.dump([tmxmin, tmymin, tmxmax+1, tmymax+1], f)\n",
        "\n",
        "        if not os.path.isdir(out_per_image + 'no_tiling/'):\n",
        "            os.mkdir(out_per_image + 'no_tiling/')\n",
        "\n",
        "        print(\"Now saving the whole image without tiling...\")\n",
        "        np.save(out_per_image + 'no_tiling/' + 'data10', data10_gt.astype(np.float32))\n",
        "        np.save(out_per_image + 'no_tiling/' + 'data20', data20_gt.astype(np.float32))\n",
        "        np.save(out_per_image + 'no_tiling/' + 'data60', data60_gt.astype(np.float32))\n",
        "\n",
        "    else:\n",
        "        if run_60:\n",
        "            out_per_image0 = save_prefix + 'train60/'\n",
        "            out_per_image = save_prefix + 'train60/' + data_file + '/'\n",
        "        else:\n",
        "            out_per_image0 = save_prefix + 'train/'\n",
        "            out_per_image = save_prefix + 'train/' + data_file + '/'\n",
        "        if not os.path.isdir(out_per_image0):\n",
        "            os.mkdir(out_per_image0)\n",
        "        if not os.path.isdir(out_per_image):\n",
        "            os.mkdir(out_per_image)\n",
        "        print('Writing files for training to:{}'.format(out_per_image))\n",
        "        if run_60:\n",
        "            save_random_patches60(data60_gt, data10_lr, data20_lr, data60_lr, out_per_image,NR_CROP)\n",
        "        else:\n",
        "            save_random_patches(data20_gt, data10_lr, data20_lr, out_per_image,NR_CROP)\n",
        "\n",
        "    print(\"Success.\")\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9iB8Kzcib2F_",
        "colab_type": "text"
      },
      "source": [
        "We run the function in a loop over our files. Since we just have two here it is a bit superflous, but we can prove that it works in principle."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oxRqL97YYG3r",
        "colab_type": "code",
        "outputId": "7f55f68e-06f3-420d-c3e9-a122f12cfd13",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 678
        }
      },
      "source": [
        "for i in range(0, len(data_files)):\n",
        "  data_file=data_files[i]\n",
        "  readS2fromFile(data_file,\n",
        "                test_data,\n",
        "                roi_x_y,\n",
        "                save_prefix,\n",
        "                write_images,\n",
        "                run_60,\n",
        "                true_data,\n",
        "                NR_CROP)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Selected UTM Zone:\n",
            "Selected pixel region: xmin=0, ymin=0, xmax=10979, ymax=10979:\n",
            "Selected pixel region: tmxmin=0, tmymin=0, tmxmax=10979, tmymax=10979:\n",
            "Image size: width=10980 x height=10980\n",
            "Selected 10m bands: B4 B3 B2 B8\n",
            "Selected 20m bands: B5 B6 B7 B8A B11 B12\n",
            "Selected 60m bands: B1 B9\n",
            "Loading selected data from: Bands B2, B3, B4, B8 with 10m resolution, UTM 32N\n",
            "Loading selected data from: Bands B5, B6, B7, B8A, B11, B12 with 20m resolution, UTM 32N\n",
            "Loading selected data from: Bands B1, B9, B10 with 60m resolution, UTM 32N\n",
            "S2B_MSIL1C_20190823T103029_N0208_R108_T32UNB_20190823T124349.SAFE\n",
            "Writing files for training to:Data/Training/train60/S2B_MSIL1C_20190823T103029_N0208_R108_T32UNB_20190823T124349.SAFE/\n",
            "(1500, 2, 96, 96)\n",
            "(1500, 4, 96, 96)\n",
            "(1500, 6, 48, 48)\n",
            "(1500, 2, 16, 16)\n",
            "Done!\n",
            "Success.\n",
            "Selected UTM Zone:\n",
            "Selected pixel region: xmin=0, ymin=0, xmax=10979, ymax=10979:\n",
            "Selected pixel region: tmxmin=0, tmymin=0, tmxmax=10979, tmymax=10979:\n",
            "Image size: width=10980 x height=10980\n",
            "Selected 10m bands: B4 B3 B2 B8\n",
            "Selected 20m bands: B5 B6 B7 B8A B11 B12\n",
            "Selected 60m bands: B1 B9\n",
            "Loading selected data from: Bands B2, B3, B4, B8 with 10m resolution, UTM 32N\n",
            "Loading selected data from: Bands B5, B6, B7, B8A, B11, B12 with 20m resolution, UTM 32N\n",
            "Loading selected data from: Bands B1, B9, B10 with 60m resolution, UTM 32N\n",
            "S2B_MSIL1C_20190823T103029_N0208_R108_T32UMB_20190823T124349.SAFE\n",
            "Writing files for training to:Data/Training/train60/S2B_MSIL1C_20190823T103029_N0208_R108_T32UMB_20190823T124349.SAFE/\n",
            "(1500, 2, 96, 96)\n",
            "(1500, 4, 96, 96)\n",
            "(1500, 6, 48, 48)\n",
            "(1500, 2, 16, 16)\n",
            "Done!\n",
            "Success.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RDxGL4anQ8PM",
        "colab_type": "text"
      },
      "source": [
        "####Indexing and Splitting the Dataset\n",
        "\n",
        "We create a simple index of the patches we have, and split it randomly.\n",
        "\n",
        "Normally this is done by the `create_random.py ` which uses hardcoded paths. To control the paths, we execute it here as a codeblock, and change the path to where we saved previously.\n",
        "\n",
        "Variables which we must use here are:\n",
        "\n",
        "\n",
        "*   **size**: The number of files we processed (one, in our example) multiplied by the number of patches per tile (same as NR_CROP in the previous section)\n",
        "*   **ratio**: The ratio of validation files. We set it to the suggested default of 0.1\n",
        "* **path**: This is the path to the directory which contains the files we created in the previous block. So for us, it is the same as our previous `save_prefix`, which was `Data/Training/`\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DoJ0xjfEHrH2",
        "colab_type": "code",
        "outputId": "3cad78d4-3898-416c-b165-2cc28a6f4599",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 72
        }
      },
      "source": [
        "from random import randrange\n",
        "# Size: number of S2 tiles (times) patches per tile\n",
        "size = len(data_files)*NR_CROP\n",
        "ratio = .1\n",
        "nb = int(size * ratio)\n",
        "path = save_prefix\n",
        "\n",
        "\n",
        "index = np.zeros(size).astype(np.bool)\n",
        "i = 0\n",
        "while np.sum(index.astype(np.int)) < nb:\n",
        "    x = randrange(0, size)\n",
        "    index[x] = True\n",
        "    i += 1\n",
        "    \n",
        "if run_60:\n",
        "  path=path+\"train60/\"\n",
        "else:\n",
        "  path=path+\"train/\"\n",
        "np.save(path + 'val_index.npy', index)\n",
        "\n",
        "print('Full no of samples: {}'.format(size))\n",
        "print('Validation samples: {}'.format(np.sum(index.astype(np.int))))\n",
        "\n",
        "print(\"Number of iterations: {}\".format(i))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Full no of samples: 3000\n",
            "Validation samples: 300\n",
            "Number of iterations: 320\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z1Al1KSPVHXL",
        "colab_type": "text"
      },
      "source": [
        "We should now see how many samples we have that we can use for training. Two images with 450 patches each results in 900 patches sampled."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y_XZSTVScPOr",
        "colab_type": "text"
      },
      "source": [
        "####Training\n",
        "Now we have the resources we need for the training itself.\n",
        "We can run from the `supres_train.py`, but again, we crack open the code to examine it and to reduce the number of epochs."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PR-Pe47WLGcz",
        "colab_type": "text"
      },
      "source": [
        "Setting variables\n",
        "\n",
        "*   **predict_file**: If we provide a file here, we can do a single prediction. The purpose for this is presumably for testing.\n",
        "*   **resume_file**:  Path to the previously created weights file to resume from. We want to train the 030 model for 60m data, which is why we use the `s2_030_lr_1e-05.hdf5` from the model folder - otherwise we would get the 032 model `s2_032_lr_1e-04.hdf5`\n",
        "*   **true**: Use true scale data? No simulation or different resolutions.\n",
        "*   **run_60**: Whether to run a 60->10m network. Default 20->10m.\n",
        "*   **deep**: train a deep network. Takes too long, lets not do that here.\n",
        "\n",
        "We added here one more variable:\n",
        "\n",
        "* **n_epochs**: The number of epochs used for training. Fixed to 8192 in the provided implementation. For our purposes of testing, 6 is enough.\n",
        "\n",
        "At this point, we could also adjust the following parameters:\n",
        "\n",
        "\n",
        "*  **model_nr**: The name of the output model. Must be 7 characters. Here, we go with the default of \"s2_038_\"\n",
        "*   **SCALE**: A factor by which the raw reflectance values are divided by for \"numerical stability\".\n",
        "*   **lr**: The base learning rate.\n",
        "\n",
        "We will leave those unchanged."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IpYbxFTNepqS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "predict_file=None\n",
        "if(run_60):\n",
        "  resume_file=\"DSen2/models/s2_030_lr_1e-05.hdf5\"\n",
        "else:\n",
        "  resume_file=\"DSen2/models/s2_032_lr_1e-04.hdf5\"\n",
        "true=True\n",
        "run_60=run_60 #from previous code block\n",
        "deep=False\n",
        "path=\"Data/Training/\"\n",
        "n_epochs=5"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YAk765JEeoKq",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "from __future__ import division\n",
        "import numpy as np\n",
        "import datetime\n",
        "import glob\n",
        "import time\n",
        "import argparse\n",
        "import os\n",
        "import sys\n",
        "import matplotlib as mpl\n",
        "mpl.use('Agg')\n",
        "from keras.optimizers import Nadam\n",
        "from keras.callbacks import ModelCheckpoint, Callback, ReduceLROnPlateau\n",
        "from keras.utils import plot_model\n",
        "import keras.backend as K\n",
        "\n",
        "sys.path.append('../')\n",
        "from DSen2.utils.patches import recompose_images, OpenDataFilesTest, OpenDataFiles\n",
        "from DSen2.utils.DSen2Net import s2model\n",
        "\n",
        "K.set_image_data_format('channels_first')\n",
        "\n",
        "# Define file prefix for new training, must be 7 characters of this form:\n",
        "model_nr = 's2_038_'\n",
        "SCALE = 2000\n",
        "lr = 1e-4\n",
        "\n",
        "\n",
        "path = path\n",
        "if not os.path.isdir(path):\n",
        "    os.mkdir(path)\n",
        "out_path = path+\"network_data/\"\n",
        "if not os.path.isdir(out_path):\n",
        "    os.mkdir(out_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Sq5_AJRMxMll",
        "colab_type": "text"
      },
      "source": [
        "We also create a callback which allows us to track how our losses are progressing - this will help us judge the performance of our network."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0tAVWBDFJgo9",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "#@title\n",
        "class PlotLosses(Callback):\n",
        "    def __init__(self, model_nr, lr):\n",
        "        self.model_nr = model_nr\n",
        "        self.lr = lr\n",
        "\n",
        "    def on_train_begin(self, logs=None):\n",
        "        self.losses = []\n",
        "        self.val_losses = []\n",
        "        self.i = 0\n",
        "        self.x = []\n",
        "        self.filename = out_path + self.model_nr + '_lr_{:.1e}.txt'.format(self.lr)\n",
        "        open(self.filename, 'w').close()\n",
        "\n",
        "    def on_epoch_end(self, epoch, logs=None):\n",
        "        import matplotlib.pyplot as plt\n",
        "        plt.ioff()\n",
        "\n",
        "        lr = float(K.get_value(self.model.optimizer.lr))\n",
        "        # data = np.loadtxt(\"training.log\", skiprows=1, delimiter=',')\n",
        "        self.losses.append(logs.get('loss'))\n",
        "        self.val_losses.append(logs.get('val_loss'))\n",
        "        self.x.append(self.i)\n",
        "        self.i += 1\n",
        "        try:\n",
        "            with open(self.filename, 'a') as self.f:\n",
        "                self.f.write('Finished epoch {:5d}: loss {:.3e}, valid: {:.3e}, lr: {:.1e}\\n'\n",
        "                             .format(epoch, logs.get('loss'), logs.get('val_loss'), lr))\n",
        "\n",
        "            if epoch > 500:\n",
        "                plt.clf()\n",
        "                plt.plot(self.x[475:], self.losses[475:], label='loss')\n",
        "                plt.plot(self.x[475:], self.val_losses[475:], label='val_loss')\n",
        "                plt.legend()\n",
        "                plt.xlabel('epochs')\n",
        "                # plt.waitforbuttonpress(0)\n",
        "                plt.savefig(out_path + self.model_nr + '_loss4.png')\n",
        "            elif epoch > 250:\n",
        "                plt.clf()\n",
        "                plt.plot(self.x[240:], self.losses[240:], label='loss')\n",
        "                plt.plot(self.x[240:], self.val_losses[240:], label='val_loss')\n",
        "                plt.legend()\n",
        "                plt.xlabel('epochs')\n",
        "                # plt.waitforbuttonpress(0)\n",
        "                plt.savefig(out_path + self.model_nr + '_loss3.png')\n",
        "            elif epoch > 100:\n",
        "                plt.clf()\n",
        "                plt.plot(self.x[85:], self.losses[85:], label='loss')\n",
        "                plt.plot(self.x[85:], self.val_losses[85:], label='val_loss')\n",
        "                plt.legend()\n",
        "                plt.xlabel('epochs')\n",
        "                # plt.waitforbuttonpress(0)\n",
        "                plt.savefig(out_path + self.model_nr + '_loss2.png')\n",
        "            elif epoch > 50:\n",
        "                plt.clf()\n",
        "                plt.plot(self.x[50:], self.losses[50:], label='loss')\n",
        "                plt.plot(self.x[50:], self.val_losses[50:], label='val_loss')\n",
        "                plt.legend()\n",
        "                plt.xlabel('epochs')\n",
        "                # plt.waitforbuttonpress(0)\n",
        "                plt.savefig(out_path + self.model_nr + '_loss1.png')\n",
        "            else:\n",
        "                plt.clf()\n",
        "                plt.plot(self.x[0:], self.losses[0:], label='loss')\n",
        "                plt.plot(self.x[0:], self.val_losses[0:], label='val_loss')\n",
        "                plt.legend()\n",
        "                plt.xlabel('epochs')\n",
        "                # plt.waitforbuttonpress(0)\n",
        "                plt.savefig(out_path + self.model_nr + '_loss0.png')\n",
        "        except IOError:\n",
        "            print('Network drive unavailable.')\n",
        "            print(datetime.datetime.now().time())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vVRYp-7xqBiS",
        "colab_type": "text"
      },
      "source": [
        "First we create the model using the `s2model` function. \n",
        "We choose the [nadam](https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/Nadam) optimizer and compile our model.\n",
        "\n",
        "We load the weights from the previously set `resume_file`, and then the data.\n",
        "\n",
        "Then we can start the training! This should take about 3-5 minutes with our chosen settings, so now is a good time to take a break.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "deWof-2fgDgC",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "if path is not None:\n",
        "    path = path\n",
        "\n",
        "# input_shape = ((4,32,32),(6,16,16))\n",
        "if run_60:\n",
        "    input_shape = ((4, None, None), (6, None, None), (2, None, None))\n",
        "else:\n",
        "    input_shape = ((4, None, None), (6, None, None))\n",
        "# create model\n",
        "if deep:\n",
        "    model = s2model(input_shape, num_layers=32, feature_size=256)\n",
        "    batch_size = 8\n",
        "else:\n",
        "    model = s2model(input_shape, num_layers=6, feature_size=128)\n",
        "    batch_size = 128\n",
        "print('Symbolic Model Created.')\n",
        "\n",
        "nadam = Nadam(lr=lr,\n",
        "              beta_1=0.9,\n",
        "              beta_2=0.999,\n",
        "              epsilon=1e-8,\n",
        "              schedule_decay=0.004)\n",
        "              # clipvalue=0.000005)\n",
        "\n",
        "model.compile(optimizer=nadam, loss='mean_absolute_error', metrics=['mean_squared_error'])\n",
        "print('Model compiled.')\n",
        "model.count_params()\n",
        "# model.summary()\n",
        "\n",
        "if predict_file:\n",
        "    if true:\n",
        "        folder = 'true/'\n",
        "        border = 12\n",
        "    elif run_60:\n",
        "        folder = 'test60/'\n",
        "        border = 12\n",
        "    else:\n",
        "        folder = 'test/'\n",
        "        border = 4\n",
        "    model_nr = predict_file[-20:-13]\n",
        "    print('Changing the model number to: {}'.format(model_nr))\n",
        "    model.load_weights(args.predict_file)\n",
        "    print(\"Predicting using file: {}\".format(predict_file))\n",
        "    fileList = [os.path.basename(x) for x in sorted(glob.glob(path + folder + '*SAFE'))]\n",
        "    print(\"Using patches from the files:\")\n",
        "    print(fileList)\n",
        "    for dset in fileList:\n",
        "        start = time.time()\n",
        "        print(\"Timer started.\")\n",
        "        print(\"Predicting: {}.\".format(dset))\n",
        "        train, image_size = OpenDataFilesTest(path + folder + dset, run_60, SCALE, true)\n",
        "        prediction = model.predict(train,\n",
        "                                    batch_size=8,\n",
        "                                    verbose=1)\n",
        "        prediction_file = model_nr + '-predict'\n",
        "        # np.save(path + 'test/' + dset + '/' + prediction_file + 'pat', prediction * SCALE)\n",
        "        images = recompose_images(prediction, border=border, size=image_size)\n",
        "        print('Writing to file...')\n",
        "        np.save(path + folder + dset + '/' + prediction_file, images * SCALE)\n",
        "        end = time.time()\n",
        "        print('Elapsed time: {}.'.format(end - start))\n",
        "    sys.exit(0)\n",
        "\n",
        "if resume_file:\n",
        "    print(\"Will resume from the weights {}\".format(resume_file))\n",
        "    model.load_weights(resume_file)\n",
        "    model_nr = resume_file[-20:-13]\n",
        "    print('Changing the model number to: {}'.format(model_nr))\n",
        "\n",
        "else:\n",
        "    print('Model number is {}'.format(model_nr))\n",
        "    plot_model(model, to_file=out_path + model_nr+'model.png', show_shapes=True, show_layer_names=True)\n",
        "\n",
        "    model_yaml = model.to_yaml()\n",
        "    with open(out_path + model_nr + \"model.yaml\", 'w') as yaml_file:\n",
        "        yaml_file.write(model_yaml)\n",
        "\n",
        "filepath = out_path + model_nr + 'lr_{:.0e}.hdf5'.format(lr)\n",
        "checkpoint = ModelCheckpoint(filepath,\n",
        "                              monitor='val_loss',\n",
        "                              verbose=1,\n",
        "                              save_best_only=True,\n",
        "                              save_weights_only=False,\n",
        "                              mode='auto')\n",
        "plot_losses = PlotLosses(model_nr, lr)\n",
        "LRreducer = ReduceLROnPlateau(monitor='val_loss',\n",
        "                              factor=0.5,\n",
        "                              patience=5,\n",
        "                              verbose=1,\n",
        "                              epsilon=1e-6,\n",
        "                              cooldown=20,\n",
        "                              min_lr=1e-5)\n",
        "\n",
        "callbacks_list = [checkpoint, plot_losses, LRreducer]\n",
        "\n",
        "print('Loading the training data...from...')\n",
        "print(path)\n",
        "train, label, val_tr, val_lb = OpenDataFiles(path,  run_60, SCALE)\n",
        "print('Training starts...')\n",
        "\n",
        "model.fit(x=train,\n",
        "          y=label,\n",
        "          batch_size=batch_size,\n",
        "          epochs=n_epochs,\n",
        "          verbose=1,\n",
        "          callbacks=callbacks_list,\n",
        "          validation_split=0.,\n",
        "          validation_data=(val_tr, val_lb),\n",
        "          shuffle=True,\n",
        "          class_weight=None,\n",
        "          sample_weight=None,\n",
        "          initial_epoch=0,\n",
        "          validation_steps=None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q-en1lI_mQg3",
        "colab_type": "text"
      },
      "source": [
        "Via the callbacks, statistics are collected during training and finally plotted to an image. We can take a look at it to examing the progress."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RxsCeGu5nFVj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%matplotlib inline\n",
        "from IPython.display import Image\n",
        "Image(out_path+'s2_030__loss0.png')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KDn4Lr1_2icZ",
        "colab_type": "text"
      },
      "source": [
        "Hopefully, the validation loss is decreasing nicely, and then slowly flattening out. But since we use few epochs and little training data, it is not unlikely that some weirdness could also be happening."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2KFr8yLKndXb",
        "colab_type": "text"
      },
      "source": [
        "####Replacing the original\n",
        "\n",
        "If we are happy with the result, we replace the original weights file with the newly created one. Not very elegant - but it allows us to run the next section without any changes, regardless of whether we skipped the training or not."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JGHC0cB-n4Lg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from shutil import copyfile\n",
        "if(run_60):\n",
        "  copyfile(out_path+\"s2_030_lr_1e-04.hdf5\", \"DSen2/models/s2_030_lr_1e-05.hdf5\")\n",
        "else:\n",
        "  copyfile(out_path+\"s2_030_lr_1e-04.hdf5\", \"DSen2/models/s2_032_lr_1e-04.hdf5\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1HYe9VLhriwn",
        "colab_type": "text"
      },
      "source": [
        "Now we are done with training.\n",
        "\n",
        "What we have achieved in this section is this: We have tuned our model to perform better in the 60->10 for our downloaded tiles, and, therefore, our AOI. At least we hope that this is the case. In any case, now we can move on towards the superresolution process itself!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QQkULQv31zaC",
        "colab_type": "text"
      },
      "source": [
        "### Chapter 4: Superresolution\n",
        "\n",
        "Now we are ready to execute the superresolution algorithm.\n",
        "\n",
        "We want to superresolve both the 60m and the 20m bands, as well as keep the original 10m bands, to get a full stack at 10m resolution. To keep the processing time short, we will restrict ourselves to a small subset of the data.\n",
        "\n",
        "By default, DSen2 will be called via command line, with the command line arguments getting parsed to the function `s2_tiles_supres`.\n",
        "\n",
        "This does not tell us much about the inner workings, so we will dissect that function and execute its components step-by-step."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6dwyWpiDhwd8",
        "colab_type": "text"
      },
      "source": [
        "#### Setting up DSen2\n",
        "\n",
        "We import the DSen2 implementation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w6OmeNRqkeeX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import sys\n",
        "sys.path.append(\"DSen2/testing\")\n",
        "sys.path.append(\"DSen2\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YIorEUJpxQVr",
        "colab_type": "text"
      },
      "source": [
        "We import a couple of libraries which we need. DSen2-Supres itself depends on tensorflow, so this gets imported as well."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AMF4jIpZIrPk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from __future__ import division\n",
        "import argparse\n",
        "import numpy as np\n",
        "import os\n",
        "import re\n",
        "import sys\n",
        "from osgeo import gdal, osr\n",
        "from collections import defaultdict\n",
        "from supres import DSen2_20, DSen2_60"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nq16txaBxZQB",
        "colab_type": "text"
      },
      "source": [
        "#### Setting the parameters\n",
        "\n",
        "Here, we set a couple of parameters.\n",
        "Essential are:\n",
        "\n",
        "\n",
        "\n",
        "*   **data_file**: The file on which will be predicted. For S2 imagery, this filename should point to the xml file in the highest directory.\n",
        "*   **output_file**: The output filename. Tif is recommended, but other common types are supposedly implemented as well *(I did not test that so far)*. See **output_file_formats** below.\n",
        "\n",
        "Optional are:\n",
        "\n",
        "* **roi_lon_lat**: Sets the region of interest to extract, WGS84, decimal notation. lon_1,lat_1,lon_2,lat_2.\n",
        "* **roi_x_y**: Sets the region of interest to extract as pixels locations on the 10m bands. Use this syntax: x_1,y_1,x_2,y_2.  To speed things up, we use this parameter to only work on a small subset of the data (1500 by 1500 pixels).\n",
        "* **list_bands**: If `True`, lists bands in the input file subdata set matching the selected UTM zone. When run via command line, it then exits, but thats not what we do here.\n",
        "* **run_60**: If `True`, super-resolve the 20m **and** 60m bands(B1,B2,B3,B4,B5,B6,B7,B8,B8A,B9,B11,B12). Otherwise, only super-resolve the 20m bands (B2,B3,B4,B5,B6,B7,B8,B8A,B11,B12). Note: \"Band B10 is too noisy and is not super-resolved\".\n",
        "* **list_UTM**: If `True`, list all UTM zones present in the input file, together with their coverage of the ROI in 10m x 10m pixels.\n",
        "* **select_UTM**: Manually select a UTM zone to use - otherwise the algorithm will select the UTM zone with the best coverage.\n",
        "* **list_output_file_formats**: If `true`, lists all output file formats supported by GDAL. \n",
        "* **output_file_format**: Specifies the name of a GDAL driver that supports file creation, like ENVI or GTiff. We use the proven GTiff here.\n",
        "\n",
        "* **copy_original_bands**: Copy the original 10m bands into the new, predicted file? We set this to `true`, because then we have more bands to compare, but it obviously means the output files will be much larger. \n",
        "* **save_prefix**: We can use this to add a prefix to the output file. This could also be a directory, as in `\"Test/Outputs/\"`. Here, we use the previously created output folder.\n",
        "\n",
        "\n",
        "In this example, we use mostly defaults to explore the most important aspects of the algorithm, first only for one file.\n",
        "\n",
        "We also add one more switch that is not provided by the implementation:\n",
        "* **deep**: If we have acquired the VDSen2 models, and have placed them in the `DSen2/models` directory we can use this switch to use them for prediction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-uWxVkcaxXnM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "data_file=downloaded_files[0]\n",
        "output_file=\"GuterOutput.Tif\"\n",
        "roi_lon_lat=\"\"\n",
        "roi_x_y=\"1,1,1500,1500\"\n",
        "list_bands=True\n",
        "run_60=True\n",
        "list_UTM=False\n",
        "select_UTM=\"\"\n",
        "list_output_file_formats=True\n",
        "output_file_format=\"GTiff\"\n",
        "copy_original_bands=True\n",
        "save_prefix=\"Outputs/\"\n",
        "\n",
        "deep= False"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7KJyd4L6GSJx",
        "colab_type": "text"
      },
      "source": [
        "#### Listing possible output file formats\n",
        "\n",
        "Originally, this option will list the possible formats in which the result can be output and then quits. We can execute the codeblock to see what options we have. In this experiment we have chosen GTiff already, so this does not matter much, but it is good to know the alternatives. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9pbRiiaNQnrQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if list_output_file_formats:\n",
        "    dcount = gdal.GetDriverCount()\n",
        "    for didx in range(dcount):\n",
        "        driver = gdal.GetDriver(didx)\n",
        "        if driver:\n",
        "            metadata = driver.GetMetadata()\n",
        "        if (gdal.DCAP_CREATE in (driver and metadata) and metadata[gdal.DCAP_CREATE] == 'YES' and\n",
        "        gdal.DCAP_RASTER in metadata and metadata[gdal.DCAP_RASTER] == 'YES'):\n",
        "            name = driver.GetDescription()\n",
        "            if \"DMD_LONGNAME\" in metadata:\n",
        "                name += \": \" + metadata[\"DMD_LONGNAME\"]\n",
        "            else:\n",
        "                name = driver.GetDescription()\n",
        "            if \"DMD_EXTENSIONS\" in metadata: name += \" (\" + metadata[\"DMD_EXTENSIONS\"] + \")\"\n",
        "            print(name)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A7B1PXOFoE2c",
        "colab_type": "text"
      },
      "source": [
        "As we can see, there is quite a long list of alternatives."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YSzsvPjGI8uH",
        "colab_type": "text"
      },
      "source": [
        "#### Choosing the bands for predictions\n",
        "\n",
        "By the parameter **run60** to True, we have chosen to superresolve also the 60m bands in addition to the 20m bands. We create a list of all those bands.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DzL4aRCqGOXA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if run_60:\n",
        "    select_bands = 'B1,B2,B3,B4,B5,B6,B7,B8,B8A,B9,B11,B12'\n",
        "else:\n",
        "    select_bands = 'B2,B3,B4,B5,B6,B7,B8,B8A,B11,B12'\n",
        "\n",
        "# convert comma separated band list into a list\n",
        "select_bands = [x for x in re.split(',', select_bands)]\n",
        "select_bands"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9D1r1AZN-cR2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if roi_lon_lat:\n",
        "    roi_lon1, roi_lat1, roi_lon2, roi_lat2 = [float(x) for x in re.split(',', roi_lon_lat)]\n",
        "else:\n",
        "    roi_lon1, roi_lat1, roi_lon2, roi_lat2 = -180, -90, 180, 90\n",
        "\n",
        "if roi_x_y:\n",
        "    roi_x1, roi_y1, roi_x2, roi_y2 = [float(x) for x in re.split(',', roi_x_y)]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zHYYOFvpJYlO",
        "colab_type": "text"
      },
      "source": [
        "#### Loading the Data\n",
        "We load in the data into a couple of lists, one for each of the resolutions\n",
        "\n",
        "*   10m\n",
        "*   20m\n",
        "*  60m\n",
        "* Unknown\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "0scvyOqVkuJn",
        "colab": {}
      },
      "source": [
        "raster = gdal.Open(data_file)\n",
        "datasets = raster.GetSubDatasets();\n",
        "tenMsets = []\n",
        "twentyMsets = []\n",
        "sixtyMsets = []\n",
        "unknownMsets = []\n",
        "for (dsname, dsdesc) in datasets:\n",
        "    if '10m resolution' in dsdesc:\n",
        "        tenMsets += [(dsname, dsdesc)]\n",
        "    elif '20m resolution' in dsdesc:\n",
        "        twentyMsets += [(dsname, dsdesc)]\n",
        "    elif '60m resolution' in dsdesc:\n",
        "        sixtyMsets += [(dsname, dsdesc)]\n",
        "    else:\n",
        "        unknownMsets += [(dsname, dsdesc)]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7hUQfK1dJtid",
        "colab_type": "text"
      },
      "source": [
        "#### Choosing the subset correct UTM zone\n",
        "\n",
        "Choosing the utm zone is more complex than it seems.\n",
        "First, we must account for the possible restriction of the AOI with the **roi_lon_lat** and **roi_x_y** arguments.\n",
        "Then, we choose the UTM zone with the best coverage.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sGZ9jg3xT0dR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# case where we have several UTM in the data set\n",
        "# => select the one with maximal coverage of the study zone\n",
        "utm_idx = 0\n",
        "utm = select_UTM\n",
        "all_utms = defaultdict(int)\n",
        "xmin, ymin, xmax, ymax = 0, 0, 0, 0\n",
        "largest_area = -1\n",
        "# process even if there is only one 10m set, in order to get roi -> pixels\n",
        "for (tmidx, (dsname, dsdesc)) in enumerate(tenMsets + unknownMsets):\n",
        "    ds = gdal.Open(dsname)\n",
        "    if roi_x_y:\n",
        "        tmxmin = max(min(roi_x1, roi_x2, ds.RasterXSize - 1), 0)\n",
        "        tmxmax = min(max(roi_x1, roi_x2, 0), ds.RasterXSize - 1)\n",
        "        tmymin = max(min(roi_y1, roi_y2, ds.RasterYSize - 1), 0)\n",
        "        tmymax = min(max(roi_y1, roi_y2, 0), ds.RasterYSize - 1)\n",
        "        # enlarge to the nearest 60 pixel boundary for the super-resolution\n",
        "        tmxmin = int(tmxmin / 6) * 6\n",
        "        tmxmax = int((tmxmax + 1) / 6) * 6 - 1\n",
        "        tmymin = int(tmymin / 6) * 6\n",
        "        tmymax = int((tmymax + 1) / 6) * 6 - 1\n",
        "    elif not roi_lon_lat:\n",
        "        tmxmin = 0\n",
        "        tmxmax = ds.RasterXSize - 1\n",
        "        tmymin = 0\n",
        "        tmymax = ds.RasterYSize - 1\n",
        "    else:\n",
        "        xoff, a, b, yoff, d, e = ds.GetGeoTransform()\n",
        "        srs = osr.SpatialReference()\n",
        "        srs.ImportFromWkt(ds.GetProjection())\n",
        "        srsLatLon = osr.SpatialReference()\n",
        "        srsLatLon.SetWellKnownGeogCS(\"WGS84\");\n",
        "        ct = osr.CoordinateTransformation(srsLatLon, srs)\n",
        "\n",
        "\n",
        "        def to_xy(lon, lat):\n",
        "            (xp, yp, h) = ct.TransformPoint(lon, lat, 0.)\n",
        "            xp -= xoff\n",
        "            yp -= yoff\n",
        "            # matrix inversion\n",
        "            det_inv = 1. / (a * e - d * b)\n",
        "            x = (e * xp - b * yp) * det_inv\n",
        "            y = (-d * xp + a * yp) * det_inv\n",
        "            return (int(x), int(y))\n",
        "\n",
        "\n",
        "        x1, y1 = to_xy(roi_lon1, roi_lat1)\n",
        "        x2, y2 = to_xy(roi_lon2, roi_lat2)\n",
        "        tmxmin = max(min(x1, x2, ds.RasterXSize - 1), 0)\n",
        "        tmxmax = min(max(x1, x2, 0), ds.RasterXSize - 1)\n",
        "        tmymin = max(min(y1, y2, ds.RasterYSize - 1), 0)\n",
        "        tmymax = min(max(y1, y2, 0), ds.RasterYSize - 1)\n",
        "        # enlarge to the nearest 60 pixel boundary for the super-resolution\n",
        "        tmxmin = int(tmxmin / 6) * 6\n",
        "        tmxmax = int((tmxmax + 1) / 6) * 6 - 1\n",
        "        tmymin = int(tmymin / 6) * 6\n",
        "        tmymax = int((tmymax + 1) / 6) * 6 - 1\n",
        "    area = (tmxmax - tmxmin + 1) * (tmymax - tmymin + 1)\n",
        "    current_utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "    if area > all_utms[current_utm]:\n",
        "        all_utms[current_utm] = area\n",
        "    if current_utm == select_UTM:\n",
        "        xmin, ymin, xmax, ymax = tmxmin, tmymin, tmxmax, tmymax\n",
        "        utm_idx = tmidx\n",
        "        utm = current_utm\n",
        "        break\n",
        "    if area > largest_area:\n",
        "        xmin, ymin, xmax, ymax = tmxmin, tmymin, tmxmax, tmymax\n",
        "        largest_area = area\n",
        "        utm_idx = tmidx\n",
        "        utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "\n",
        "if list_UTM:\n",
        "    print(\"List of UTM zones (with ROI coverage in pixels):\")\n",
        "    for u in all_utms:\n",
        "        print(\"%s (%d)\" % (u, all_utms[u]))\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jhMJjy7GLcVu",
        "colab_type": "text"
      },
      "source": [
        "If successful, we now have chosen the correct subset of our AOI and, from that, best UTM zone for this prediction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JDm0z4uSLLkN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print(\"Selected UTM Zone:\", utm)\n",
        "print(\"Selected pixel region: xmin=%d, ymin=%d, xmax=%d, ymax=%d:\" % (xmin, ymin, xmax, ymax))\n",
        "print(\"Image size: width=%d x height=%d\" % (xmax - xmin + 1, ymax - ymin + 1))\n",
        "\n",
        "if xmax < xmin or ymax < ymin:\n",
        "    print(\"Invalid region of interest / UTM Zone combination\")\n",
        "    sys.exit(0)\n",
        "else:\n",
        "    print(\"all good\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jgaonaCq5pWe",
        "colab_type": "text"
      },
      "source": [
        "#### Opening the datasets bands\n",
        "\n",
        "We now use gdal to open the file and just load in the bands into three gdal datasets, one for each resolution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "caFTnPNoUKHK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "selected_10m_data_set = None\n",
        "if not tenMsets:\n",
        "    selected_10m_data_set = unknownMsets[0]\n",
        "else:\n",
        "    selected_10m_data_set = tenMsets[utm_idx]\n",
        "selected_20m_data_set = None\n",
        "for (dsname, dsdesc) in enumerate(twentyMsets):\n",
        "    if utm in dsdesc:\n",
        "        selected_20m_data_set = (dsname, dsdesc)\n",
        "# if not found, assume the listing is in the same order\n",
        "# => OK if only one set\n",
        "if not selected_20m_data_set: selected_20m_data_set = twentyMsets[utm_idx]\n",
        "selected_60m_data_set = None\n",
        "for (dsname, dsdesc) in enumerate(sixtyMsets):\n",
        "    if utm in dsdesc:\n",
        "        selected_60m_data_set = (dsname, dsdesc)\n",
        "if not selected_60m_data_set: selected_60m_data_set = sixtyMsets[utm_idx]\n",
        "\n",
        "ds10 = gdal.Open(selected_10m_data_set[0])\n",
        "ds20 = gdal.Open(selected_20m_data_set[0])\n",
        "ds60 = gdal.Open(selected_60m_data_set[0])\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OzBi4_Mf58de",
        "colab_type": "text"
      },
      "source": [
        "Since we set the **list_bands** to True, we now can get an output of all the bands, not just the bands we selected, as well as a description with their wavelength. Note the order, which is different from the original structure. This is also the order in which the results will be written into the output file.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i88O_1xwUfLa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def validate_description(description):\n",
        "    m = re.match(\"(.*?), central wavelength (\\d+) nm\", description)\n",
        "    if m:\n",
        "        return m.group(1) + \" (\" + m.group(2) + \" nm)\"\n",
        "    # Some HDR restrictions... ENVI band names should not include commas\n",
        "    if output_file_format == 'ENVI' and ',' in description:\n",
        "        pos = description.find(',')\n",
        "        return description[:pos] + description[(pos + 1):]\n",
        "    return description\n",
        "\n",
        "\n",
        "if list_bands:\n",
        "    print(\"\\n10m bands:\")\n",
        "    for b in range(0, ds10.RasterCount):\n",
        "        print(\"- \" + validate_description(ds10.GetRasterBand(b + 1).GetDescription()))\n",
        "    print(\"\\n20m bands:\")\n",
        "    for b in range(0, ds20.RasterCount):\n",
        "        print(\"- \" + validate_description(ds20.GetRasterBand(b + 1).GetDescription()))\n",
        "    print(\"\\n60m bands:\")\n",
        "    for b in range(0, ds60.RasterCount):\n",
        "        print(\"- \" + validate_description(ds60.GetRasterBand(b + 1).GetDescription()))\n",
        "    print(\"\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ZDYKey06efw",
        "colab_type": "text"
      },
      "source": [
        "#### Validating the bands\n",
        "\n",
        "The function then \"validates\" the bands by selecting only the bands which we selected via the **select_bands** and **run_60** parameters.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RqYgp2E7V94_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def get_band_short_name(description):\n",
        "    if ',' in description:\n",
        "        return description[:description.find(',')]\n",
        "    if ' ' in description:\n",
        "        return description[:description.find(' ')]\n",
        "    return description[:3]\n",
        "\n",
        "\n",
        "validated_10m_bands = []\n",
        "validated_10m_indices = []\n",
        "validated_20m_bands = []\n",
        "validated_20m_indices = []\n",
        "validated_60m_bands = []\n",
        "validated_60m_indices = []\n",
        "validated_descriptions = defaultdict(str)\n",
        "\n",
        "sys.stdout.write(\"Selected 10m bands:\")\n",
        "for b in range(0, ds10.RasterCount):\n",
        "    desc = validate_description(ds10.GetRasterBand(b + 1).GetDescription())\n",
        "    shortname = get_band_short_name(desc)\n",
        "    if shortname in select_bands:\n",
        "        sys.stdout.write(\" \" + shortname)\n",
        "        select_bands.remove(shortname)\n",
        "        validated_10m_bands += [shortname]\n",
        "        validated_10m_indices += [b]\n",
        "        validated_descriptions[shortname] = desc\n",
        "sys.stdout.write(\"\\nSelected 20m bands:\")\n",
        "for b in range(0, ds20.RasterCount):\n",
        "    desc = validate_description(ds20.GetRasterBand(b + 1).GetDescription())\n",
        "    shortname = get_band_short_name(desc)\n",
        "    if shortname in select_bands:\n",
        "        sys.stdout.write(\" \" + shortname)\n",
        "        select_bands.remove(shortname)\n",
        "        validated_20m_bands += [shortname]\n",
        "        validated_20m_indices += [b]\n",
        "        validated_descriptions[shortname] = desc\n",
        "sys.stdout.write(\"\\nSelected 60m bands:\")\n",
        "for b in range(0, ds60.RasterCount):\n",
        "    desc = validate_description(ds60.GetRasterBand(b + 1).GetDescription())\n",
        "    shortname = get_band_short_name(desc)\n",
        "    if shortname in select_bands:\n",
        "        sys.stdout.write(\" \" + shortname)\n",
        "        select_bands.remove(shortname)\n",
        "        validated_60m_bands += [shortname]\n",
        "        validated_60m_indices += [b]\n",
        "        validated_descriptions[shortname] = desc\n",
        "sys.stdout.write(\"\\n\")\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mEDpmvox8Ao-",
        "colab_type": "text"
      },
      "source": [
        "In this case, since we selected all the bands, the list of the bands is the same as before."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aqLUtzJE8LNn",
        "colab_type": "text"
      },
      "source": [
        "#### Setting the Output Filename\n",
        "\n",
        "We now set the output filename. If we had not given one, it would now create one here, with the name identical to the input file, and the default extension being .tif.\n",
        "\n",
        "At this point we also add the prefix, so that in our case, the full path points to the *Outputs* folder."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A-f1_J2oWJfQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# All query options are processed, we now require an output file\n",
        "if not output_file:\n",
        "    print(\"Error: you must provide the name of an output file. I will set it identical to the input...\")\n",
        "    output_file = os.path.split(data_file)[1] + '.tif'\n",
        "    # sys.exit(1)\n",
        "\n",
        "\n",
        "output_file = save_prefix + output_file\n",
        "# Some HDR restrictions... ENVI file name should be the .bin, not the .hdr\n",
        "if output_file_format == 'ENVI' and (output_file[-4:] == '.hdr' or output_file[-4:] == '.HDR'):\n",
        "    output_file = output_file[:-4] + '.bin'\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JQVVPQVkxHNv",
        "colab_type": "text"
      },
      "source": [
        "#### Prediction\n",
        "Any weights are loaded from the Dsen2/models folder. If we have \n",
        "\n",
        "Now we start the prediction for the desired bands - On the VM, this can take a few hours for a single image.\n",
        "\n",
        "The core is in the `DSen2_20` and `DSen2_60` functions which can be found in DSen2/testing/supres.py . We will not execute them step-by-step here, but we can take a look at the most important functions.\n",
        "\n",
        "We begin with `DSen2_20` which takes as inputs the 10m bands and 20m bands and can be used to construct both the deep and very deep versions of the network. \n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "**DSen2_20**\n",
        "```\n",
        "def DSen2_20(d10, d20, deep=False):\n",
        "    # Input to the funcion must be of shape:\n",
        "    #     d10: [x,y,4]      (B2, B3, B4, B8)\n",
        "    #     d20: [x/2,y/4,6]  (B5, B6, B7, B8a, B11, B12)\n",
        "    #     deep: specifies whether to use VDSen2 (True), or DSen2 (False)\n",
        "\n",
        "    border = 8\n",
        "    p10, p20 = get_test_patches(d10, d20, patchSize=128, border=border)\n",
        "    p10 /= SCALE\n",
        "    p20 /= SCALE\n",
        "    test = [p10, p20]\n",
        "    input_shape = ((4, None, None), (6, None, None))\n",
        "    prediction = _predict(test, input_shape, deep=deep)\n",
        "    images = recompose_images(prediction, border=border, size=d10.shape)\n",
        "    images *= SCALE\n",
        "    return images\n",
        "```\n",
        "\n",
        "\n",
        "\n",
        "We see that the `DSen2_20` function samples test patches using `get_test_patches` and passes them to the `_predict` function, which itself \n",
        "\n",
        "\n",
        "\n",
        "1.   Builds a `s2model` (see code cell below)\n",
        "2.   Loads the appropriate model weights from the DSen2/models folder using `load_weights`\n",
        "3.   Predicts the model using `predict` and returns the prediction\n",
        "\n",
        "\n",
        "**_predict**\n",
        "```\n",
        "\n",
        "def _predict(test, input_shape, deep=False, run_60=False):\n",
        "    # create model\n",
        "    if deep:\n",
        "        model = s2model(input_shape, num_layers=32, feature_size=256)\n",
        "        predict_file = MDL_PATH+'s2_034_lr_1e-04.hdf5' if run_60 else MDL_PATH+'s2_033_lr_1e-04.hdf5'\n",
        "    else:\n",
        "        model = s2model(input_shape, num_layers=6, feature_size=128)\n",
        "        predict_file = MDL_PATH+'s2_030_lr_1e-05.hdf5' if run_60 else MDL_PATH+'s2_032_lr_1e-04.hdf5'\n",
        "    print('Symbolic Model Created.')\n",
        "\n",
        "    model.load_weights(predict_file)\n",
        "    print(\"Predicting using file: {}\".format(predict_file))\n",
        "    prediction = model.predict(test, verbose=1)\n",
        "    return prediction\n",
        "```\n",
        "\n",
        "Unfortunately, because - as we can see here - the path to the model is hardcoded as relative path name to supres.py, we have to change the working directory to the testing folder to find the model weights.\n",
        "\n",
        "The model itself is built by the s2model function from Dsen2/utils/DSen2Net.py :\n",
        "\n",
        "\n",
        "**s2model**\n",
        "```\n",
        "def s2model(input_shape, num_layers=32, feature_size=256):\n",
        "\n",
        "    input10 = Input(shape=input_shape[0])\n",
        "    input20 = Input(shape=input_shape[1])\n",
        "    if len(input_shape) == 3:\n",
        "        input60 = Input(shape=input_shape[2])\n",
        "        x = Concatenate(axis=1)([input10, input20, input60])\n",
        "    else:\n",
        "        x = Concatenate(axis=1)([input10, input20])\n",
        "\n",
        "    # Treat the concatenation\n",
        "    x = Conv2D(feature_size, (3, 3), kernel_initializer='he_uniform', activation='relu', padding='same')(x)\n",
        "\n",
        "    for i in range(num_layers):\n",
        "        x = resBlock(x, feature_size)\n",
        "\n",
        "    # One more convolution, and then we add the output of our first conv layer\n",
        "    x = Conv2D(input_shape[-1][0], (3, 3), kernel_initializer='he_uniform', padding='same')(x)\n",
        "    # x = Dropout(0.3)(x)\n",
        "    if len(input_shape) == 3:\n",
        "        x = Add()([x, input60])\n",
        "        model = Model(inputs=[input10, input20, input60], outputs=x)\n",
        "    else:\n",
        "        x = Add()([x, input20])\n",
        "        model = Model(inputs=[input10, input20], outputs=x)\n",
        "    return model\n",
        "```\n",
        "We can see that the bulk of the architecture consists of a number of `resBlock`. These residual blocks consist of\n",
        "\n",
        "\n",
        "1.   Convolutional layer\n",
        "2.   ReLU layer\n",
        "3.   Convolutional layer\n",
        "4.   Residual Scaling\n",
        "\n",
        "![alt text](https://ars.els-cdn.com/content/image/1-s2.0-S0924271618302636-gr6.sml)\n",
        "\n",
        "These layers are supplemented by a skip connection straight from the input to the Output.\n",
        "```\n",
        "def resBlock(x, channels, kernel_size=[3, 3], scale=0.1):\n",
        "    tmp = Conv2D(channels, kernel_size, kernel_initializer='he_uniform', padding='same')(x)\n",
        "    tmp = Activation('relu')(tmp)\n",
        "    tmp = Conv2D(channels, kernel_size, kernel_initializer='he_uniform', padding='same')(tmp)\n",
        "    tmp = Lambda(lambda x: x * scale)(tmp)\n",
        "\n",
        "    return Add()([x, tmp])\n",
        "\n",
        "```\n",
        "\n",
        "Now that we have a basic idea of how the code works, we can move towards the prediction!\n",
        "\n",
        "First, we read the different bands into arrays `data10` `data20` and possibly `data60`. We rearrange the dimensions of the arrays (probably to fit tensorflow syntax).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "19ecZGGMWWH5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if validated_10m_indices:\n",
        "    print(\"Loading selected data from: %s\" % selected_10m_data_set[1])\n",
        "    data10 = np.rollaxis(\n",
        "        ds10.ReadAsArray(xoff=xmin, yoff=ymin, xsize=xmax - xmin + 1, ysize=ymax - ymin + 1, buf_xsize=xmax - xmin + 1,\n",
        "                         buf_ysize=ymax - ymin + 1), 0, 3)[:, :, validated_10m_indices]\n",
        "\n",
        "if validated_20m_indices:\n",
        "    print(\"Loading selected data from: %s\" % selected_20m_data_set[1])\n",
        "    data20 = np.rollaxis(\n",
        "        ds20.ReadAsArray(xoff=xmin // 2, yoff=ymin // 2, xsize=(xmax - xmin + 1) // 2, ysize=(ymax - ymin + 1) // 2,\n",
        "                         buf_xsize=(xmax - xmin + 1) // 2, buf_ysize=(ymax - ymin + 1) // 2), 0, 3)[:, :,\n",
        "             validated_20m_indices]\n",
        "\n",
        "if validated_60m_indices:\n",
        "    print(\"Loading selected data from: %s\" % selected_60m_data_set[1])\n",
        "    data60 = np.rollaxis(\n",
        "        ds60.ReadAsArray(xoff=xmin // 6, yoff=ymin // 6, xsize=(xmax - xmin + 1) // 6, ysize=(ymax - ymin + 1) // 6,\n",
        "                         buf_xsize=(xmax - xmin + 1) // 6, buf_ysize=(ymax - ymin + 1) // 6), 0, 3)[:, :,\n",
        "             validated_60m_indices]\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "utQQ9IR3vIgw",
        "colab_type": "text"
      },
      "source": [
        "As mentioned before, we have to change the working directory just while predicting. No big deal."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bxfeBp4WvGtP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%cd DSen2/testing/"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3JAMhEnRvFb5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if validated_60m_bands and validated_20m_bands and validated_10m_bands:\n",
        "    print(\"Super-resolving the 60m data into 10m bands\")\n",
        "    sr60 = DSen2_60(data10, data20, data60, deep=deep)\n",
        "else:\n",
        "    sr60 = None\n",
        "\n",
        "if validated_10m_bands and validated_20m_bands:\n",
        "    print(\"Super-resolving the 20m data into 10m bands\")\n",
        "    sr20 = DSen2_20(data10, data20, deep=deep)\n",
        "else:\n",
        "    sr20 = None\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FAXGkUUpz9Sh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%cd -"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w3a48efBfRQO",
        "colab_type": "text"
      },
      "source": [
        "Note that at the end we change back to our real working directory, the DSen20 Folder"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EcKE9Mqjl2bI",
        "colab_type": "text"
      },
      "source": [
        "#### Writing the output\n",
        "\n",
        "Now that the prediction is done, we can write the outputs. We have created the output folder and set the path to it earlier. \n",
        "If we use a geospatial output we use GDAL to create the output file and write the bands into that file.  If we have opted for the numpy-specific *npz* format, we use numpy to write the file.\n",
        "\n",
        "\n",
        "The order in which the bands are written is as follows:\n",
        "\n",
        "\n",
        "1.   Original 10m bands, if copied with **copy_original_bands**, in order `4 3 2 8`\n",
        "2.   Superresolved 20m bands, in order `5 6 7 8A 11 12 `\n",
        "3.   Superresolved 60m bands, in order `1 9`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LoHdVGFLJInC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "if output_file_format != \"npz\":\n",
        "    revert_to_npz = True\n",
        "    driver = gdal.GetDriverByName(output_file_format)\n",
        "    if driver:\n",
        "        metadata = driver.GetMetadata()\n",
        "        if gdal.DCAP_CREATE in metadata and metadata[gdal.DCAP_CREATE] == 'YES':\n",
        "            revert_to_npz = False\n",
        "    if revert_to_npz:\n",
        "        print(\"Gdal doesn't support creating %s files\" % output_file_format)\n",
        "        print(\"Writing to npz as a fallback\")\n",
        "        output_file_format = \"npz\"\n",
        "    bands = None\n",
        "else:\n",
        "    bands = dict()\n",
        "    result_dataset = None\n",
        "\n",
        "bidx = 0\n",
        "all_descriptions = []\n",
        "source_band = dict()\n",
        "\n",
        "def write_band_data(data, description, shortname=None):\n",
        "    global all_descriptions\n",
        "    global bidx\n",
        "    all_descriptions += [description]\n",
        "    if output_file_format == \"npz\":\n",
        "        bands[description] = data\n",
        "    else:\n",
        "        bidx += 1\n",
        "        result_dataset.GetRasterBand(bidx).SetDescription(description)\n",
        "        result_dataset.GetRasterBand(bidx).WriteArray(data)\n",
        "\n",
        "\n",
        "if sr60 is not None:\n",
        "    sr = np.concatenate((sr20, sr60), axis=2)\n",
        "    validated_sr_bands = validated_20m_bands + validated_60m_bands\n",
        "else:\n",
        "    sr = sr20\n",
        "    validated_sr_bands = validated_20m_bands\n",
        "\n",
        "if copy_original_bands:\n",
        "    out_dims = data10.shape[2] + sr.shape[2]\n",
        "else:\n",
        "    out_dims = sr.shape[2]\n",
        "\n",
        "sys.stdout.write(\"Writing\")\n",
        "result_dataset = driver.Create(output_file, data10.shape[1], data10.shape[0], out_dims, gdal.GDT_Float64)\n",
        "# Translate the image upper left corner. We multiply x10 to transform from pixel position in the 10m_band to meters.\n",
        "geot = list(ds10.GetGeoTransform())\n",
        "geot[0] += xmin * 10\n",
        "geot[3] -= ymin * 10\n",
        "result_dataset.SetGeoTransform(tuple(geot))\n",
        "result_dataset.SetProjection(ds10.GetProjection())\n",
        "\n",
        "if copy_original_bands:\n",
        "    sys.stdout.write(\" the original 10m bands and\")\n",
        "    # Write the original 10m bands\n",
        "    for bi, bn in enumerate(validated_10m_bands):\n",
        "        write_band_data(data10[:, :, bi], validated_descriptions[bn])\n",
        "print(\" the super-resolved bands in %s\" % output_file)\n",
        "for bi, bn in enumerate(validated_sr_bands):\n",
        "    write_band_data(sr[:, :, bi], \"SR\" + validated_descriptions[bn], \"SR\" + bn)\n",
        "\n",
        "\n",
        "for desc in all_descriptions:\n",
        "    print(desc)\n",
        "\n",
        "if output_file_format == \"npz\":\n",
        "    np.savez(output_file, bands=bands)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YQJHxOeixO3R",
        "colab_type": "text"
      },
      "source": [
        "Now we close the dataset and finish the writing process."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wt4sKB3UkQe9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "result_dataset=None"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B3HIgcUs2iyN",
        "colab_type": "text"
      },
      "source": [
        "The result can now be downloaded from the drive. In this example, it is around 200MB in size. This is because although we predicted all possible bands, and also included the original 10m bands, we only worked on a subset of the data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vECfcoMjviaq",
        "colab_type": "text"
      },
      "source": [
        "### Chapter 5: Visualisation\n",
        "\n",
        "After doing all that hard work, we now take a look at the result. We use the rasterio package to read in the file we just created.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y5A_CgfAvkI9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import rasterio\n",
        "from rasterio.plot import show\n",
        "from rasterio.plot import show_hist"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1rqBv76kyM9A",
        "colab_type": "text"
      },
      "source": [
        "#### Data\n",
        "\n",
        "We check the file and see that it has 12 bands. If we had decided not to copy the original bands, or not to predict the 60m bands, we would have fewer."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vRxP3idLv8pS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "src = rasterio.open(output_file)\n",
        "src.count"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TSfTBjU4yeVo",
        "colab_type": "text"
      },
      "source": [
        "We can also see that the file has the same dimension as the subset we specified at the beginning of Chapter 2. At least, it should have!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s7yIN9NkybCw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "src.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bzqeO9KKyyms",
        "colab_type": "text"
      },
      "source": [
        "####Histogram\n",
        "\n",
        "We check a histogram of our create image and see that most DN for all bands are between 0 and 5000."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rE64MJQhvu2t",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "show_hist(src, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
        "      histtype='stepfilled', title=\"Histogram\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oJDR6sVAHhg1",
        "colab_type": "text"
      },
      "source": [
        "#### True Color Images\n",
        "\n",
        "First, we enjoy a true color composite of the original bands, to make sure that it worked out and to have an overview of the study area.\n",
        "\n",
        "For doing this, we first use a little helper function to normalize the data into a range of 0-10"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FD-EvJzT7-9i",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "b2_red = src.read(1)\n",
        "b3_green = src.read(2)\n",
        "b4_blue = src.read(3)\n",
        "\n",
        "# Function to normalize the grid values\n",
        "def normalize(array):\n",
        "    \"\"\"Normalizes numpy arrays into scale 0.0 - 1.0\"\"\"\n",
        "    array_min, array_max = array.min(), array.max()\n",
        "    return (10*(array - array_min)/(array_max - array_min))\n",
        "\n",
        "    \n",
        "# Normalize the bands\n",
        "b2_redn = normalize(b2_red)\n",
        "b3_greenn = normalize(b3_green)\n",
        "b4_bluen = normalize(b4_blue)\n",
        "\n",
        "print(\"Normalized bands: \\n\")\n",
        "\n",
        "print(b2_redn.min(), '-', b2_redn.max(), 'mean:', b2_redn.mean())\n",
        "print(b3_greenn.min(), '-', b3_greenn.max(), 'mean:', b3_greenn.mean())\n",
        "print(b4_bluen.min(), '-', b4_bluen.max(), 'mean:', b4_bluen.mean())\n",
        "\n",
        "# Create RGB natural color composite\n",
        "rgb = np.dstack((b2_redn, b3_greenn, b4_bluen))\n",
        "plt.title('True Color Composite of the entire study area')\n",
        "print(\"\\n\\n\")\n",
        "# Let's see how our color composite looks like\n",
        "plt.imshow(rgb)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P_aqHd_X46dy",
        "colab_type": "text"
      },
      "source": [
        "#### Superresolved Bands\n",
        "\n",
        "What we are really interested in are the superresolved bands. After applying a little histogram stretch and normalizing to 0-255, which aids visualisation, we plot a few bands side by side:\n",
        "\n",
        "B8, B8a and B9 side by side.\n",
        "\n",
        "*   **B8**, an original 10m band\n",
        "*   **B8a**, originally 20m, now superresolved to 10m\n",
        "*   **B9**, originally 60m, now superresolved to 10m\n",
        "\n",
        "As all of these bands are in the NIR range of around 900nm wavelength, they should be roughly comparable.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i5B23bYG8SHq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "b8_nir = src.read(4)\n",
        "b8a_nir = src.read(8)\n",
        "b9_wv = src.read(12)\n",
        "\n",
        "# Normalize the bands\n",
        "\n",
        "def stretch(a, lower_thresh, upper_thresh):\n",
        "    r = 255.0/(upper_thresh-lower_thresh+2) # unit of stretching\n",
        "    out = np.round(r*(a-lower_thresh+1)).astype(a.dtype) # stretched values\n",
        "    out[a<lower_thresh] = 0\n",
        "    out[a>upper_thresh] = 255\n",
        "    return out\n",
        "\n",
        "b8_nirn = stretch(b8_nir,np.quantile(b8_nir,0.1),np.quantile(b8_nir,0.9))\n",
        "b8a_nirn = stretch(b8a_nir,np.quantile(b8a_nir,0.1),np.quantile(b8a_nir,0.9))\n",
        "b9_wvn = stretch(b9_wv,np.quantile(b9_wv,0.1),np.quantile(b9_wv,0.9))\n",
        "\n",
        "print(\"Stretched bands\")\n",
        "print(b8_nirn.min(), '-', b8_nirn.max(), 'mean:', b8_nirn.mean())\n",
        "print(b8a_nirn.min(), '-', b8a_nirn.max(), 'mean:', b8a_nirn.mean())\n",
        "print(b9_wvn.min(), '-', b9_wvn.max(), 'mean:', b9_wvn.mean())\n",
        "\n",
        "\n",
        "from matplotlib import pyplot\n",
        "fig, (axr, axg, axb) = pyplot.subplots(1,3, figsize=(21,7))\n",
        "show(b8_nirn, ax=axr, cmap='Reds', title='B8 original: 833 nm (10m)')\n",
        "show(b8a_nirn, ax=axg, cmap='Greens', title='B8a superresolved: 865 nm (20m -> 10m)')\n",
        "show(b9_wvn, ax=axb, cmap='Blues', title='B9 superresolved: 945 nm (60m -> 10m)')\n",
        "pyplot.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ymGw0VZnK__F",
        "colab_type": "text"
      },
      "source": [
        "At this range, we see that the results are visually pleasing, but at this scale we see little difference. We now also try the same with a smaller subset, effectively zooming in."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dMxtrz_EI9mm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "b8_nirn_sel= b8_nirn[1000:1200,1000:1200]\n",
        "b8a_nirn_sel = b8a_nirn[1000:1200,1000:1200]\n",
        "b9_wvn_sel = b9_wvn[1000:1200,1000:1200]\n",
        "\n",
        "fig, (axr, axg, axb) = pyplot.subplots(1,3, figsize=(21,7))\n",
        "show(b8_nirn_sel, ax=axr, cmap='Oranges', title='B8 original: 833 nm (10m)')\n",
        "show(b8a_nirn_sel, ax=axg, cmap='Oranges', title='B8a superresolved: 865 nm (20m -> 10m)')\n",
        "show(b9_wvn_sel, ax=axb, cmap='Oranges', title='B9 superresolved: 945 nm (60m -> 10m)')\n",
        "pyplot.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oPfHVqzuLQTM",
        "colab_type": "text"
      },
      "source": [
        "We see that the superresolved bands do not look too bad compared to the original 10m band! Nice!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W34LhH7XzS4o",
        "colab_type": "text"
      },
      "source": [
        "### Chapter 6: Streamlining\n",
        "\n",
        "So far we have executed the code step-by step. This made it easy to understand the code, but to perform superresolution on large numbers of files, it is not practical.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aB6_kKkMYk2i",
        "colab_type": "text"
      },
      "source": [
        "#### Wrapping the process into a function\n",
        "\n",
        " To speed things up for the future, we put the essential parts of the superresolution procedure into a function `DSen2_sl`. We do not include things like query options, and intermediate outputs.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xAB3ho0ozdCp",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "#@title\n",
        "import os\n",
        "def DSen2_sl( data_file,output_file,roi_lon_lat=\"\",roi_x_y=\"\",run_60=True,copy_original_bands=True,output_file_format=\"GTiff\"):\n",
        "  \"function_docstring\"\n",
        "  print(\"=========STARTING DSen2==========\")\n",
        "  print(\"data file:\"+data_file)\n",
        "  print(\"output file:\"+output_file)\n",
        "  print(\"roi_lon_lat:\"+roi_lon_lat)\n",
        "  print(\"roi_x_y:\"+roi_x_y)\n",
        "  print(\"run_60:\"+str(run_60))\n",
        "  print(\"output_file_format:\"+str(output_file_format))\n",
        "  print(\"copy_original_bands file:\"+str(copy_original_bands))\n",
        "\n",
        "  if run_60:\n",
        "      select_bands = 'B1,B2,B3,B4,B5,B6,B7,B8,B8A,B9,B11,B12'\n",
        "  else:\n",
        "      select_bands = 'B2,B3,B4,B5,B6,B7,B8,B8A,B11,B12'\n",
        "\n",
        "  # convert comma separated band list into a list\n",
        "  select_bands = [x for x in re.split(',', select_bands)]\n",
        "\n",
        "\n",
        "  if roi_lon_lat:\n",
        "      roi_lon1, roi_lat1, roi_lon2, roi_lat2 = [float(x) for x in re.split(',', roi_lon_lat)]\n",
        "  else:\n",
        "      roi_lon1, roi_lat1, roi_lon2, roi_lat2 = -180, -90, 180, 90\n",
        "\n",
        "  if roi_x_y:\n",
        "      roi_x1, roi_y1, roi_x2, roi_y2 = [float(x) for x in re.split(',', roi_x_y)]\n",
        "\n",
        "  raster = gdal.Open(data_file)\n",
        "  datasets = raster.GetSubDatasets();\n",
        "  tenMsets = []\n",
        "  twentyMsets = []\n",
        "  sixtyMsets = []\n",
        "  unknownMsets = []\n",
        "  for (dsname, dsdesc) in datasets:\n",
        "      if '10m resolution' in dsdesc:\n",
        "          tenMsets += [(dsname, dsdesc)]\n",
        "      elif '20m resolution' in dsdesc:\n",
        "          twentyMsets += [(dsname, dsdesc)]\n",
        "      elif '60m resolution' in dsdesc:\n",
        "          sixtyMsets += [(dsname, dsdesc)]\n",
        "      else:\n",
        "          unknownMsets += [(dsname, dsdesc)]\n",
        "\n",
        "\n",
        "\n",
        "            # case where we have several UTM in the data set\n",
        "  # => select the one with maximal coverage of the study zone\n",
        "  utm_idx = 0\n",
        "  utm = select_UTM\n",
        "  all_utms = defaultdict(int)\n",
        "  xmin, ymin, xmax, ymax = 0, 0, 0, 0\n",
        "  largest_area = -1\n",
        "  # process even if there is only one 10m set, in order to get roi -> pixels\n",
        "  for (tmidx, (dsname, dsdesc)) in enumerate(tenMsets + unknownMsets):\n",
        "      ds = gdal.Open(dsname)\n",
        "      if roi_x_y:\n",
        "          tmxmin = max(min(roi_x1, roi_x2, ds.RasterXSize - 1), 0)\n",
        "          tmxmax = min(max(roi_x1, roi_x2, 0), ds.RasterXSize - 1)\n",
        "          tmymin = max(min(roi_y1, roi_y2, ds.RasterYSize - 1), 0)\n",
        "          tmymax = min(max(roi_y1, roi_y2, 0), ds.RasterYSize - 1)\n",
        "          # enlarge to the nearest 60 pixel boundary for the super-resolution\n",
        "          tmxmin = int(tmxmin / 6) * 6\n",
        "          tmxmax = int((tmxmax + 1) / 6) * 6 - 1\n",
        "          tmymin = int(tmymin / 6) * 6\n",
        "          tmymax = int((tmymax + 1) / 6) * 6 - 1\n",
        "      elif not roi_lon_lat:\n",
        "          tmxmin = 0\n",
        "          tmxmax = ds.RasterXSize - 1\n",
        "          tmymin = 0\n",
        "          tmymax = ds.RasterYSize - 1\n",
        "      else:\n",
        "          xoff, a, b, yoff, d, e = ds.GetGeoTransform()\n",
        "          srs = osr.SpatialReference()\n",
        "          srs.ImportFromWkt(ds.GetProjection())\n",
        "          srsLatLon = osr.SpatialReference()\n",
        "          srsLatLon.SetWellKnownGeogCS(\"WGS84\");\n",
        "          ct = osr.CoordinateTransformation(srsLatLon, srs)\n",
        "\n",
        "\n",
        "          def to_xy(lon, lat):\n",
        "              (xp, yp, h) = ct.TransformPoint(lon, lat, 0.)\n",
        "              xp -= xoff\n",
        "              yp -= yoff\n",
        "              # matrix inversion\n",
        "              det_inv = 1. / (a * e - d * b)\n",
        "              x = (e * xp - b * yp) * det_inv\n",
        "              y = (-d * xp + a * yp) * det_inv\n",
        "              return (int(x), int(y))\n",
        "\n",
        "\n",
        "          x1, y1 = to_xy(roi_lon1, roi_lat1)\n",
        "          x2, y2 = to_xy(roi_lon2, roi_lat2)\n",
        "          tmxmin = max(min(x1, x2, ds.RasterXSize - 1), 0)\n",
        "          tmxmax = min(max(x1, x2, 0), ds.RasterXSize - 1)\n",
        "          tmymin = max(min(y1, y2, ds.RasterYSize - 1), 0)\n",
        "          tmymax = min(max(y1, y2, 0), ds.RasterYSize - 1)\n",
        "          # enlarge to the nearest 60 pixel boundary for the super-resolution\n",
        "          tmxmin = int(tmxmin / 6) * 6\n",
        "          tmxmax = int((tmxmax + 1) / 6) * 6 - 1\n",
        "          tmymin = int(tmymin / 6) * 6\n",
        "          tmymax = int((tmymax + 1) / 6) * 6 - 1\n",
        "      area = (tmxmax - tmxmin + 1) * (tmymax - tmymin + 1)\n",
        "      current_utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "      if area > all_utms[current_utm]:\n",
        "          all_utms[current_utm] = area\n",
        "      if current_utm == select_UTM:\n",
        "          xmin, ymin, xmax, ymax = tmxmin, tmymin, tmxmax, tmymax\n",
        "          utm_idx = tmidx\n",
        "          utm = current_utm\n",
        "          break\n",
        "      if area > largest_area:\n",
        "          xmin, ymin, xmax, ymax = tmxmin, tmymin, tmxmax, tmymax\n",
        "          largest_area = area\n",
        "          utm_idx = tmidx\n",
        "          utm = dsdesc[dsdesc.find(\"UTM\"):]\n",
        "\n",
        "  selected_10m_data_set = None\n",
        "  if not tenMsets:\n",
        "      selected_10m_data_set = unknownMsets[0]\n",
        "  else:\n",
        "      selected_10m_data_set = tenMsets[utm_idx]\n",
        "  selected_20m_data_set = None\n",
        "  for (dsname, dsdesc) in enumerate(twentyMsets):\n",
        "      if utm in dsdesc:\n",
        "          selected_20m_data_set = (dsname, dsdesc)\n",
        "  # if not found, assume the listing is in the same order\n",
        "  # => OK if only one set\n",
        "  if not selected_20m_data_set: selected_20m_data_set = twentyMsets[utm_idx]\n",
        "  selected_60m_data_set = None\n",
        "  for (dsname, dsdesc) in enumerate(sixtyMsets):\n",
        "      if utm in dsdesc:\n",
        "          selected_60m_data_set = (dsname, dsdesc)\n",
        "  if not selected_60m_data_set: selected_60m_data_set = sixtyMsets[utm_idx]\n",
        "\n",
        "  ds10 = gdal.Open(selected_10m_data_set[0])\n",
        "  ds20 = gdal.Open(selected_20m_data_set[0])\n",
        "  ds60 = gdal.Open(selected_60m_data_set[0])\n",
        "\n",
        "\n",
        "\n",
        "  def validate_description(description):\n",
        "      m = re.match(\"(.*?), central wavelength (\\d+) nm\", description)\n",
        "      if m:\n",
        "          return m.group(1) + \" (\" + m.group(2) + \" nm)\"\n",
        "      # Some HDR restrictions... ENVI band names should not include commas\n",
        "      if output_file_format == 'ENVI' and ',' in description:\n",
        "          pos = description.find(',')\n",
        "          return description[:pos] + description[(pos + 1):]\n",
        "      return description\n",
        "\n",
        "  def get_band_short_name(description):\n",
        "      if ',' in description:\n",
        "          return description[:description.find(',')]\n",
        "      if ' ' in description:\n",
        "          return description[:description.find(' ')]\n",
        "      return description[:3]\n",
        "\n",
        "\n",
        "  validated_10m_bands = []\n",
        "  validated_10m_indices = []\n",
        "  validated_20m_bands = []\n",
        "  validated_20m_indices = []\n",
        "  validated_60m_bands = []\n",
        "  validated_60m_indices = []\n",
        "  validated_descriptions = defaultdict(str)\n",
        "\n",
        "  for b in range(0, ds10.RasterCount):\n",
        "      desc = validate_description(ds10.GetRasterBand(b + 1).GetDescription())\n",
        "      shortname = get_band_short_name(desc)\n",
        "      if shortname in select_bands:\n",
        "          sys.stdout.write(\" \" + shortname)\n",
        "          select_bands.remove(shortname)\n",
        "          validated_10m_bands += [shortname]\n",
        "          validated_10m_indices += [b]\n",
        "          validated_descriptions[shortname] = desc\n",
        "  for b in range(0, ds20.RasterCount):\n",
        "      desc = validate_description(ds20.GetRasterBand(b + 1).GetDescription())\n",
        "      shortname = get_band_short_name(desc)\n",
        "      if shortname in select_bands:\n",
        "          sys.stdout.write(\" \" + shortname)\n",
        "          select_bands.remove(shortname)\n",
        "          validated_20m_bands += [shortname]\n",
        "          validated_20m_indices += [b]\n",
        "          validated_descriptions[shortname] = desc\n",
        "  for b in range(0, ds60.RasterCount):\n",
        "      desc = validate_description(ds60.GetRasterBand(b + 1).GetDescription())\n",
        "      shortname = get_band_short_name(desc)\n",
        "      if shortname in select_bands:\n",
        "          sys.stdout.write(\" \" + shortname)\n",
        "          select_bands.remove(shortname)\n",
        "          validated_60m_bands += [shortname]\n",
        "          validated_60m_indices += [b]\n",
        "          validated_descriptions[shortname] = desc\n",
        "  sys.stdout.write(\"\\n\")\n",
        "\n",
        "  # Some HDR restrictions... ENVI file name should be the .bin, not the .hdr\n",
        "  if output_file_format == 'ENVI' and (output_file[-4:] == '.hdr' or output_file[-4:] == '.HDR'):\n",
        "      output_file = output_file[:-4] + '.bin'\n",
        "\n",
        "\n",
        "  if validated_10m_indices:\n",
        "      data10 = np.rollaxis(\n",
        "          ds10.ReadAsArray(xoff=xmin, yoff=ymin, xsize=xmax - xmin + 1, ysize=ymax - ymin + 1, buf_xsize=xmax - xmin + 1,\n",
        "                          buf_ysize=ymax - ymin + 1), 0, 3)[:, :, validated_10m_indices]\n",
        "\n",
        "  if validated_20m_indices:\n",
        "      data20 = np.rollaxis(\n",
        "          ds20.ReadAsArray(xoff=xmin // 2, yoff=ymin // 2, xsize=(xmax - xmin + 1) // 2, ysize=(ymax - ymin + 1) // 2,\n",
        "                          buf_xsize=(xmax - xmin + 1) // 2, buf_ysize=(ymax - ymin + 1) // 2), 0, 3)[:, :,\n",
        "              validated_20m_indices]\n",
        "\n",
        "  if validated_60m_indices:\n",
        "      data60 = np.rollaxis(\n",
        "          ds60.ReadAsArray(xoff=xmin // 6, yoff=ymin // 6, xsize=(xmax - xmin + 1) // 6, ysize=(ymax - ymin + 1) // 6,\n",
        "                          buf_xsize=(xmax - xmin + 1) // 6, buf_ysize=(ymax - ymin + 1) // 6), 0, 3)[:, :,\n",
        "              validated_60m_indices]\n",
        "\n",
        "  maindir = os.getcwd()\n",
        "  subdir = maindir+\"/DSen2/testing/\"\n",
        "  os.chdir(subdir)\n",
        "  ##CODE START\n",
        "  print(\"Predicting.\")\n",
        "  if validated_60m_bands and validated_20m_bands and validated_10m_bands:\n",
        "      sr60 = DSen2_60(data10, data20, data60, deep=False)\n",
        "  else:\n",
        "      sr60 = None\n",
        "\n",
        "  if validated_10m_bands and validated_20m_bands:\n",
        "      sr20 = DSen2_20(data10, data20, deep=False)\n",
        "  else:\n",
        "      sr20 = None\n",
        "  ##CODE END\n",
        "  os.chdir(maindir)\n",
        "\n",
        "\n",
        "\n",
        "  if output_file_format != \"npz\":\n",
        "      revert_to_npz = True\n",
        "      driver = gdal.GetDriverByName(output_file_format)\n",
        "      if driver:\n",
        "          metadata = driver.GetMetadata()\n",
        "          if gdal.DCAP_CREATE in metadata and metadata[gdal.DCAP_CREATE] == 'YES':\n",
        "              revert_to_npz = False\n",
        "      if revert_to_npz:\n",
        "          print(\"Gdal doesn't support creating %s files\" % output_file_format)\n",
        "          print(\"Writing to npz as a fallback\")\n",
        "          output_file_format = \"npz\"\n",
        "      bands = None\n",
        "  else:\n",
        "      bands = dict()\n",
        "      result_dataset = None\n",
        "  global bidx\n",
        "  bidx = 0\n",
        "  all_descriptions = []\n",
        "  source_band = dict()\n",
        "  def write_band_data(data, description, shortname=None):\n",
        "      global all_descriptions\n",
        "      global bidx\n",
        "      all_descriptions += [description]\n",
        "      if output_file_format == \"npz\":\n",
        "          bands[description] = data\n",
        "      else:\n",
        "          bidx += 1\n",
        "          result_dataset.GetRasterBand(bidx).SetDescription(description)\n",
        "          result_dataset.GetRasterBand(bidx).WriteArray(data)\n",
        "\n",
        "\n",
        "  if sr60 is not None:\n",
        "      sr = np.concatenate((sr20, sr60), axis=2)\n",
        "      validated_sr_bands = validated_20m_bands + validated_60m_bands\n",
        "  else:\n",
        "      sr = sr20\n",
        "      validated_sr_bands = validated_20m_bands\n",
        "\n",
        "  if copy_original_bands:\n",
        "      out_dims = data10.shape[2] + sr.shape[2]\n",
        "  else:\n",
        "      out_dims = sr.shape[2]\n",
        "\n",
        "  print(\"Writing to:\"+str(output_file))\n",
        "  result_dataset = driver.Create(output_file, data10.shape[1], data10.shape[0], out_dims, gdal.GDT_Float64)\n",
        "  # Translate the image upper left corner. We multiply x10 to transform from pixel position in the 10m_band to meters.\n",
        "  geot = list(ds10.GetGeoTransform())\n",
        "  geot[0] += xmin * 10\n",
        "  geot[3] -= ymin * 10\n",
        "  result_dataset.SetGeoTransform(tuple(geot))\n",
        "  result_dataset.SetProjection(ds10.GetProjection())\n",
        "  if copy_original_bands:\n",
        "      # Write the original 10m bands\n",
        "      for bi, bn in enumerate(validated_10m_bands):\n",
        "        write_band_data(data10[:, :, bi], validated_descriptions[bn])\n",
        "  for bi, bn in enumerate(validated_sr_bands):\n",
        "      write_band_data(sr[:, :, bi], \"SR\" + validated_descriptions[bn], \"SR\" + bn)\n",
        "\n",
        "\n",
        "  if output_file_format == \"npz\":\n",
        "      np.savez(output_file, bands=bands)\n",
        "\n",
        "  result_dataset=None\n",
        "  return(output_file)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZodU-arw82uc",
        "colab_type": "text"
      },
      "source": [
        "\n",
        " We still repeat some steps like the creation of a new model, which is not optimal. Still, it allows us to semi-automate the superresolution."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OdNagbYfYx2b",
        "colab_type": "text"
      },
      "source": [
        "#### Processing multiple files\n",
        "\n",
        "Now with this function, it is easy to loop over a list of previously downloaded files, applying the superresolution algorithm to all of them and writing their results to the drive."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9BBHVKYv-W4d",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "downloaded_files"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9khPfgndXfDo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for i in range(0, len(downloaded_files)):\n",
        "    DSen2_sl(downloaded_files[i],\"Outputs/DSen2_Output\"+str(i)+\".Tif\",roi_x_y=\"1,1,300,300\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E9KAKjvcZAKD",
        "colab_type": "text"
      },
      "source": [
        "That's it! We can now do this for as many files as we can get our hands on. When they are done, we just download them from our drive!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ceSF2rUu7EDn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4FEAO8uD7Ecn",
        "colab_type": "text"
      },
      "source": [
        "## Final Note: Emulating the publication\n",
        "\n",
        "In this script, we have made a good number of changes to the source code. Nevertheless, it is possible to closely follow the procedure as described in the publication. We just have to make sure a few settings match:\n",
        "\n",
        "In the Training section:\n",
        "\n",
        "* **NR_CROP** should be set to `8000`.\n",
        "* **SCALE** should be set to `2000`\n",
        "* **n_epochs** should be set to `8192`.\n",
        "\n",
        "In the Superresolution section:\n",
        "\n",
        "* **roi_lon_lat** should be an empty string `\"\"` to predict on whole images.\n",
        "* **roi_x_y** should be an empty string `\"\"` to predict on whole images.\n",
        "* **select_UTM** should also be an empty string to let the algorithm choose the optimal UTM zone.\n",
        "* **run_60** and **copy_original_bands** should be set to `True` to get a full stack of 10m resolution.\n",
        "\n",
        "The publication compares the performances of the deep DSen2 model and the very deep VDSen2 model. If we want to use the latter, we  have to do two things:\n",
        "\n",
        "* Acquire the weights for the VCSen2 model and place them in the `DSen2/models` directory.\n",
        "* Set the switch **deep** to `True`.\n",
        "\n",
        "Using these settings, the algorithm can be executed as described in the publication.\n",
        "\n",
        "To recreate the study itself, we can train from scratch and on the same Sentinel-2 images as the study, a [list](https://github.com/lanha/DSen2/blob/master/S2_tiles_training.txt) of which is available on the Github. Another [list](https://github.com/lanha/DSen2/blob/master/S2_tiles_testing.txt) is available for the images which were used for testing. Using these images it should be possible to closely recreate the study."
      ]
    }
  ]
}