<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%env GDAL_DATA = /usr/share/gdal/2.2/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Water Across Synthetic Aperature Radar (WASARD)\n",
    "WASARD is a general purpose transfer model between optical and SAR imagery for water classification.   \n",
    "A trained WASARD model can be executed on SAR imagery to create water classification maps over a region.  \n",
    "\n",
    "This notebook is inspired by an IGARSS publication titled [Water Across Synthetic Aperture Radar Data (WASARD): SAR Water Body Classification for the Open Data Cube](https://www.igarss2018.org/Papers/viewpapers.asp?papernum=3380) authored by <NAME>, <NAME>, <NAME>.  \n",
    "\n",
    "\n",
    "WASARD is trained using water classifications on optical imagery as a point of reference.  A machine learning model is used to approximate a transfer function between SAR data and Optical Classifications.  \n",
    "\n",
    "\n",
    ">### Details regarding WASARD in this notebook  \n",
    "> `Transfer Model:` Linear SVM  \n",
    "> `Optical water classifier:` WOFS  \n",
    "> `Optical Source:` Landsat 8  \n",
    "> `SAR Target`: Sentinel 1a  \n",
    "\n",
    "\n",
    "# Training process  \n",
    "![](./diagrams/WASARD/flow_chart.png)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# In this notebook\n",
    "* Existing coordinates are loaded for `Sentinel1` and `Landsat8` Imagery\n",
    "* Upon comparison,coordinates are used to inform how data should be loaded.  \n",
    "* `Sentinel1` imagery is loaded.\n",
    "* `Landsat8`  imagery is loaded, reprojected, and upsampled to match the resolution of `Sentine1` \n",
    "* an existing water classifier is run on optical data.  \n",
    "* optical water classifications are averaged across time and visualized\n",
    "* a wasard model is created and trained  \n",
    "* a wasard model is run on sentinel data\n",
    "* the output of wasard is averaged across time and visualized \n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Retreiving/Cleaning Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connect to the data cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datacube\n",
    "dc = datacube.Datacube(app = \"[notebook][wasard][samoa]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Extents of Analysis  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sar_product_name = \"s1g_gamma0_samoa\"\n",
    "optical_product_name = \"ls8_lasrc_samoa\"\n",
    "\n",
    "# #Sa'anapu\n",
    "# longitude_extent = (-171.904492, -171.790327)\n",
    "# latitude_extent = (-14.0  , -13.962341)\n",
    "\n",
    "#Apia\n",
    "latitude_extent = -13.853425,-13.815715\n",
    "longitude_extent =-171.787842, -171.681356\n",
    "\n",
    "date_range = ('2016-8-1','2018-3-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_display_map import display_map\n",
    "display_map(latitude = latitude_extent, longitude = longitude_extent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine Coordinate Extents "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_coordinates =  dc.load(product = sar_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                measurements = [])\n",
    "landsat_coordinates  =  dc.load(product = optical_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                measurements = [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils.data_cube_utilities.xarray_bokeh_plotting as xr_bokeh\n",
    "xr_bokeh.init_notebook() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "xr_bokeh.dim_alignement( sentinel_coordinates.isel(latitude = slice(0,100), longitude = slice(0,100)), \" Sentinel 1\",\n",
    "                          landsat_coordinates.isel(latitude = slice(0,100), longitude = slice(0,100)),  \"Landsat 8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and reproject  \n",
    "\n",
    "Add bit about reprojecting landsat the way I do with sentinel.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Sentinel CRS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sentinel_details = dc.list_products()[dc.list_products()[\"name\"].str.contains(sar_product_name)]\n",
    "sentinel_details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_crs        = str(sentinel_details['crs'].values[0])\n",
    "sentinel_resolution = tuple(sentinel_details['resolution'].values[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (sentinel_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load both SAR and Optical imagery "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_dataset     =  dc.load(product = sar_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "landsat_dataset     =    dc.load(product = optical_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                output_crs = sentinel_crs,\n",
    "                                resolution = sentinel_resolution\n",
    "                                )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### compare reprojected coordinates for alignment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_landsat_coords = landsat_dataset.isel(latitude = slice(0,100),\n",
    "                                        longitude = slice(0,100),\n",
    "                                        time = 0).coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_sentinel_coords = sentinel_dataset.isel(latitude = slice(0,100),\n",
    "                                        longitude = slice(0,100),\n",
    "                                        time = 0).coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "xr_bokeh.dim_alignement(subset_of_sentinel_coords, \"S1\",\n",
    "                        subset_of_landsat_coords, \"LS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean SAR Data  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr  \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def remove_all_zero(dataset):\n",
    "    return dataset.drop([c[0].values \n",
    "        for c in [(t,np.count_nonzero(dataset.sel(time=t).vv)) \n",
    "                  for t in dataset.time] if c[1]==0],dim='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List  \n",
    "import itertools\n",
    "\n",
    "has_time_dimension = lambda x: \"time\" in dict(x.dims).keys()\n",
    "get_first = lambda x: x[0]\n",
    "get_last =  lambda x: x[-1]\n",
    "\n",
    "\n",
    "def group_dates_by_day( dates: List[np.datetime64]) -> List[List[np.datetime64]]:\n",
    "    generate_key = lambda b: ((b - np.datetime64('1970-01-01T00:00:00Z')) / (np.timedelta64(1, 'h')*24)).astype(int)\n",
    "    return [list(group) for key, group in itertools.groupby(dates, key=generate_key)]\n",
    "\n",
    "def reduce_on_day(ds: xr.Dataset,\n",
    "                  reduction_func: np.ufunc = np.nanmean) -> xr.Dataset:\n",
    "    #Group dates by day into date_groups\n",
    "    day_groups = group_dates_by_day(ds.time.values)\n",
    "    \n",
    "    #slice large dataset into many smaller datasets by date_group\n",
    "    group_chunks = (ds.sel(time = t) for t in day_groups)\n",
    "    \n",
    "    #reduce each dataset using something like \"average\" or \"median\" such that many values for a day become one value   \n",
    "    group_slices = (_ds.reduce(reduction_func, dim = \"time\") for _ds in group_chunks if has_time_dimension(_ds))\n",
    "\n",
    "    #recombine slices into larger dataset\n",
    "    new_dataset  = xr.concat(group_slices, dim = \"time\") \n",
    "    \n",
    "    #rename times values using the first time in each date_group  \n",
    "    new_times = list(map(get_first, day_groups))    \n",
    "    new_dataset = new_dataset.reindex(dict(time = np.array(new_times)))\n",
    "    \n",
    "    return new_dataset\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_dataset = remove_all_zero(sentinel_dataset)\n",
    "sentinel_dataset = reduce_on_day(sentinel_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_landsat_coords = landsat_dataset.isel(latitude = 0, longitude = 0).coords\n",
    "subset_of_sentinel_coords = sentinel_dataset.isel(latitude = 0, longitude = 0).coords\n",
    "xr_bokeh.dim_alignement(subset_of_sentinel_coords, \"S1\",\n",
    "                        subset_of_landsat_coords, \"LS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examine Optical Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_mosaic import ls8_unpack_qa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_mosaic import ls8_unpack_qa\n",
    "from utils.data_cube_utilities.dc_mosaic import create_median_mosaic\n",
    "import xarray as xr  \n",
    "import numpy as np\n",
    "\n",
    "def clean_mask_ls8(ds:xr.Dataset) -> np.array:\n",
    "    clear_xarray  = ls8_unpack_qa(ds.pixel_qa, \"clear\")  \n",
    "    water_xarray  = ls8_unpack_qa(ds.pixel_qa, \"water\")\n",
    "\n",
    "    cloud_free_boolean_mask = np.logical_or(clear_xarray, water_xarray)\n",
    "    return cloud_free_boolean_mask\n",
    "\n",
    "def median_mosaic_ls8(dataset):\n",
    "    # The mask here is based on pixel_qa products. It comes bundled in with most Landsat Products.\n",
    "        return create_median_mosaic(dataset, clean_mask = clean_mask_ls8(dataset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mosaic = median_mosaic_ls8(landsat_dataset.isel(time = slice(0,20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_rgb import rgb\n",
    "rgb(mosaic, bands=['red', 'green', 'blue'], width= 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Classification: Optical Imagery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Import water classifier\n",
    "> This notebooks uses WOFS. A water classification algorithm for landsat imagery, developed by Geoscience Australia. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_water_classifier import wofs_classify"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Classify\n",
    "> `wofs_classify` is used. An additional mask is passed in that indicates the presence of clouds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "water_classifications = wofs_classify(landsat_dataset, clean_mask= clean_mask_ls8(landsat_dataset), no_data= np.nan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Visualize optical water classifciaton\n",
    ">A percentage (between the range of [0-1]) denotes the frequency of water in any given area, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "def aspect_ratio_helper(x,y, fixed_width = 20):\n",
    "    width = fixed_width\n",
    "    height = y * (fixed_width / x)\n",
    "    return (width, height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  \n",
    "\n",
    "plt.figure(figsize = aspect_ratio_helper(*reversed(list(water_classifications.wofs.shape)[1:])))\n",
    "water_classifications.mean(dim = 'time').wofs.plot(cmap = \"jet_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train WASARD water classifier using optical data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Import wasard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities import wasard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Train WASARD classifier\n",
    "> In this case only data from samoa will be used. In the general case, It is suggested that a geographicall varied dataset is used to train the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samoa_classifier = wasard.wasard_classifier(sar_dataset=sentinel_dataset,\n",
    "                                              landsat_dataset=landsat_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Classifcation: WASARD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Classify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samoa_classified = samoa_classifier.wasard_classify(sentinel_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Visualize WASARD SAR water classifciaton\n",
    ">A percentage (between the range of [0-1]) denotes the frequency of water in any given area, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = aspect_ratio_helper(*reversed(list(water_classifications.wofs.shape)[1:])))\n",
    "samoa_classified.wasard.mean(dim = \"time\").plot(cmap = \"jet_r\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}<commit_msg>Changed S1_WASARD.ipynb to use data available on the sandbox.<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Across Synthetic Aperature Radar (WASARD)\n",
    "WASARD is a general purpose transfer model between optical and SAR imagery for water classification.   \n",
    "A trained WASARD model can be executed on SAR imagery to create water classification maps over a region.  \n",
    "\n",
    "This notebook is inspired by an IGARSS publication titled [Water Across Synthetic Aperture Radar Data (WASARD): SAR Water Body Classification for the Open Data Cube](https://www.igarss2018.org/Papers/viewpapers.asp?papernum=3380) authored by <NAME>, <NAME>, <NAME>.  \n",
    "\n",
    "\n",
    "WASARD is trained using water classifications on optical imagery as a point of reference.  A machine learning model is used to approximate a transfer function between SAR data and Optical Classifications.  \n",
    "\n",
    "\n",
    ">### Details regarding WASARD in this notebook  \n",
    "> `Transfer Model:` Linear SVM  \n",
    "> `Optical water classifier:` WOFS  \n",
    "> `Optical Source:` Landsat 8  \n",
    "> `SAR Target`: Sentinel 1a  \n",
    "\n",
    "# In this notebook\n",
    "* Existing coordinates are loaded for `Sentinel1` and `Landsat8` Imagery\n",
    "* Upon comparison,coordinates are used to inform how data should be loaded.  \n",
    "* `Sentinel1` imagery is loaded.\n",
    "* `Landsat8`  imagery is loaded, reprojected, and upsampled to match the resolution of `Sentine1` \n",
    "* an existing water classifier is run on optical data.  \n",
    "* optical water classifications are averaged across time and visualized\n",
    "* a wasard model is created and trained  \n",
    "* a wasard model is run on sentinel data\n",
    "* the output of wasard is averaged across time and visualized \n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Retreiving/Cleaning Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connect to the data cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datacube\n",
    "dc = datacube.Datacube(app = \"[notebook][wasard][samoa]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Extents of Analysis  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sar_product_name = \"s1g_gamma0\"\n",
    "optical_product_name = \"ls8_usgs_sr_scene\"\n",
    "\n",
    "# #Sa'anapu\n",
    "# longitude_extent = (-171.904492, -171.790327)\n",
    "# latitude_extent = (-14.0  , -13.962341)\n",
    "\n",
    "#Apia\n",
    "latitude_extent = -13.853425,-13.815715\n",
    "longitude_extent =-171.787842, -171.681356\n",
    "\n",
    "date_range = ('2016-8-1','2018-3-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_display_map import display_map\n",
    "display_map(latitude = latitude_extent, longitude = longitude_extent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine Coordinate Extents "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_coordinates =  dc.load(product = sar_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                measurements = [])\n",
    "landsat_coordinates  =  dc.load(product = optical_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                measurements = [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import utils.data_cube_utilities.xarray_bokeh_plotting as xr_bokeh\n",
    "xr_bokeh.init_notebook() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr_bokeh.dim_alignement( sentinel_coordinates.isel(latitude = slice(0,100), longitude = slice(0,100)), \" Sentinel 1\",\n",
    "                          landsat_coordinates.isel(latitude = slice(0,100), longitude = slice(0,100)),  \"Landsat 8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and reproject  \n",
    "\n",
    "Add bit about reprojecting landsat the way I do with sentinel.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Sentinel CRS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sentinel_details = dc.list_products()[dc.list_products()[\"name\"].str.contains(sar_product_name)]\n",
    "sentinel_details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_crs        = str(sentinel_details['crs'].values[0])\n",
    "sentinel_resolution = tuple(sentinel_details['resolution'].values[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (sentinel_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load both SAR and Optical imagery "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_dataset     =  dc.load(product = sar_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "landsat_dataset     =    dc.load(product = optical_product_name,\n",
    "                                latitude = latitude_extent,\n",
    "                                longitude = longitude_extent,\n",
    "                                time = date_range,\n",
    "                                output_crs = sentinel_crs,\n",
    "                                resolution = sentinel_resolution\n",
    "                                )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### compare reprojected coordinates for alignment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_landsat_coords = landsat_dataset.isel(latitude = slice(0,100),\n",
    "                                        longitude = slice(0,100),\n",
    "                                        time = 0).coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_sentinel_coords = sentinel_dataset.isel(latitude = slice(0,100),\n",
    "                                        longitude = slice(0,100),\n",
    "                                        time = 0).coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr_bokeh.dim_alignement(subset_of_sentinel_coords, \"S1\",\n",
    "                        subset_of_landsat_coords, \"LS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean SAR Data  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr  \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def remove_all_zero(dataset):\n",
    "    return dataset.drop([c[0].values \n",
    "        for c in [(t,np.count_nonzero(dataset.sel(time=t).vv)) \n",
    "                  for t in dataset.time] if c[1]==0],dim='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List  \n",
    "import itertools\n",
    "\n",
    "has_time_dimension = lambda x: \"time\" in dict(x.dims).keys()\n",
    "get_first = lambda x: x[0]\n",
    "get_last =  lambda x: x[-1]\n",
    "\n",
    "\n",
    "def group_dates_by_day( dates: List[np.datetime64]) -> List[List[np.datetime64]]:\n",
    "    generate_key = lambda b: ((b - np.datetime64('1970-01-01T00:00:00Z')) / (np.timedelta64(1, 'h')*24)).astype(int)\n",
    "    return [list(group) for key, group in itertools.groupby(dates, key=generate_key)]\n",
    "\n",
    "def reduce_on_day(ds: xr.Dataset,\n",
    "                  reduction_func: np.ufunc = np.nanmean) -> xr.Dataset:\n",
    "    #Group dates by day into date_groups\n",
    "    day_groups = group_dates_by_day(ds.time.values)\n",
    "    \n",
    "    #slice large dataset into many smaller datasets by date_group\n",
    "    group_chunks = (ds.sel(time = t) for t in day_groups)\n",
    "    \n",
    "    #reduce each dataset using something like \"average\" or \"median\" such that many values for a day become one value   \n",
    "    group_slices = (_ds.reduce(reduction_func, dim = \"time\") for _ds in group_chunks if has_time_dimension(_ds))\n",
    "\n",
    "    #recombine slices into larger dataset\n",
    "    new_dataset  = xr.concat(group_slices, dim = \"time\") \n",
    "    \n",
    "    #rename times values using the first time in each date_group  \n",
    "    new_times = list(map(get_first, day_groups))    \n",
    "    new_dataset = new_dataset.reindex(dict(time = np.array(new_times)))\n",
    "    \n",
    "    return new_dataset\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sentinel_dataset = remove_all_zero(sentinel_dataset)\n",
    "sentinel_dataset = reduce_on_day(sentinel_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_of_landsat_coords = landsat_dataset.isel(latitude = 0, longitude = 0).coords\n",
    "subset_of_sentinel_coords = sentinel_dataset.isel(latitude = 0, longitude = 0).coords\n",
    "xr_bokeh.dim_alignement(subset_of_sentinel_coords, \"S1\",\n",
    "                        subset_of_landsat_coords, \"LS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examine Optical Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_mosaic import ls8_unpack_qa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_mosaic import ls8_unpack_qa\n",
    "from utils.data_cube_utilities.dc_mosaic import create_median_mosaic\n",
    "import xarray as xr  \n",
    "import numpy as np\n",
    "\n",
    "def clean_mask_ls8(ds:xr.Dataset) -> np.array:\n",
    "    clear_xarray  = ls8_unpack_qa(ds.pixel_qa, \"clear\")  \n",
    "    water_xarray  = ls8_unpack_qa(ds.pixel_qa, \"water\")\n",
    "\n",
    "    cloud_free_boolean_mask = np.logical_or(clear_xarray, water_xarray)\n",
    "    return cloud_free_boolean_mask\n",
    "\n",
    "def median_mosaic_ls8(dataset):\n",
    "    # The mask here is based on pixel_qa products. It comes bundled in with most Landsat Products.\n",
    "        return create_median_mosaic(dataset, clean_mask = clean_mask_ls8(dataset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mosaic = median_mosaic_ls8(landsat_dataset.isel(time = slice(0,20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_rgb import rgb\n",
    "rgb(mosaic, bands=['red', 'green', 'blue'], width= 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Classification: Optical Imagery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Import water classifier\n",
    "> This notebooks uses WOFS. A water classification algorithm for landsat imagery, developed by Geoscience Australia. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities.dc_water_classifier import wofs_classify"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Classify\n",
    "> `wofs_classify` is used. An additional mask is passed in that indicates the presence of clouds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "water_classifications = wofs_classify(landsat_dataset, clean_mask= clean_mask_ls8(landsat_dataset), no_data= np.nan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Visualize optical water classifciaton\n",
    ">A percentage (between the range of [0-1]) denotes the frequency of water in any given area, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "def aspect_ratio_helper(x,y, fixed_width = 20):\n",
    "    width = fixed_width\n",
    "    height = y * (fixed_width / x)\n",
    "    return (width, height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  \n",
    "\n",
    "plt.figure(figsize = aspect_ratio_helper(*reversed(list(water_classifications.wofs.shape)[1:])))\n",
    "water_classifications.mean(dim = 'time').wofs.plot(cmap = \"jet_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train WASARD water classifier using optical data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Import wasard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.data_cube_utilities import wasard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Train WASARD classifier\n",
    "> In this case only data from samoa will be used. In the general case, It is suggested that a geographicall varied dataset is used to train the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samoa_classifier = wasard.wasard_classifier(sar_dataset=sentinel_dataset,\n",
    "                                              landsat_dataset=landsat_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Classifcation: WASARD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Classify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samoa_classified = samoa_classifier.wasard_classify(sentinel_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Visualize WASARD SAR water classifciaton\n",
    ">A percentage (between the range of [0-1]) denotes the frequency of water in any given area, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = aspect_ratio_helper(*reversed(list(water_classifications.wofs.shape)[1:])))\n",
    "samoa_classified.wasard.mean(dim = \"time\").plot(cmap = \"jet_r\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
{'fraction_non_alphanumeric': 0.17594038115997526, 'fraction_numerical': 0.011752923502901412, 'mean_word_length': 2.4310257705912077, 'pattern_counts': {'":': 477, '<': 11, '<?xml version=': 0, '>': 51, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 64, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2172113', 'n_tokens_mistral': 12512, 'n_tokens_neox': 11898, 'n_words': 2371}