Combined_UBurn_WorkBook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combined (Backscattering and Coherence) Analysis Notebook\n",
    "\n",
    "When initially deploying the notebook, you will be asked to choose the dataset (either Coherence or Backscatter) and then the case study area (Skye, Cairngorms or Peak District)\n",
    "\n",
    "Once these have been chosen, the third drop down menu will automatically be populated with a list of available polygons (there may be a brief delay). These are listed by a unique reference number, followed by the date on which the burn occured. \n",
    "\n",
    "The last dropdown menu, is the dataset version you want to investigate.\n",
    "\n",
    "**Note:** Upon initially processing the coherence dataset, the tiles were paired based on the slice number from Sentinel-1. However, upon analysing the results at the end of the project, it became apparent that this didn't always provide full coverage of the area due to slight staggering in the area they covered. Therefore, additional combinations would need to be processed to fill these gaps in some images / orbits. Therfore, some orbits were reprocessed for the Peak District National Park and Cairngorms but due to the way the slices lay, we already had a full dataset for the Isle of Skye so no v2 dataset was created. For Backscatter there is only a single dataset. \n",
    "\n",
    "**Note:** For the test dataset limited data has been provided. Skye v1 relative orbit 125, Cairngorms v2 relative orbit 30 and Peak District v2 relative orbit 154. The same orbits have been provided for the backscatter data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded configuration: {'sharedfolder': 'my_shared_data_folder'}.\n",
      "Datasets in /home/slavender/my_shared_data_folder/datasets\n"
     ]
    }
   ],
   "source": [
    "#Reload modules without shutting notebook\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.windows import Window\n",
    "import datetime\n",
    "import pandas as pd\n",
    "import geojson as gj\n",
    "import matplotlib.pyplot as plt \n",
    "from ipywidgets import Dropdown, interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "import yaml\n",
    "\n",
    "# Import upburn utils\n",
    "from utils.array_creation import generate_array, generate_array_backscatter, get_window\n",
    "from utils.get_configuration import get_config\n",
    "from utils.functions import order_files_by_date, get_data_dict, get_aoi_date_dict, get_polygon_list\n",
    "from utils.load_corine import subset_corine\n",
    "from utils.plot_jsons import plot_jsons\n",
    "from utils.extract_s1ard import extract_s1ard\n",
    "from utils.calc_era_api import calc_era_api\n",
    "\n",
    "current_burn = None\n",
    "home = os.path.expanduser(\"~\")\n",
    "utils_dir = os.path.join(os.getcwd(), \"utils\")\n",
    "with open(os.path.join(utils_dir, \"configuration.yml\")) as yaml_file:\n",
    "    config = yaml.safe_load(yaml_file)\n",
    "    print(\"Loaded configuration: {:}.\".format(config))\n",
    "    sharedfolder = os.path.expanduser(config.get(\"sharedfolder\"))\n",
    "shared_folder = os.path.join(home, sharedfolder)\n",
    "datasets = os.path.join(shared_folder, \"datasets\")\n",
    "print(\"Datasets in %s\"%(datasets))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose the type of data to look at\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cdb1829a90a4441c96f37860c7be9886",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Dataset:', options=('', 'Coherence', 'Backscatter'), value='')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose a Case study area\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "895f15463a28486384519bab46deaa70",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Case study:', options=('', 'skye', 'cairngorms', 'pdistrict'), value='')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose burn area\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "513693c062084543abb1d61a4b39d506",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Polygon:', options=(), value=None)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose dataset version\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a88e11aacee8469294f16150c4858bc9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Version:', options=('v2', 'v1'), value='v2')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a widget that allows us to choose our analysis type\n",
    "workbook_widget = Dropdown(options = [\"\", \"Coherence\", \"Backscatter\"], description = \"Dataset:\", value = \"\")\n",
    "def change_wbook(*args): \n",
    "    global workbook \n",
    "    workbook = workbook_widget.value\n",
    "    if workbook == \"\":\n",
    "        return\n",
    "\n",
    "workbook_widget.observe(change_wbook, 'value')\n",
    "print(\"Choose the type of data to look at\")\n",
    "display(workbook_widget)\n",
    "\n",
    "# Create a widget that can be used the choose the case study area\n",
    "cstudy_widget = Dropdown(options = [\"\", \"skye\", \"cairngorms\", \"pdistrict\"], description = \"Case study:\", value = \"\")\n",
    "def change_cstudy(*args):\n",
    "    global cstudy \n",
    "    global aoi\n",
    "    global poly_list\n",
    "    global poly_dict\n",
    "    global polygons\n",
    "    cstudy = cstudy_widget.value\n",
    "    length = 3 # Number of characters in the id number\n",
    "    # Set the aoi depending on the chosen cstudy\n",
    "    if cstudy == \"\":\n",
    "        poly_list = []\n",
    "        return\n",
    "    elif cstudy == \"skye\":\n",
    "        aoi = os.path.join(datasets, \"Skye_extent_OSGB36.geojson\")\n",
    "    elif cstudy == \"cairngorms\":\n",
    "        aoi = os.path.join(datasets, \"Cairngorms_extent_OSGB36-extended.geojson\")    \n",
    "    else:\n",
    "        aoi = os.path.join(datasets, \"PDistrict_extent_OSGB36.geojson\")\n",
    "        length = 5\n",
    "    # get a list of polygons and a dictionary linking the polygons to the burn date\n",
    "    poly_list = get_polygon_list(datasets, aoi)\n",
    "    date_dict = get_aoi_date_dict(cstudy, shared_folder)\n",
    "    \n",
    "    ###  Below needs updating when we get dates for burns ###    \n",
    "    fstart = os.path.basename(aoi).split(\"_\")[0] \n",
    "    searchstr = fstart + \"_burn_extent_*.geojson\"\n",
    "    polygons = glob.glob(os.path.join(datasets, searchstr)) # get the polygons relevant to the chosen cstudy \n",
    "    if len(polygons) == 0:\n",
    "        print(\"No polygons exist for {}\".format(searchstr))\n",
    "    poly_list = []\n",
    "    for poly in polygons:\n",
    "        fpoly = os.path.splitext(os.path.basename(poly))[0]\n",
    "        poly_list.append(fpoly.split(\"_\")[3])\n",
    "\n",
    "    # Convert string list to int list, sort then convert back\n",
    "    poly_int = [i for i in poly_list]\n",
    "    poly_int.sort()\n",
    "    poly_list = [str(i).zfill(3) for i in poly_int]\n",
    "    poly_list2 = [pol + \" (%s)\"%(date_dict[int(pol[:length])]) for pol in poly_list]\n",
    "\n",
    "    poly_dict = {pol_string : pol for pol, pol_string in zip(poly_list, poly_list2)}\n",
    "    b_widget.options = poly_dict.keys() # set the polygon widget options \n",
    "    del poly_int\n",
    "    \n",
    "cstudy_widget.observe(change_cstudy, 'value')\n",
    "print(\"Choose a Case study area\")\n",
    "display(cstudy_widget)\n",
    "\n",
    "poly_dict = {}\n",
    "\n",
    "# Create a widget that can be used the choose the polygon of interest\n",
    "b_widget = Dropdown(options = poly_dict.keys(), description = \"Polygon:\")\n",
    "def change_x(*args):\n",
    "    print(\"Set to {}\".format(b_widget.value))\n",
    "    \n",
    "b_widget.observe(change_x, 'value')\n",
    "print(\"Choose burn area\")\n",
    "display(b_widget)\n",
    "\n",
    "# Create a widget that can be used the choose the polygon of interest\n",
    "version_widget = Dropdown(options = [\"v2\", \"v1\"], description = \"Version:\")\n",
    "def change_version(*args):\n",
    "    print(\"Set to {}\".format(version_widget.value))\n",
    "    \n",
    "version_widget.observe(change_version, 'value')\n",
    "print(\"Choose dataset version\")\n",
    "display(version_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display polygons\n",
    "Once you have made the selections above, you are able to visualise the locations of each burn by generating an interactive map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Displaying 39 polygons\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d25042151c904b949686c24e82569588",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='background', options=('map', 'image'), value='map'), Output(layout…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Show all the case study polygons on a map\n",
    "if len(polygons) < 1:\n",
    "        raise SystemExit(\"Could not find any polygons for {}\".format(searchstr))\n",
    "else:\n",
    "    print(\"Displaying {} polygons\".format(len(polygons)))\n",
    "\n",
    "# Plot on map\n",
    "wdown = widgets.Dropdown(\n",
    "    options=['map', 'image'],\n",
    "    value='map',\n",
    "    description='Background:'\n",
    ")\n",
    "\n",
    "def display_map(background):\n",
    "    m = plot_jsons(background, polygons, cstudy)\n",
    "    display(m)\n",
    "\n",
    "# display the dropdown and map\n",
    "interactive_map = interactive(display_map, background=['map','image'])\n",
    "output = interactive_map.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull in additional data based on the chosen polygon\n",
    "The following cells set the configuration based on the selections above, and pull in data from ERA5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chosen polygon: Skye_burn_extent_006-01\n",
      "Non-burn polygon available\n"
     ]
    }
   ],
   "source": [
    "# Set polygon from drop-down above\n",
    "endstr = \".geojson\"\n",
    "fstart = os.path.basename(aoi).split(\"_\")[0]\n",
    "burn = fstart + \"_burn_extent_\" + poly_dict[b_widget.value]\n",
    "non_burn = fstart + \"_baseline_extent_\" + poly_dict[b_widget.value].split(\"-\")[0]\n",
    "if os.path.exists(os.path.join(datasets, burn + endstr)):\n",
    "    print(\"Chosen polygon: {}\".format(burn))\n",
    "    if os.path.exists(os.path.join(datasets, non_burn + endstr)):\n",
    "        print(\"Non-burn polygon available\")\n",
    "        bdata = True\n",
    "    else:\n",
    "        print(\"No non-burn data available\")\n",
    "        bdata = False\n",
    "else:\n",
    "    print(\"ERROR: chosen polygon {} not found\".format(burn))\n",
    "chosen_burn = b_widget.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input dataset from: /home/slavender/my_shared_data_folder/datasets\n",
      "Already downloaded: /home/slavender/my_shared_data_folder/datasets/era-downloaded/Skye_burn_extent_006-01-download.nc \n",
      "Loading: odict_keys(['longitude', 'latitude', 'time', 'u10', 'v10', 'lai_hv', 'lai_lv', 'snowc', 'sp', 'sro', 'ssrd', 'tp'])  \n",
      "Extracted data shape:  (730, 9)\n",
      "Input dataset from: /home/slavender/my_shared_data_folder/datasets\n",
      "Already downloaded: /home/slavender/my_shared_data_folder/datasets/era-downloaded/Skye_baseline_extent_006-download.nc \n",
      "Loading: odict_keys(['longitude', 'latitude', 'time', 'u10', 'v10', 'lai_hv', 'lai_lv', 'snowc', 'sp', 'sro', 'ssrd', 'tp'])  \n",
      "Extracted data shape:  (730, 9)\n"
     ]
    }
   ],
   "source": [
    "# Load ERA5Land data\n",
    "graphics = False # do not display graphs\n",
    "verbose = False # do not show extra print statements\n",
    "polygon = burn\n",
    "%run ./Pull-ERA5.ipynb $polygon $verbose $graphics $cstudy\n",
    "if bdata:  # run for non burn polygon:\n",
    "    polygon = non_burn\n",
    "    %run ./Pull-ERA5.ipynb $polygon $verbose $graphics $cstudy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded configuration: {'sharedfolder': 'my_shared_data_folder'}.\n",
      "Loaded configuration: {'sharedfolder': 'my_shared_data_folder'}.\n",
      "Found 111 potential files and 0 baselines for v2\n"
     ]
    }
   ],
   "source": [
    "if version_widget.value == \"v2\":\n",
    "    v2=True\n",
    "    ofiles, outfolder, baselines = get_config(cstudy, workbook, v2=v2)\n",
    "    if len(ofiles) == 0:\n",
    "        print(\"No v2 data, switching to v1\")\n",
    "        v2=False\n",
    "else:\n",
    "    v2=False\n",
    "    ofiles, outfolder, baselines = get_config(cstudy, workbook, v2=v2)\n",
    "    if len(ofiles) == 0:\n",
    "        print(\"No v1 data, switching to v2\")\n",
    "        v2=True\n",
    "ofiles, outfolder, baselines = get_config(cstudy, workbook, v2=v2)\n",
    "if len(ofiles) == 0:\n",
    "     raise SystemExit(\"No valid data found\")\n",
    "else:\n",
    "    print(\"Found %s potential files and %s baselines for %s\"%(len(ofiles), len(baselines), version_widget.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract the data of interest\n",
    "Once the above selection has been completed, running the cell below will produce an interactive dataframe, where you can choose the files you want to investigate based on a combination of date range, orbit number and or orbit direction. \n",
    "When you are happy with your selection you can the following two cells will load the files you have chosen into a 3D array and display them. These can then be investigated visually by scrolling between them / switching between polarisation.\n",
    "\n",
    "You should see either 2 or 4 images (depending if a non burn area was available in the close vicinity for comparison). The first 2 images will show the burn area side by side, with the right hand image showing an extended area around it and the edge of the burn area drawn in red. The left image will show just the extracted burn area, with the surrounding area \"blanked out\".\n",
    "If a non burn area is available, the next 2 images will show the same for that area, respectively. \n",
    "\n",
    "The raw images / arrays (without the border drawn on or the area outside \"blanked\") can be output in a variety of ways. \n",
    "* \"Save\" will save a projected version of the image you are current viewing. \n",
    "* \"Save all burn geotiffs\" will save a projected version of every image chosen into a zip folder. \n",
    "* \"Save burn array\" will save the whole array as a python pickle object that can be loaded and analysed elsewhere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 111 input files for v2\n",
      "Polygon being extracted 2906 3618 2276 2816\n",
      "Removed 3 of 111 files that don't have data for this polygon\n"
     ]
    }
   ],
   "source": [
    "# get a list of only VV files and only VH files\n",
    "only_VV = [file for file in ofiles if \"VV\" in file]\n",
    "only_VH = [file for file in ofiles if \"VH\" in file]\n",
    "bands = set([file.strip(\"_%s.tif\"%(cstudy))[-2:] for file in ofiles])\n",
    "polar_dict = {'VV' : 0, 'VH' : 1}\n",
    "pol, other_pol = \"VV\", \"VH\"\n",
    "if workbook == \"Coherence\": # determine if we have more VV or VH as we will use this for the dataframe\n",
    "    if len(only_VV) > len(only_VH):\n",
    "        most_files = only_VV\n",
    "    else:\n",
    "        most_files = only_VH\n",
    "        pol, other_pol = \"VH\", \"VV\" \n",
    "else:\n",
    "    most_files = ofiles\n",
    "print(\"Found %s input files for %s\"%(len(most_files), version_widget.value))\n",
    "\n",
    "# Setup data frame\n",
    "most_files = order_files_by_date(most_files, workbook)\n",
    "data_dict = get_data_dict(most_files, workbook, baselines)\n",
    "pd.set_option('display.max_rows', 500)\n",
    "df = pd.DataFrame.from_dict(data_dict)\n",
    "wanted_files = df['filename']\n",
    "unwanted_files = []\n",
    "xmin, xmax, ymin, ymax, pixelWidth, pixelHeight = get_window(os.path.join(datasets, burn + \".geojson\"), wanted_files[0], just_coords= True)\n",
    "if xmin < 0: # polygon overlapping edge of area of interest\n",
    "    xmin = 0\n",
    "    print(\"Warning: Polygon overlaps the edge of the processed area\")\n",
    "else:\n",
    "    print(\"Polygon being extracted %d %d %d %d\"%(xmin, xmax, ymin, ymax))\n",
    "\n",
    "# loop through the files and remove any that don't have data as well as loading in a matching file for the other polarisation\n",
    "for n, file in enumerate(wanted_files): \n",
    "    win = Window.from_slices((ymin, ymax), (xmin, xmax))\n",
    "    with rasterio.open(file) as src:\n",
    "        arr = src.read(window = win)\n",
    "    if workbook == \"Backscatter\":\n",
    "        if np.mean(arr) == 0:\n",
    "            unwanted_files.append(n)\n",
    "            continue\n",
    "    matching_file = file.replace(pol, other_pol)\n",
    "    if os.path.exists(matching_file):\n",
    "        with rasterio.open(matching_file) as src:\n",
    "            arr2 = src.read(window = win)\n",
    "    else:\n",
    "        arr2 = np.zeros(arr.shape)\n",
    "    both_arrays = np.concatenate([arr, arr2])\n",
    "    if np.mean(both_arrays) == 0:\n",
    "        unwanted_files.append(n)\n",
    "        continue\n",
    "        \n",
    "# Drop unwanted files from the dataframe\n",
    "unwanted_files = sorted(unwanted_files, reverse=True)\n",
    "print(\"Removed %s of %s files that don't have data for this polygon\"%(len(unwanted_files), len(wanted_files)))\n",
    "\n",
    "df = df.drop(unwanted_files).reset_index(drop=True) # remove any images which have no data for either polarisation\n",
    "if df.empty:\n",
    "    raise SystemExit('DataFrame is empty after filtering')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0041eabeeb7e413eb7d9dab1df339ddb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(DatePicker(value=Timestamp('2018-01-02 00:00:00'), description='start_date'), DatePicker…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.get_files_by_date(start_date, end_date, rel_orbit, orbit_dir)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chosen_orbit = 'All'\n",
    "chosen_direction = 'both'        \n",
    "# create an interactive dataframe for filtering\n",
    "def get_files_by_date(start_date, end_date, rel_orbit, orbit_dir):\n",
    "    global chosen_orbit\n",
    "    global chosen_direction\n",
    "    global time_filtered_df\n",
    "    # filter it by date\n",
    "    time_filtered_df = df[(df['date'] >= pd.to_datetime(start_date)) & (df['date'] <= pd.to_datetime(end_date))]\n",
    "    if rel_orbit != \"\":\n",
    "        # filter by relative orbit\n",
    "        time_filtered_df = time_filtered_df[time_filtered_df[\"rorbit\"] == rel_orbit]\n",
    "        chosen_orbit = \"All\" \n",
    "    if orbit_dir != \"\":\n",
    "        # filter by orbit direction\n",
    "        time_filtered_df = time_filtered_df[time_filtered_df[\"direction\"] == orbit_dir]\n",
    "        chosen_direction = \"Both\"\n",
    "    if len(set(time_filtered_df.direction.values)) == 1:\n",
    "        # if we only have one orbit direction we can set this now\n",
    "        orbit_dir = time_filtered_df.direction.values[0]\n",
    "    time_filtered_df_vis = time_filtered_df.copy(deep=True)\n",
    "    # set the chosen values\n",
    "    chosen_orbit = rel_orbit\n",
    "    chosen_direction = orbit_dir\n",
    "    # return the filtered df without filename and polarisation for visualisation\n",
    "    return time_filtered_df_vis.drop(['filename', 'polarisation'], axis = 1)\n",
    "\n",
    "orbit_opts = list(set(df['rorbit'].values))\n",
    "orbit_opts.append(\"\")\n",
    "orbit_dir_opts = list(set(df['direction'].values))\n",
    "orbit_dir_opts.append(\"\")\n",
    "interact(get_files_by_date,\n",
    "        start_date=widgets.DatePicker(value=pd.to_datetime(min(df['date']))),\n",
    "        end_date=widgets.DatePicker(value=pd.to_datetime(max(df['date']))),\n",
    "        rel_orbit=widgets.Dropdown(options=orbit_opts, value = \"\"),\n",
    "        orbit_dir=widgets.Dropdown(options=orbit_dir_opts, value = \"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn some of the column info into lists we can use during plotting\n",
    "wanted_files = time_filtered_df.filename.to_list()\n",
    "list_of_dates = time_filtered_df.date.to_list()\n",
    "list_of_directions = time_filtered_df.direction.to_list()\n",
    "# lets load the correct data\n",
    "if workbook == \"Coherence\":\n",
    "    list_of_baselines = time_filtered_df.perpbaseline.to_list()\n",
    "    burn_subarray, burn_mask, burn_transform = generate_array(wanted_files, bands, polar_dict, pol, other_pol, os.path.join(datasets, burn + \".geojson\"))\n",
    "    if bdata:\n",
    "        non_burn_subarray, non_burn_mask, non_burn_transform = generate_array(wanted_files, bands, polar_dict, pol, other_pol, os.path.join(datasets, polygon + \".geojson\"))\n",
    "elif workbook == \"Backscatter\":\n",
    "    burn_subarray, burn_mask, burn_transform = generate_array_backscatter(wanted_files, 2, os.path.join(datasets, burn + \".geojson\"))\n",
    "    if bdata:\n",
    "        non_burn_subarray, non_burn_mask, non_burn_transform = generate_array_backscatter(wanted_files, 2, os.path.join(datasets, polygon + \".geojson\"))\n",
    "\n",
    "if burn_subarray is None:\n",
    "    raise SystemExit(\"No valid {} data found intersecting {}\".format(workbook, os.path.splitext(os.path.basename(polygon_file))[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ab33922f0fd14646824b3b55fb1574f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='polarisation', options=('VV', 'VH'), value='VV'), IntSlider(value=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "461d3c286b524d4ea4a6c1ad0b3cf2d9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save', style=ButtonStyle()), Output()), _dom_classes=('widget-intera…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "39e960bcaf6d434fb2ea7d57fdc783ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save all burn geotiffs', style=ButtonStyle()), Output()), _dom_class…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b9e322ff52984f4b87521893bd0341b9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save burn array', style=ButtonStyle()), Output()), _dom_classes=('wi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2a7a8d3a474a46d9aff69a1aba368634",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='polarisation_non_burn', options=('VV', 'VH'), value='VV'), IntSlid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "83074dfb2e1343629d6e72cc7391b07b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save', style=ButtonStyle()), Output()), _dom_classes=('widget-intera…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a9676c85b99f4b60975c313ebeba174e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save all non-burn geotiffs', style=ButtonStyle()), Output()), _dom_c…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e072e92a32f34e619683f6f25110c21a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Button(description='Save non burn array', style=ButtonStyle()), Output()), _dom_classes=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.ticker as ticker \n",
    "\n",
    "def save_image():\n",
    "    # create the filename for saving\n",
    "    dst_filename = os.path.join(outfolder, \"%s_%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], date, pol_text, direct, \"burn_%s\"%(workbook)))\n",
    "    if not os.path.exists(outfolder):\n",
    "        os.mkdir(outfolder)\n",
    "    if os.path.exists(dst_filename):\n",
    "        return\n",
    "    driver = gdal.GetDriverByName( 'GTiff' )\n",
    "    dst_ds=driver.Create(dst_filename, curr_arr.shape[1],curr_arr.shape[0], 1, gdal.GDT_Float32) # is shape the wrong way around?\n",
    "    dst_ds.SetGeoTransform(burn_transform)\n",
    "    srs = osr.SpatialReference()\n",
    "    srs.ImportFromEPSG(27700)\n",
    "    dest_wkt = srs.ExportToWkt()\n",
    "    dst_ds.SetProjection (dest_wkt) # is this right?\n",
    "\n",
    "    dst_ds.GetRasterBand(1).WriteArray( curr_arr )\n",
    "    print(\"Saved to: {}\".format(dst_filename))\n",
    "\n",
    "\n",
    "@interact(img_num=(1,len(burn_subarray),1), polarisation=['VV','VH'])\n",
    "def plot_polygons(polarisation='VV', img_num=0):\n",
    "    # get max and min vals to normalise plots with\n",
    "    max_val= np.max(burn_subarray)\n",
    "    min_val = np.min(burn_subarray)\n",
    "    if workbook == \"Coherence\":\n",
    "        max_val = 1\n",
    "        min_val = 0\n",
    "    global date\n",
    "    global pol_text\n",
    "    global curr_arr\n",
    "    global direct\n",
    "    pol_text = polarisation\n",
    "    pol = polar_dict[polarisation]\n",
    "    date = list_of_dates[img_num-1]\n",
    "    if workbook == \"Coherence\":\n",
    "        baseline = list_of_baselines[img_num-1]\n",
    "    if not isinstance(date, str):\n",
    "        date = date.strftime(\"%Y%m%d\")\n",
    "    direct = list_of_directions[img_num-1]\n",
    "    # set the current array so the user can save it if needed \n",
    "    curr_arr = np.copy(burn_subarray[(img_num-1), pol, :, :])\n",
    "    curr_arr[burn_mask==0] = np.nan # set anything outside the polygon to np.nan\n",
    "    fig = plt.figure(figsize=((20,10)))\n",
    "    ax = fig.add_subplot(121)\n",
    "    ax.imshow(curr_arr, vmax=max_val, vmin=min_val) # show the chosen array \n",
    "    if workbook == \"Coherence\":\n",
    "        plt.title(\"Burn %s: %s direction %s perpbaseline %s\"%(workbook, date, direct, baseline))\n",
    "    else:\n",
    "        plt.title(\"Burn %s: %s direction %s\"%(workbook, date, direct))\n",
    "    # Do not plot ticks\n",
    "    plt.xticks([]) \n",
    "    plt.yticks([]) \n",
    "              \n",
    "    curr_arr2= np.copy(burn_subarray[(img_num-1), pol, :, ])\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    im = ax2.imshow(curr_arr2, vmax=max_val, vmin=min_val) # show the current array \n",
    "    # this time we will keep the data outside the polygon and just draw a contour around it\n",
    "    ax2.contour(burn_mask > 0, levels=[2], colors=['r'])\n",
    "    ax2.xaxis.set_major_locator(ticker.NullLocator())\n",
    "    ax2.yaxis.set_major_locator(ticker.NullLocator())\n",
    "    \n",
    "    fig.subplots_adjust(right=0.85)# add space for colour bar\n",
    "    cbar_ax = fig.add_axes([0.88, 0.25, 0.04, 0.5])\n",
    "    fig.colorbar(im, cax=cbar_ax)\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "savebutton = interactive(save_image, {'manual' : True, 'manual_name' : 'Save'})\n",
    "display(savebutton)\n",
    "\n",
    "def save_all_images():\n",
    "    \"\"\"\n",
    "    Saves all the chosen images individually and zips them for downloading\n",
    "    \"\"\"\n",
    "    fold_name = list_of_dates[0].strftime(\"%Y%m%d\") + \"_\" + list_of_dates[-1].strftime(\"%Y%m%d\") + \"_\" + burn +\"_\"+ \"orb%s\"%(chosen_orbit)+ \"_\" +\"dir%s\"%(chosen_direction)+ \"_\" + workbook\n",
    "    fold_name = os.path.join(outfolder, fold_name)\n",
    "    if not os.path.exists(fold_name):\n",
    "        os.mkdir(fold_name)\n",
    "    for n, subarr in enumerate(burn_subarray, start=0):\n",
    "        # loop through and save each image to the folder\n",
    "        for pol in polar_dict.keys():\n",
    "            wanted_arr = subarr[polar_dict[pol], :, :]\n",
    "            arrdate = list_of_dates[n]\n",
    "            if not isinstance(arrdate, str):\n",
    "                arrdate = arrdate.strftime(\"%Y%m%d\")\n",
    "            direct = list_of_directions[n]\n",
    "            dst_filename = os.path.join(fold_name, \"%s_%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], arrdate, pol, list_of_directions[n], \"burn_%s\"%(workbook)))\n",
    "            if os.path.exists(dst_filename):\n",
    "                continue\n",
    "            driver = gdal.GetDriverByName( 'GTiff' )\n",
    "            dst_ds=driver.Create(dst_filename, wanted_arr.shape[1], wanted_arr.shape[0], 1, gdal.GDT_Float32) # is shape the wrong way around?\n",
    "            dst_ds.SetGeoTransform(burn_transform)\n",
    "            srs = osr.SpatialReference()\n",
    "            srs.ImportFromEPSG(27700)\n",
    "            dest_wkt = srs.ExportToWkt()\n",
    "            dst_ds.SetProjection (dest_wkt) # is this right?\n",
    "            dst_ds.GetRasterBand(1).WriteArray( wanted_arr )\n",
    "    import shutil\n",
    "    shutil.make_archive(fold_name, 'zip', fold_name)\n",
    "    shutil.rmtree(fold_name) # remove the non zipped folder\n",
    "    \n",
    "    print(\"Saved to: {}\".format(fold_name + \".zip\"))\n",
    "\n",
    "saveallbutton = interactive(save_all_images, {'manual' : True, 'manual_name' : 'Save all burn geotiffs'})\n",
    "display(saveallbutton)\n",
    "\n",
    "def save_whole_array():\n",
    "    import pickle \n",
    "    print(burn)\n",
    "    if not isinstance(list_of_dates[0], str):\n",
    "        start = list_of_dates[0].strftime(\"%Y%m%d\")\n",
    "        end = list_of_dates[-1].strftime(\"%Y%m%d\")\n",
    "        print(start, end)\n",
    "    else: \n",
    "        start = list_of_dates[0]\n",
    "        end = list_of_dates[-1]\n",
    "    fname = burn +\"_\"+ \"orb%s\"%(chosen_orbit)+ \"_\" +\"dir%s\"%(chosen_direction)+ \"_\" + start + \"_\" + end +  \".p\"\n",
    "    pickle.dump(burn_subarray, open(os.path.join(outfolder, fname), \"wb\"))\n",
    "    print(\"Saved to: {}\".format(os.path.join(outfolder, fname)))\n",
    "    \n",
    "save_whole_burn = interactive(save_whole_array, {'manual' : True, 'manual_name' : 'Save burn array'})\n",
    "display(save_whole_burn)\n",
    "\n",
    "if bdata:\n",
    "    def save_image2():\n",
    "        dst_filename = os.path.join(outfolder, \"%s_%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], date2, pol_text2, direct, \"burn_%s\"%(workbook)))\n",
    "        if not os.path.exists(outfolder):\n",
    "            os.mkdir(outfolder)\n",
    "        if os.path.exists(dst_filename):\n",
    "            return\n",
    "        driver = gdal.GetDriverByName( 'GTiff' )\n",
    "        dst_ds=driver.Create(dst_filename, curr_arr3.shape[1],curr_arr3.shape[0], 1, gdal.GDT_Float32) # is shape the wrong way around?\n",
    "        dst_ds.SetGeoTransform(non_burn_transform)\n",
    "        srs = osr.SpatialReference()\n",
    "        srs.ImportFromEPSG(27700)\n",
    "        dest_wkt = srs.ExportToWkt()\n",
    "        dst_ds.SetProjection (dest_wkt) # is this right?\n",
    "\n",
    "        dst_ds.GetRasterBand(1).WriteArray( curr_arr3 )\n",
    "        print(\"Saved to: {}\".format(dst_filename))\n",
    "\n",
    "\n",
    "    @interact(img_num_non_burn=(1,len(non_burn_subarray),1), polarisation_non_burn=['VV','VH'])\n",
    "    def plot_polygons2(polarisation_non_burn='VV', img_num_non_burn=0):\n",
    "        max_val= np.max(burn_subarray)\n",
    "        min_val = np.min(burn_subarray)\n",
    "        if workbook == \"Coherence\":\n",
    "            max_val = 1\n",
    "            min_val = 0\n",
    "        global date2\n",
    "        global pol_text2\n",
    "        global curr_arr3\n",
    "        global direct2\n",
    "        pol_text2 = polarisation_non_burn\n",
    "        pol = polar_dict[polarisation_non_burn]\n",
    "        date2 = list_of_dates[img_num_non_burn-1]\n",
    "        if workbook == \"Coherence\":\n",
    "            baseline2 = list_of_baselines[img_num_non_burn-1]\n",
    "        if not isinstance(date2, str):\n",
    "            date2 = date2.strftime(\"%Y%m%d\")\n",
    "        direct2 = list_of_directions[img_num_non_burn-1]\n",
    "        curr_arr3 = np.copy(non_burn_subarray[(img_num_non_burn-1), pol, :, :])\n",
    "        curr_arr3[non_burn_mask==0] = np.nan # set anything outside the polygon to np.nan\n",
    "        fig2 = plt.figure(figsize=((20,10)))\n",
    "        ax3 = fig2.add_subplot(121)\n",
    "        ax3.imshow(curr_arr3, vmax=max_val, vmin=min_val)   \n",
    "        if workbook == \"Coherence\":\n",
    "            plt.title(\"Non-Burn %s: %s direction %s perpbaseline %s\"%(workbook, date, direct, baseline2))\n",
    "        else:\n",
    "            plt.title(\"Non-Burn %s: %s direction %s\"%(workbook, date, direct))\n",
    "        # Do not plot axes ticks\n",
    "        plt.xticks([]) \n",
    "        plt.yticks([]) \n",
    "\n",
    "        curr_arr4= np.copy(non_burn_subarray[(img_num_non_burn-1), pol, :, ])\n",
    "        ax4 = fig2.add_subplot(122)\n",
    "        im2 =ax4.imshow(curr_arr4, vmax=max_val, vmin=min_val)\n",
    "        ax4.contour(non_burn_mask > 0, levels=[2], colors=['r'])\n",
    "        ax4.xaxis.set_major_locator(ticker.NullLocator())\n",
    "        ax4.yaxis.set_major_locator(ticker.NullLocator())\n",
    "        # add space for colour bar\n",
    "        fig2.subplots_adjust(right=0.85)# add space for colour bar\n",
    "        cbar_ax = fig2.add_axes([0.88, 0.25, 0.04, 0.5])\n",
    "        fig2.colorbar(im2, cax=cbar_ax)\n",
    "        plt.show()\n",
    "\n",
    "    savebutton2 = interactive(save_image2, {'manual' : True, 'manual_name' : 'Save'})\n",
    "    display(savebutton2)\n",
    "\n",
    "    def save_all_images2():\n",
    "        fold_name = list_of_dates[0].strftime(\"%Y%m%d\") + \"_\" + list_of_dates[-1].strftime(\"%Y%m%d\") + \"_\" + non_burn +\"_\"+ \"orb%s\"%(chosen_orbit)+ \"_\" +\"dir%s\"%(chosen_direction) + \"_\" + workbook\n",
    "        fold_name = os.path.join(outfolder, fold_name)\n",
    "        if not os.path.exists(fold_name):\n",
    "            os.mkdir(fold_name)\n",
    "        for n, subarr in enumerate(non_burn_subarray, start=0):\n",
    "            for pol in polar_dict.keys():\n",
    "                wanted_arr = subarr[polar_dict[pol], :, :]\n",
    "                arrdate = list_of_dates[n]\n",
    "                if not isinstance(arrdate, str):\n",
    "                    arrdate = arrdate.strftime(\"%Y%m%d\")\n",
    "                direct = list_of_directions[n]\n",
    "                dst_filename = os.path.join(fold_name, \"%s_%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], arrdate, pol, list_of_directions[n], \"burn_%s\"%(workbook)))\n",
    "                if os.path.exists(dst_filename):\n",
    "                    continue\n",
    "                driver = gdal.GetDriverByName( 'GTiff' )\n",
    "                dst_ds=driver.Create(dst_filename, wanted_arr.shape[1], wanted_arr.shape[0], 1, gdal.GDT_Float32) # is shape the wrong way around?\n",
    "                dst_ds.SetGeoTransform(non_burn_transform)\n",
    "                srs = osr.SpatialReference()\n",
    "                srs.ImportFromEPSG(27700)\n",
    "                dest_wkt = srs.ExportToWkt()\n",
    "                dst_ds.SetProjection (dest_wkt) # is this right?\n",
    "                dst_ds.GetRasterBand(1).WriteArray( wanted_arr )\n",
    "        import shutil\n",
    "        shutil.make_archive(fold_name, 'zip', fold_name)\n",
    "        shutil.rmtree(fold_name)\n",
    "        print(\"Saved to: {}\".format(fold_name + \".zip\"))\n",
    "\n",
    "    saveallbutton2 = interactive(save_all_images2, {'manual' : True, 'manual_name' : 'Save all non-burn geotiffs'})\n",
    "    display(saveallbutton2)\n",
    "\n",
    "    def save_whole_array2():\n",
    "        import pickle \n",
    "        if not isinstance(list_of_dates[0], str):\n",
    "            start = list_of_dates[0].strftime(\"%Y%m%d\")\n",
    "            end = list_of_dates[-1].strftime(\"%Y%m%d\")\n",
    "            print(start, end)\n",
    "        else: \n",
    "            start = list_of_dates[0]\n",
    "            end = list_of_dates[-1]\n",
    "        fname = non_burn+\"_\"+ \"orb%s\"%(chosen_orbit)+ \"_\" +\"dir%s\"%(chosen_direction)+ \"_\" + start + \"_\" + end +  \".p\"\n",
    "        pickle.dump(non_burn_subarray, open(os.path.join(outfolder, fname), \"wb\"))\n",
    "        print(\"Saved to: {}\".format(os.path.join(outfolder, fname)))\n",
    "\n",
    "    save_whole_burn2 = interactive(save_whole_array2, {'manual' : True, 'manual_name' : 'Save non burn array'})\n",
    "    display(save_whole_burn2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import CORINE land cover\n",
    "Running the next cell will allow you to chose a specific land cover classification you wish to further analyse. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<open Collection '/home/slavender/my_shared_data_folder/datasets/u2018_clc2018_v2020_20u1_geoPackage/DATA/U2018_CLC2018_V2020_20u1.gpkg:U2018_CLC2018_V2020_20u1', mode 'r' at 0x7fae44e36cf8>\n",
      "Extracted 3 polygons\n",
      "Choose skye land cover type: \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4855c7d6593441a0a571354981c74e1f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Landcover types:', options=('All', 'Peat bogs', 'Transitional woodland-shrub', 'Conifero…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "current_burn = chosen_burn \n",
    "burn_df = time_filtered_df\n",
    "string_dates = [pd.to_datetime(d).strftime(\"%Y%m%d\") for d in list_of_dates]\n",
    "burn_df['date'] = string_dates\n",
    "\n",
    "# Filter the S1ARD data\n",
    "filtered_burndf, filtered_burnarray, filtered_dates_burn = time_filtered_df, burn_subarray, time_filtered_df.date.to_list()\n",
    "if bdata:\n",
    "    filtered_basedf, filtered_basearray, filtered_dates_base = time_filtered_df, non_burn_subarray, time_filtered_df.date.to_list()\n",
    "#print(filtered_burndf, filtered_burnarray.shape, filtered_dates_burn)\n",
    "\n",
    "# Load CORINE land cover for polygon burn area subset\n",
    "feature_list = ['All']\n",
    "cpolygons, legend_df = subset_corine(xymin[0], xymin[1], xymax[0], xymax[1], pixelWidth, pixelHeight, datasets, verb = verbose)\n",
    "if cpolygons == None:\n",
    "    print(\"CORINE data is not available, see utils/load_corine for details of what is expected\")\n",
    "    feature_list = ['Whole-Polygon']\n",
    "    # Extract burn area using original polygon\n",
    "    mean_s1burn, std_s1burn, pnames, features = extract_s1ard(cpolygons, filtered_burnarray, verb = verbose)\n",
    "    if bdata:\n",
    "        mean_s1base, std_s1base, pnames, features = extract_s1ard(cpolygons, filtered_basearray, verb = verbose)\n",
    "else:\n",
    "    # Extract subsets of burn area using CORINE polygons\n",
    "    mean_s1burn, std_s1burn, pnames, features = extract_s1ard(cpolygons, filtered_burnarray, verb = verbose)\n",
    "    if bdata:\n",
    "        mean_s1base, std_s1base, pnames, features = extract_s1ard(cpolygons, filtered_basearray, verb = verbose)\n",
    "\n",
    "    for feature in features:\n",
    "        feature_list.append(feature)\n",
    "\n",
    "c_widget = Dropdown(options = feature_list, description = \"Landcover types:\", value = feature_list[0])\n",
    "def change_x(*args):\n",
    "    print(\"Set to {}\".format(c_widget.value))\n",
    "\n",
    "c_widget.observe(change_x, 'value')\n",
    "print(\"Choose {} land cover type: \".format(cstudy))\n",
    "display(c_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate precipitation index and generate plots\n",
    "\n",
    "Finally, the last 3 cells will determine precipitation, create plots and then finallyt display plots that can be analysed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "108 S1 dates: [Timestamp('2018-01-02 00:00:00'), Timestamp('2018-01-05 00:00:00'), Timestamp('2018-01-08 00:00:00'), Timestamp('2018-01-11 00:00:00'), Timestamp('2018-01-14 00:00:00')] ...\n",
      "S1 Backscatter array (3, 6, 108)\n",
      "   Burn Peat bogs max VV -8.884 max VH -16.842 VHrVV 2.230 RFDI -0.195 NVHI 0.674 NVVI 0.402\n",
      "  Non-burn Peat bogs max VV -8.946 max VH -16.655 VHrVV 1.704 RFDI -0.191 NVHI 0.675 NVVI 0.404\n",
      "   Burn Transitional woodland-shrub max VV -6.597 max VH -14.389 VHrVV 2.521 RFDI -0.213 NVHI 0.699 NVVI 0.393\n",
      "  Non-burn Transitional woodland-shrub max VV -9.657 max VH -15.838 VHrVV 1.955 RFDI -0.181 NVHI 0.658 NVVI 0.409\n",
      "   Burn Coniferous forest max VV -6.786 max VH -15.998 VHrVV 2.726 RFDI -0.185 NVHI 0.718 NVVI 0.408\n",
      "  Non-burn Coniferous forest max VV -9.363 max VH -15.785 VHrVV 1.955 RFDI -0.180 NVHI 0.658 NVVI 0.410\n",
      "ERA5 Precip: 108 values from 20180102 to 20190614 max 0.055\n"
     ]
    }
   ],
   "source": [
    "subpolygons = {}\n",
    "if c_widget.value == \"Whole-Polygon\":\n",
    "    print(\"Using the whole polygon as CORINE data not available\")\n",
    "elif c_widget.value != \"All\":\n",
    "    for feat, p in cpolygons.items():\n",
    "        if feat == c_widget.value:  # Matching names\n",
    "            subpolygons[feat] = p\n",
    "                    \n",
    "    # Run subset again for just the chosen feature\n",
    "    mean_s1burn, std_s1burn, pnames, features = extract_s1ard(subpolygons, filtered_burnarray, verb = verbose)\n",
    "    if bdata:\n",
    "        mean_s1base, std_s1base, pnames, features = extract_s1ard(subpolygons, filtered_basearray, verb = verbose)\n",
    "\n",
    "# Extract ERA5 data for matching dates and calculate Antecedent Precipitation Index\n",
    "filtered_dates = burn_df['date'].values\n",
    "mean_edata, std_edata = calc_era_api(dt_num, mean_parameters, filtered_burnarray, filtered_dates, verb = True)\n",
    "filtered_edata = filtered_dates\n",
    "eraname = 'Antecedent Precipitation Index'\n",
    "\n",
    "# Print initial arrays\n",
    "print(\"{} S1 dates: {} ...\".format(len(list_of_dates),list_of_dates[0:5]))\n",
    "print(\"S1 {} array {}\".format(workbook, mean_s1burn.shape))\n",
    "for fnum in range(len(features)):\n",
    "    print(\"   Burn {} max VV {:.3f} max VH {:.3f} VHrVV {:.3f} RFDI {:.3f} NVHI {:.3f} NVVI {:.3f}\"\n",
    "        .format(features[fnum],np.nanmax(mean_s1burn[fnum,0,:]),np.nanmax(mean_s1burn[fnum,1,:]),np.nanmax(mean_s1burn[fnum,2,:]),\n",
    "        np.nanmax(mean_s1burn[fnum,3,:]),np.nanmax(mean_s1burn[fnum,4,:]),np.nanmax(mean_s1burn[fnum,5,:])))\n",
    "    if bdata:\n",
    "        print(\"  Non-burn {} max VV {:.3f} max VH {:.3f} VHrVV {:.3f} RFDI {:.3f} NVHI {:.3f} NVVI {:.3f}\"\n",
    "            .format(features[fnum],np.nanmax(mean_s1base[fnum,0,:]),np.nanmax(mean_s1base[fnum,1,:]),np.nanmax(mean_s1base[fnum,2,:]),\n",
    "            np.nanmax(mean_s1base[fnum,3,:]),np.nanmax(mean_s1base[fnum,4,:]),np.nanmax(mean_s1base[fnum,5,:])))\n",
    "print(\"ERA5 Precip: {} values from {} to {} max {:.3f}\".format(len(mean_edata), filtered_edata[0], filtered_edata[-1], np.nanmax(mean_edata)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering: between 20180102 and 20190614 for rorbit  and direction \n",
      "Running for 3 features with 6 products\n",
      "Generated: /home/slavender/notebooks/output_folder/plot-Backscatter-era-non-burn.png\n",
      "Running for 3 features with 6 products\n",
      "Generated: /home/slavender/notebooks/output_folder/plot-Backscatter-era-burn.png\n"
     ]
    }
   ],
   "source": [
    "burn_start = current_burn.split(\" (\")[1][:-1]\n",
    "burn_start = datetime.strptime(burn_start, '%d/%m/%Y')\n",
    "# Create plots for ERA5 data in comparison to Sentinel-1 ARD dates for each CORINE polygon\n",
    "import utils.plot_data as plotd\n",
    "plt.set_loglevel(\"critical\")\n",
    "\n",
    "if not os.path.exists(outfolder):\n",
    "    os.mkdir(outfolder)\n",
    "\n",
    "print(\"Filtering: between {} and {} for rorbit {} and direction {}\".format(filtered_dates[0], filtered_dates[-1], chosen_orbit, chosen_direction)) \n",
    "if bdata:\n",
    "    outbase = os.path.join(outfolder,\"plot-\" + workbook + \"-era-non-burn.png\")\n",
    "    plotd.plot_lines(filtered_dates, mean_s1base, std_s1base, filtered_edata, mean_edata, std_edata, features, pnames, eraname, outbase, burn_start, graphics = False)\n",
    "    print(\"Generated: {}\".format(outbase))\n",
    "\n",
    "outburn = os.path.join(outfolder,\"plot-\" + workbook + \"-era-burn.png\")\n",
    "plotd.plot_lines(filtered_dates, mean_s1burn, std_s1burn, filtered_edata, mean_edata, std_edata, features, pnames, eraname, outburn, burn_start, graphics = False)\n",
    "print(\"Generated: {}\".format(outburn))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1175bb236b764cf1b276feaa950f32ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(Text(value='Non-burn Backscatter polygon'), Text(value='Plotting for rorbit  and…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Display plots side by side if have both non-burn and burn data extracted\n",
    "\n",
    "# read plot images\n",
    "if bdata:\n",
    "    img1 = open(outbase, 'rb').read()\n",
    "    img2 = open(outburn, 'rb').read()\n",
    "\n",
    "    # Text widgets\n",
    "    t1 = widgets.Text(value='Non-burn %s polygon'%(workbook))\n",
    "    t2 = widgets.Text(value='Plotting for rorbit {} and direction {}'.format(chosen_orbit, chosen_direction))\n",
    "    t3 = widgets.Text(value='Burn %s polygon'%(workbook))\n",
    "    ## Create image widgets.\n",
    "    i1 = widgets.Image(value=img1, width='50%')\n",
    "    i2 = widgets.Image(value=img2, width='50%')\n",
    "    ## Side by side thanks to HBox widget\n",
    "    titles = widgets.HBox([t1, t2, t3], width='100%')\n",
    "    plots = widgets.HBox([i1, i2], width='100%')\n",
    "\n",
    "else:\n",
    "    img1 = open(outburn, 'rb').read()\n",
    "    \n",
    "    # Text widgets\n",
    "    t1 = widgets.Text(value='Burn %s polygon'%(workbook), width='20%')\n",
    "    t2 = widgets.Text(value='Plotting for rorbit {} and direction {}'.format(chosen_orbit, chosen_direction), width='80%')\n",
    "    ## Create image widgets.\n",
    "    i1 = widgets.Image(value=img1, width='100%')\n",
    "    ## Side by side thanks to HBox widget\n",
    "    titles = widgets.HBox([t1, t2], width='100%')\n",
    "    plots = widgets.HBox([i1], width='100%')\n",
    "\n",
    "# User interface\n",
    "ui = widgets.VBox([titles,plots])\n",
    "## Finally, show.\n",
    "display(ui)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Python [conda env:uburn_env]",
   "language": "python",
   "name": "conda-env-uburn_env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}