OPERA-RTC-S1-time_series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false
   },
   "source": [
    "# Time Series Analysis with OPERA RTC-S1\n",
    "\n",
    "### Franz J Meyer; University of Alaska Fairbanks & Josef Kellndorfer, [Earth Big Data, LLC](http://earthbigdata.com/)\n",
    "\n",
    "This notebook will introduce you to the analysis of deep multi-temporal SAR image data stacks in the framework of *Jupyter Notebooks*. The Jupyter Notebook environment is easy to launch in any web browser for interactive data exploration with provided or new training data. Notebooks are comprised of text written in a combination of executable python code and markdown formatting including latex style mathematical equations. Another advantage of Jupyter Notebooks is that they can easily be expanded, changed, and shared with new data sets or newly available time series steps. Therefore, they provide an excellent basis for collaborative and repeatable data analysis.\n",
    "\n",
    "**We introduce the following data analysis concepts:**\n",
    "\n",
    "- How to load your own SAR data into Jupyter Notebooks and create a time series stack \n",
    "- How to apply calibration constants to covert initial digital number (DN) data into calibrated radar cross section information.\n",
    "- How to subset images and create a time series of your subset data.\n",
    "- How to explore the time-series information in SAR data stacks for environmental analysis.\n",
    "\n",
    "<div class=\"alert alert-info\" style=\"display: flex; align-items: center; font-family: 'Times New Roman', Times, serif; background-color: #d1ecf1;\">\n",
    "  <div style=\"display: flex; align-items: center; width: 10%;\">\n",
    "    <a href=\"https://github.com/ASFOpenSARlab/opensarlab_OPERA-CSLC_Recipe_Book/issues\">\n",
    "      <img src=\"https://opensarlab-docs.asf.alaska.edu/opensarlab-notebook-assets/logos/github_issues.png\" alt=\"GitHub logo over the word Issues\" style=\"width: 100%;\">\n",
    "    </a>\n",
    "  </div>\n",
    "  <div style=\"width: 95%;\">\n",
    "    <b>Did you find a bug? Do you have a feature request?</b>\n",
    "    <br/>\n",
    "    Explore GitHub Issues on this Jupyter Book's GitHub repository. Find solutions, add to the discussion, or start a new bug report or feature request: <a href=\"https://github.com/ASFOpenSARlab/opensarlab_OPERA-RTC-S1_Recipe_Book/issues\">opensarlab_OPERA-RTC-S1_Recipe_Book Issues</a>\n",
    "  </div>\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\" style=\"display: flex; align-items: center; justify-content: space-between; font-family: 'Times New Roman', Times, serif; background-color: #d1ecf1;\">\n",
    "  <div style=\"display: flex; align-items: center; width: 10%; margin-right: 10px;\">\n",
    "    <a href=\"mailto:uso@asf.alaska.edu\">\n",
    "      <img src=\"https://opensarlab-docs.asf.alaska.edu/opensarlab-notebook-assets/logos/ASF_support_logo.png\" alt=\"ASF logo\" style=\"width: 100%\">\n",
    "    </a>\n",
    "  </div>\n",
    "  <div style=\"width: 95%;\">\n",
    "    <b>Have a question related to SAR, OPERA RTCs, or ASF data access?</b>\n",
    "    <br/>\n",
    "    Contact ASF User Support: <a href=\"mailto:uso@asf.alaska.edu\">uso@asf.alaska.edu</a>\n",
    "  </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 0. Importing Relevant Python Packages\n",
    "\n",
    "In this notebook we will use the following scientific libraries:\n",
    "\n",
    "- [Pandas](https://pandas.pydata.org/) is a Python library that provides high-level data structures and a vast variety of tools for analysis. The great feature of this package is the ability to translate rather complex operations with data into one or two commands. Pandas contains many built-in methods for filtering and combining data, as well as the time-series functionality.\n",
    "- [GDAL](https://www.gdal.org/) is a software library for reading and writing raster and vector geospatial data formats. It includes a collection of programs tailored for geospatial data processing. Most modern GIS systems (such as ArcGIS or QGIS) use GDAL in the background.\n",
    "- [NumPy](http://www.numpy.org/) is one of the principal packages for scientific applications of Python. It is intended for processing large multidimensional arrays and matrices, and an extensive collection of high-level mathematical functions and implemented methods makes it possible to perform various operations with these objects.\n",
    "- [XArray](https://docs.xarray.dev/en/stable/index.html) provides a more intuitive way to load and view multi-dimensional geospatial arrays. It builds on Numpy/GDAL and we will use it to load in our geospatial data.\n",
    "- [Matplotlib](https://matplotlib.org/index.html) is a low-level library for creating two-dimensional diagrams and graphs. With its help, you can build diverse charts, from histograms and scatterplots to non-Cartesian coordinates graphs. Moreover, many popular plotting libraries are designed to work in conjunction with matplotlib.\n",
    "- [SciPy](https://www.scipy.org/about.html) is a library that provides functions for numerical integration, interpolation, optimization, linear algebra and statistics.\n",
    "\n",
    "**Our first step is to import them:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "import re\n",
    "from pathlib import Path\n",
    "import json # for loads\n",
    "import math # for ceil\n",
    "from datetime import datetime\n",
    "\n",
    "import pandas as pd # for DatetimeIndex\n",
    "from osgeo import gdal # for Info\n",
    "gdal.UseExceptions()\n",
    "import numpy as np # for copy, isnan, log10, ma.masked_where, max, mean, min, percentile, power, unique, var, where \n",
    "import scipy.signal\n",
    "import xarray as xr\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pylab as plb # for figure, grid, rcParams, savefig\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "from matplotlib import rc\n",
    "\n",
    "from ipyfilechooser import FileChooser\n",
    "\n",
    "from IPython.display import HTML\n",
    "\n",
    "import opensarlab_lib as asfn\n",
    "asfn.jupytertheme_matplotlib_format()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 1. Load Your Prepared Data Stack Into the Notebook\n",
    "\n",
    "**To download a set of data you can use to run this notebook, we recommend working through the [`search_download_OPERA-RTC-S1.ipynb`](../search_download_OPERA-RTC-S1.ipynb) notebook found in `data_access` folder in this repository!**\n",
    "\n",
    "You can however use any set of OPERA RTCs as long as they have the same footprint (e.g., the same burst ID)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Select the directory holding your tiffs**\n",
    "- Click the `Select` button\n",
    "- Navigate to your data directory\n",
    "- Click the `Select` button\n",
    "- Confirm that the desired path appears in green text\n",
    "- Click the `Change` button to alter your selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fc = FileChooser(Path.home())\n",
    "display(fc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Determine the path to the analysis directory containing the tiff directory:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_dir = Path(fc.selected_path)\n",
    "paths = list(analysis_dir.glob('*/*VV.tif')) + list(analysis_dir.glob('*/*VH.tif'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Load the data using Xarray**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function is used to correctly set the date and band properties of the data when using xr.open_mfdataset\n",
    "def add_time_and_pol(ds):\n",
    "    name = Path(ds.encoding['source']).with_suffix('').name\n",
    "    time = pd.to_datetime(name.split('_')[4])\n",
    "    polarization = name.split('_')[-1].upper()\n",
    "    ds = ds.assign_coords(time=time)\n",
    "    ds = ds.expand_dims(dim={'time':1})\n",
    "    ds.coords['band'] = [polarization]\n",
    "    return ds\n",
    "\n",
    "rtc = xr.open_mfdataset(paths, engine='rasterio', preprocess=add_time_and_pol).to_dataarray().squeeze()\n",
    "rtc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 3. Now You Can Work With Your Data\n",
    "\n",
    "Now you are ready to perform time series analysis on your data stack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 3.1 Open Your Data Stack and Visualize Some Layers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot images and histograms for VV and VH polarizations:**\n",
    "\n",
    "Note: Depending the histograms plotted by this cell, you may wish to adjust vmax when calling imshow() on ax1 and ax3. Increase the vmax value if the histogram cuts off much of the end of the peak, making your image too bright to see features well. Decrease vmax if the histogram extends much beyond the end of the peak, which will make your image appear dark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the pyplot plots\n",
    "fig = plb.figure(figsize=(18,10)) # Initialize figure with a size\n",
    "ax1 = fig.add_subplot(221)  # 221 determines: 2 rows, 2 plots, first plot\n",
    "ax2 = fig.add_subplot(222)  # 222 determines: 2 rows, 2 plots, second plot\n",
    "ax3 = fig.add_subplot(223)  # 223 determines: 2 rows, 2 plots, third plot\n",
    "ax4 = fig.add_subplot(224)  # 224 determines: 2 rows, 2 plots, fourth plot\n",
    "\n",
    "# Set vmin and vmax\n",
    "vmin=0\n",
    "vmax=0.2\n",
    "\n",
    "# Access the data\n",
    "time = rtc.time[2].values\n",
    "time_name = pd.to_datetime(str(time)).strftime('%Y-%m-%d')\n",
    "\n",
    "vv_array = rtc.sel(time=time, band='VV').squeeze().to_numpy()\n",
    "vh_array = rtc.sel(time=time, band='VH').squeeze().to_numpy()\n",
    "\n",
    "# Plot the VV image\n",
    "ax1.imshow(vv_array, cmap='gray', vmin=vmin, vmax=vmax) # see note above regarding vmax adjustments\n",
    "ax1.set_title(f'Image VV {time_name}')\n",
    "\n",
    "# Flatten the VV image into a 1 dimensional vector and plot the histogram:\n",
    "h = ax2.hist(vv_array.flatten(), bins=200, range=(vmin, vmax))\n",
    "ax2.xaxis.set_label_text('Amplitude? (Uncalibrated DN Values)')\n",
    "ax2.set_title(f'Histogram VV {time_name}')\n",
    "\n",
    "# Plot the VH image\n",
    "ax3.imshow(vh_array,cmap='gray', vmin=vmin, vmax=vmax) # see note above regarding vmax adjustments\n",
    "ax3.set_title(f'Image VH {time_name}')\n",
    "\n",
    "# Flatten the VH image into a 1 dimensional vector and plot the histogram:\n",
    "h = ax4.hist(vh_array.flatten(), bins=200, range=(vmin,vmax))\n",
    "ax4.xaxis.set_label_text('Amplitude? (Uncalibrated DN Values)')\n",
    "ax4.set_title(f'Histogram VH {time_name}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 3.2 Calibration and Data Conversion between dB and Power Scales\n",
    "\n",
    "**Note, that if your data were generated by HyP3, this step is not necessary!** HyP3 performs the full data calibration and provides you with calibrated data in power scale.\n",
    "    \n",
    "If, your data is from a different source, however, calibration may be necessary to ensure that image gray values correspond to proper radar cross section information. \n",
    "\n",
    "Calibration coefficients for SAR data are often defined in the decibel (dB) scale due to the high dynamic range of the imaging system. For the L-band ALOS PALSAR data at hand, the conversion from uncalibrated DN values to calibrated radar cross section values in dB scale is performed by applying a standard **calibration factor of -83 dB**. \n",
    "\n",
    "$\\gamma^0_{dB} = 20 \\cdot log10(DN) -83$\n",
    "\n",
    "The data at hand are radiometrically terrain corrected images, which are often expressed as terrain flattened $\\gamma^0$ backscattering coefficients. For forest and land cover monitoring applications $\\gamma^o$ is the preferred metric.\n",
    "\n",
    "**To apply the calibration constant for your data and export in *dB* scale, uncomment the following code cell:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# caldB = 20*np.log10(rtc)-83"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While **dB**-scaled images are often \"visually pleasing\", they are often not a good basis for mathematical operations on data. For instance, when we compute the mean of observations, it makes a difference whether we do that in power or dB scale. Since dB scale is a logarithmic scale, we cannot simply average data in that scale. \n",
    "    \n",
    "Please note that the **correct scale** in which operations need to be performed **is the power scale.** This is critical, e.g. when speckle filters are applied, spatial operations like block averaging are performed, or time series are analyzed.\n",
    "\n",
    "To **convert from dB to power**, apply: $\\gamma^o_{pwr} = 10^{\\frac{\\gamma^o_{dB}}{10}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calPwr = np.power(10,caldB/10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 3.3 Create a Time Series Animation\n",
    "\n",
    "Now we are ready to create a time series animation from the calibrated SAR data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create and move into a directory in which to store our plots and animations:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "product_path = analysis_dir/f'plots_and_animations'\n",
    "if not product_path.exists():\n",
    "    product_path.mkdir()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create a masked raster stack:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vv = rtc.sel(band='VV').squeeze().values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get dates of stack:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dates = [pd.to_datetime(time).strftime('%Y-%m-%d') for time in rtc.time.values]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generate a matplotlib time-series animation:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "fig = plt.figure(figsize=(14, 8))\n",
    "ax = fig.subplots()\n",
    "ax.axis('off')\n",
    "vmin = np.nanpercentile(vv.flatten(), 5)\n",
    "vmax = np.nanpercentile(vv.flatten(), 95)\n",
    "\n",
    "r0dB = 20 * np.log10(vv) - 83\n",
    "\n",
    "im = ax.imshow(vv[0], cmap='gray', vmin=vmin, vmax=vmax)\n",
    "ax.set_title(f'{dates[0]}')\n",
    "\n",
    "def animate(i):\n",
    "    ax.set_title(f'{dates[i]}')\n",
    "    im.set_data(vv[i])\n",
    "\n",
    "# Interval is given in milliseconds\n",
    "ani = animation.FuncAnimation(fig, animate, frames=vv.shape[0], interval=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Configure matplotlib's RC settings for the animation:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc('animation', embed_limit=40971520.0)  # We need to increase the limit maybe to show the entire animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create a javascript animation of the time-series running inline in the notebook:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HTML(ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Delete the dummy png** that was saved to the current working directory while generating the javascript animation in the last code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    product_path/Path('None0000000.png').unlink()\n",
    "except FileNotFoundError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Save the animation (animation.gif):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ani.save(f'{product_path}/animation.gif', writer='pillow', fps=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "### 3.4 Plot the Time Series of Means Calculated Across the Subset\n",
    "\n",
    "To create the time series of means, we will go through the following steps:\n",
    "1. Ensure that you use the data in **power scale** ($\\gamma^o_{pwr}$) for your mean calculations.\n",
    "1. compute means.\n",
    "1. convert the resulting mean values into dB scale for visualization.\n",
    "1. plot time series of means.\n",
    "\n",
    "**Compute the means:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data = rtc.sel(band='VV').squeeze().values\n",
    "rs_means_pwr = np.nanmean(input_data, axis=(1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convert resulting mean value time-series to dB scale for visualization:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs_means_dB = 10.*np.log10(rs_means_pwr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot and save the time series of means (RCSoverTime.png):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    fig = plt.figure(figsize=(16, 4))\n",
    "    ax1 = fig.subplots()\n",
    "    window_length = len(rs_means_pwr)-1\n",
    "    if window_length % 2 == 0:\n",
    "        window_length -= 1\n",
    "    polyorder = math.ceil(window_length*0.1)\n",
    "    yhat = scipy.signal.savgol_filter(rs_means_pwr, window_length, polyorder) \n",
    "    ax1.plot(dates, yhat, color='red', marker='o', markerfacecolor='white', linewidth=3, markersize=6)\n",
    "    ax1.plot(dates, rs_means_pwr, color='black', linewidth=1)\n",
    "    plt.grid()\n",
    "    ax1.set_xlabel('Date')\n",
    "    ax1.tick_params(axis='x', labelrotation=45)\n",
    "    ax1.set_ylabel(r'$\\overline{\\gamma^o}$ [power]')\n",
    "    plt.savefig(f'{product_path}/RCSoverTime.png', dpi=72, transparent='true')\n",
    "except ValueError as e:\n",
    "    print(f'Error: polyorder: {polyorder} >= window_length: {window_length}')\n",
    "    raise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 3.6 Calculate Coefficient of Variance\n",
    "\n",
    "The coefficient of variance describes how much the $\\sigma_{0}$ or $\\gamma_{0}$ measurements in a pixel vary over time. Hence, the coefficient of variance can indicate different vegetation cover and soil moisture regimes in your area.\n",
    "\n",
    "**Write a function to convert our plots into GeoTiffs:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def geotiff_from_plot(source_image, out_filename, extent, utm, cmap=None, vmin=None, vmax=None, interpolation=None, dpi=300):\n",
    "    assert \".\" not in out_filename, 'Error: Do not include the file extension in out_filename'\n",
    "    assert type(extent) == list and len(extent) == 2 and len(extent[0]) == 2 and len(\n",
    "        extent[1]) == 2, 'Error: extent must be a list in the form [[upper_left_x, upper_left_y], [lower_right_x, lower_right_y]]'\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.axis('off')\n",
    "    plt.imshow(source_image, cmap=cmap, vmin=vmin, vmax=vmax, interpolation=interpolation)\n",
    "    temp = f\"{out_filename}_temp.png\"\n",
    "    plt.savefig(temp, dpi=dpi, transparent='true', bbox_inches='tight', pad_inches=0)\n",
    "\n",
    "    cmd = f\"gdal_translate -of Gtiff -a_ullr {extent[0][0]} {extent[0][1]} {extent[1][0]} {extent[1][1]} -a_srs EPSG:{utm} {temp} {out_filename}.tiff\"\n",
    "    !{cmd}\n",
    "    try:\n",
    "        Path(temp).unlink()\n",
    "    except FileNotFoundError:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot the Coefficient of Variance Map and save it as a png (Coeffvar.png):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "input_data = rtc.sel(band='VV').squeeze().values\n",
    "input_data[np.isnan(input_data)] = 0\n",
    "\n",
    "test = np.var(input_data, axis=0) ** 0.5\n",
    "mtest = np.mean(input_data, axis=0)\n",
    "coeffvar = test/(mtest+0.001)\n",
    "\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "fig = plt.figure(figsize=(13, 10))\n",
    "ax = fig.subplots()\n",
    "ax.axis('off')\n",
    "vmin = np.percentile(coeffvar.flatten(), 1)\n",
    "vmax = np.percentile(coeffvar.flatten(), 99)\n",
    "ax.set_title('Coefficient of Variance Map')\n",
    "\n",
    "im = ax.imshow(coeffvar, cmap='jet', vmin=vmin, vmax=vmax)\n",
    "fig.colorbar(im, ax=ax, shrink=0.5)\n",
    "plt.savefig(f\"{product_path}/Coeffvar.png\", dpi=300, transparent='true')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Save the coefficient of variance map as a GeoTiff (Coeffvar.tiff):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyproj import CRS\n",
    "\n",
    "crs = CRS.from_wkt(rtc.coords['spatial_ref'].attrs['crs_wkt'])\n",
    "utm = crs.to_epsg()\n",
    "\n",
    "minx, x_size, _, maxy, _, y_size = [float(x) for x in rtc.coords['spatial_ref'].attrs['GeoTransform'].split(' ')]\n",
    "_, _, n_y, n_x = rtc.shape\n",
    "coords = [(minx, maxy), (minx + (x_size * n_x), maxy + (y_size * n_y))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "geotiff_from_plot(coeffvar, f\"{product_path}/Coeffvar\", coords, utm, cmap='jet', vmin=vmin, vmax=vmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 3.7 Threshold Coefficient of Variance Map\n",
    "\n",
    "This is an example how to threshold the derived coefficient of variance map. This can be useful, e.g., to detect areas of active agriculture.\n",
    "\n",
    "**Plot and save the coefficient of variance histogram and CDF (thresh_coeff_var_histogram.png):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 14})\n",
    "fig = plt.figure(figsize=(14, 6)) # Initialize figure with a size\n",
    "ax1 = fig.add_subplot(121)  # 121 determines: 2 rows, 2 plots, first plot\n",
    "ax2 = fig.add_subplot(122)\n",
    "# Second plot: Histogram\n",
    "# IMPORTANT: To get a histogram, we first need to *flatten* \n",
    "# the two-dimensional image into a one-dimensional vector.\n",
    "h = ax1.hist(coeffvar.flatten(), bins=200, range=(1e-6, 1.5))\n",
    "ax1.xaxis.set_label_text('Coefficient of Variation')\n",
    "ax1.set_title('Coeffvar Histogram')\n",
    "plt.grid()\n",
    "n, bins, patches = ax2.hist(coeffvar.flatten(), bins=200, range=(1e-6, 1.5), cumulative='True', density='True', histtype='step', label='Empirical')\n",
    "ax2.xaxis.set_label_text('Coefficient of Variation')\n",
    "ax2.set_title('Coeffvar CDF')\n",
    "plt.grid()\n",
    "plt.savefig(f\"{product_path}/thresh_coeff_var_histogram.png\", dpi=72, transparent='true')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot the Threshold Coefficient of Variance Map and save it as a png (Coeffvarthresh.png):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 14})\n",
    "\n",
    "outind = np.where(n > 0.90)\n",
    "threshind = np.min(outind)\n",
    "thresh = bins[threshind]\n",
    "coeffvarthresh = np.copy(coeffvar)\n",
    "coeffvarthresh[coeffvarthresh < thresh] = 0\n",
    "coeffvarthresh[coeffvarthresh >= thresh] = 1\n",
    "fig = plt.figure(figsize=(13, 10))\n",
    "ax = fig.subplots()\n",
    "ax.axis('off')\n",
    "ax.set_title(r'Thresholded Coeffvar Map $\\alpha=90\\%$')\n",
    "im = ax.imshow(coeffvarthresh, cmap='Blues', vmin=0, vmax=1, interpolation='none')\n",
    "bar = fig.colorbar(im, ax=ax, shrink=0.5)\n",
    "plt.savefig(f\"{product_path}/Coeffvarthresh.png\", dpi=300, transparent='true')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Save the Threshold Coefficient of Variance Map as a GeoTiff (Coeffvarthresh.tiff):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "geotiff_from_plot(coeffvarthresh, f\"{product_path}/Coeffvarthresh\", coords, utm, cmap='jet', vmin=vmin, vmax=vmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "osl_opera_rtc_s1_jbook_env [conda env:.local-osl_opera_rtc_s1_jbook_env]",
   "language": "python",
   "name": "conda-env-.local-osl_opera_rtc_s1_jbook_env-osl_opera_rtc_s1_jbook_env"
  },
  "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}