Analyse-Sentinel1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# To be run once to setup modules and paths\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n",
      "Defining variables as not called externally\n",
      "Loaded configuration: {'sharedfolder': 'my_shared_data_folder'}.\n"
     ]
    }
   ],
   "source": [
    "#Reload modules without shutting notebook\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import shutil\n",
    "import glob\n",
    "import numpy as np\n",
    "import h5py\n",
    "import pickle\n",
    "from pandas import DataFrame\n",
    "import gdal\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry.polygon import Polygon\n",
    "from PIL import Image, ImageDraw\n",
    "\n",
    "# Do not show warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Upland burn modules\n",
    "import utils.extract_aoi as extract\n",
    "import utils.plot_data as plotd\n",
    "from utils.get_configuration import get_configuration\n",
    "\n",
    "# Define input variables if not calling externally\n",
    "try:\n",
    "    print(\"Input dataset from: {}\".format(basefolder))     # check if defined\n",
    "    polygon_file = os.path.join(datasets, polygon + \".geojson\")\n",
    "\n",
    "# if not defined raises an error and control shifts to except block. \n",
    "except:\n",
    "    print (\"Defining variables as not called externally\") \n",
    "    \n",
    "    # Get configuration information\n",
    "    cstudy = \"skye\" \n",
    "\n",
    "    # Get configuration information\n",
    "    basefolder, s1ardfolder, datasets, outfolder, tmpfolder, ofiles, hdfile, pfile, verbose, graphics = get_configuration(cstudy)\n",
    "    \n",
    "    # Burn Polygon\n",
    "    polygon = \"Skye_burn_extent_163\"\n",
    "    polygon_file = os.path.join(datasets, polygon + \".geojson\")\n",
    "    \n",
    "# Output array\n",
    "burn_subarray = None\n",
    "non_burn_subarray = None\n",
    "#print(\"Check: \",ofiles[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Display generated Skye subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load pickle file and extract statistics\n",
    "infile = open(pfile,'rb')\n",
    "data_dict = pickle.load(infile)\n",
    "infile.close()\n",
    "\n",
    "# Convert to dataframe\n",
    "df = DataFrame(data_dict, columns=['layer', 'sensor', 'date', 'stime', 'etime','polarisation', 'rorbit', 'direction'])\n",
    "layers = len(df)\n",
    "if verbose:\n",
    "    print(\"Loading {} layers\".format(layers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "# Read hdf file and display chosen band: VV or VH\n",
    "if verbose:\n",
    "    hf_object = h5py.File(hdfile, 'r')\n",
    "    d = hf_object['s1ard']\n",
    "    print(\"Displaying {} layers\".format(layers))\n",
    "    plotd.plot_images(d, df['date'].values, verb = verbose, hdf = 'VV')\n",
    "    hf_object.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create subset for a polygon of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import geojson as gj\n",
    "def create_subset(polygon_file, subarray, ofiles, hdfile):\n",
    "    # Get coords for bounding box\n",
    "    with open (polygon_file, 'r') as f:\n",
    "        loadg = gj.loads(f.read())\n",
    "    x, y = zip(*gj.utils.coords(loadg))\n",
    "    min_x, max_x, min_y, max_y = min(x), max(x), min(y), max(y)\n",
    "    if verbose:\n",
    "        print(\"Input for GeoJSON ULeft {}:{} LRight {}:{} \".format(min_x, max_y, max_x, min_y))\n",
    "\n",
    "    # Get coordinate information from GeoTIFF file\n",
    "    if verbose:\n",
    "        print(\"Loading {} into HDF file\".format(ofiles[0]))\n",
    "    ds = gdal.Open(ofiles[0], gdal.GA_ReadOnly)\n",
    "    bands = ds.RasterCount\n",
    "    image = ds.GetRasterBand(1).ReadAsArray()\n",
    "    xdim, ydim = image.shape\n",
    "    del image\n",
    "\n",
    "    global pixelWidth\n",
    "    global pixelHeight\n",
    "    \n",
    "    # Getting georeference info\n",
    "    transform = ds.GetGeoTransform()\n",
    "    projection = ds.GetProjection()\n",
    "    xOrigin = transform[0] # top left x \n",
    "    yOrigin = transform[3] # top left y \n",
    "    pixelWidth = transform[1] # w-e pixel resolution \n",
    "    pixelHeight = -transform[5] # n-s pixel resolution (negative value) \n",
    "\n",
    "    # Close HDF file\n",
    "    ds = None\n",
    "\n",
    "    # Computing [top left] [lower right]\n",
    "    i1 = int((min_x - xOrigin) / pixelWidth)\n",
    "    i2 = int((max_x - xOrigin) / pixelWidth) + 1 #we need to add 1 as int rounds down\n",
    "    j1 = int((yOrigin - max_y) / pixelHeight)\n",
    "    j2 = int((yOrigin - min_y) / pixelHeight) + 1 # we need to add 1 as int rounds down\n",
    "    if verbose:\n",
    "        print(\"Pixel coordinates for subset ULeft {}:{} LRight {}:{}\".format(i1, j1, i2, j2))\n",
    "\n",
    "    hf_object = h5py.File(hdfile, 'r')\n",
    "\n",
    "    vect = []\n",
    "    # first lets get the coordinates of the shape outline\n",
    "    for x1, y1 in zip(x, y):\n",
    "        xpix = ((x1 - xOrigin) / pixelWidth) - i1\n",
    "        ypix = ((yOrigin - y1) / pixelHeight) - j1\n",
    "        #vect.append([int(ypix), int(xpix)])\n",
    "        vect.append(int(xpix))\n",
    "        vect.append(int(ypix))\n",
    "    img = Image.new('L',  (i2-i1,j2-j1), 0)                                                                                                                                                                                                                  \n",
    "    draw = ImageDraw.Draw(img)                                                                                                                                                                                                                                        \n",
    "    draw.polygon(vect, fill=1)\n",
    "    mask = np.array(img)\n",
    "\n",
    "    if verbose:\n",
    "        print(\"Loading layers...\")\n",
    "    if subarray is None:\n",
    "        # Using 4GB of chunk_cache_mem here (\"rdcc_nbytes\")\n",
    "        with h5py.File(hdfile, 'r', rdcc_nbytes =1024**2*4000, rdcc_nslots=1e7) as h5f: \n",
    "            for layer in range(layers):\n",
    "                print(\"Loading layer {} of {}\".format(layer+1,layers), end='\\r')\n",
    "                if subarray is None:\n",
    "                    subarray = h5f['s1ard'][layer, :, j1:j2, i1:i2] # only extract the layer and subset\n",
    "                    #if np.mean(subarray) == 0:\n",
    "                    #    subarray = None\n",
    "                    #    df.drop(layer, inplace=True)\n",
    "                    #    continue\n",
    "                    subarray = np.expand_dims(subarray, axis=0)\n",
    "                    continue\n",
    "                else:\n",
    "                    ind_arr = h5f['s1ard'][layer, :, j1:j2, i1:i2]\n",
    "                    #if np.mean(ind_arr) == 0:\n",
    "                    #    df.drop(layer, inplace=True)\n",
    "                    #    continue\n",
    "                    ind_arr = np.expand_dims(ind_arr, axis=0)\n",
    "                subarray = np.vstack([subarray, ind_arr]) # stack them \n",
    "    else:\n",
    "        print(\"Layers already loaded: \",subarray.shape)\n",
    "\n",
    "    print(\"\\n Extracting only polygon area...\")\n",
    "    mask =  np.broadcast_to(mask, subarray.shape) # make the mask the same shape\n",
    "    #subarray[mask==0] = np.nan # set anything outside the polygon to np.nan\n",
    "\n",
    "    # New upper-left X,Y values\n",
    "    start_x = xOrigin + (i1 * pixelWidth)\n",
    "    start_y = yOrigin - (j1 * pixelHeight)\n",
    "    end_x = xOrigin + (i2 * pixelWidth)\n",
    "    end_y = yOrigin - (j2 * pixelHeight)\n",
    "    new_transform = (start_x, transform[1], transform[2], start_y, transform[4], transform[5])\n",
    "    if verbose:\n",
    "        print(\"Origin: Old-UL {:.3f} {:.3f} New-UL {:.3f} {:.3f} New-LR {:.3f} {:.3f}\".\n",
    "              format(xOrigin, yOrigin, start_x, start_y, end_x, end_y))\n",
    "    return subarray, new_transform, mask[0,0,:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6b3f2097cb3744d383fbd82223b5df5e",
       "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": "efb5d25835444a64805da60c9e4fbc1c",
       "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"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "polar_dict = {'VV' : 0, 'VH' : 1}\n",
    "\n",
    "list_of_dates = df['date'].values\n",
    "infoW = widgets.Label(\"\")\n",
    "\n",
    "#subarray = None\n",
    "burn_subarray, new_transform, mask = create_subset(polygon_file, burn_subarray, ofiles, hdfile)\n",
    "\n",
    "if bdata:\n",
    "    polygon = non_burn\n",
    "    polygon_non_burn_file = os.path.join(datasets, polygon + \".geojson\")\n",
    "    non_burn_subarray, non_burn_transform, non_burn_mask = create_subset(polygon_non_burn_file, non_burn_subarray, ofiles, hdfile)\n",
    "    non_burn_df = df\n",
    "    \n",
    "def save_image():\n",
    "    dst_filename = os.path.join(outfolder, \"%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], date, pol_text, \"burn-backscatter\"))\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(new_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",
    "    global date\n",
    "    global pol_text\n",
    "    global curr_arr\n",
    "    pol_text = polarisation\n",
    "    pol = polar_dict[polarisation]\n",
    "    date = list_of_dates[img_num-1]\n",
    "    curr_arr = np.copy(burn_subarray[(img_num-1), pol, :, :])\n",
    "    curr_arr[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)\n",
    "    plt.title(\"Image Date: %s\"%(date))\n",
    "    \n",
    "    curr_arr2= np.copy(burn_subarray[(img_num-1), pol, :, ])\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    ax2.imshow(curr_arr2)\n",
    "    ax2.contour(mask > 0, levels=[2], colors=['r'])\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "savebutton = interactive(save_image, {'manual' : True, 'manual_name' : 'Save'})\n",
    "\n",
    "display(savebutton)\n",
    "\n",
    "def save_image2():\n",
    "    dst_filename = os.path.join(outfolder, \"%s_%s_%s_%s_%s.tif\"%(polygon.split(\"_\")[0], polygon.split(\"_\")[3], date2, pol_text, \"non-burn-backscatter\"))\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",
    "    global date2\n",
    "    global pol_text\n",
    "    global curr_arr3\n",
    "    pol_text = polarisation_non_burn\n",
    "    pol = polar_dict[polarisation_non_burn]\n",
    "    date2 = list_of_dates[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)\n",
    "    plt.title(\"Image Date: %s\"%(date))\n",
    "    \n",
    "    curr_arr4= np.copy(non_burn_subarray[(img_num_non_burn-1), pol, :, ])\n",
    "    ax4 = fig2.add_subplot(122)\n",
    "    ax4.imshow(curr_arr4)\n",
    "    ax4.contour(non_burn_mask > 0, levels=[2], colors=['r'])\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "savebutton2 = interactive(save_image2, {'manual' : True, 'manual_name' : 'Save'})\n",
    "display(savebutton2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Python [conda env:rsgislib_dev]",
   "language": "python",
   "name": "conda-env-rsgislib_dev-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
}