{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "![OpenSARlab notebook banner](NotebookAddons/blackboard-banner.png)\n", "\n", "# Exploring SAR Data and SAR Time Series Analysis with Supplied Data\n", "\n", "### Franz J Meyer; University of Alaska Fairbanks & Josef Kellndorfer, [Earth Big Data, LLC](http://earthbigdata.com/)\n", "\n", "\n", "\n", "This notebook introduces 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", "**This notebook covers the following data analysis concepts:**\n", "\n", "- How to load time series stacks into Jupyter Notebooks and how to explore image content using basic functions such as mean value calculation and histogram analysis.\n", "- How to apply calibration constants to covert initial digital number (DN) data into calibrated radar cross section information.\n", "- How to subset images create time series information of calibrated SAR amplitude values.\n", "- How to explore the time-series information in SAR data stacks for environmental analysis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Important Notes about JupyterHub**\n", "\n", "Your JupyterHub server will automatically shutdown when left idle for more than 1 hour. Your notebooks will not be lost but you will have to restart their kernels and re-run them from the beginning. You will not be able to seamlessly continue running a partially run notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import url_widget as url_w\n", "notebookUrl = url_w.URLWidget()\n", "display(notebookUrl)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "from IPython.display import Markdown\n", "from IPython.display import display\n", "\n", "notebookUrl = notebookUrl.value\n", "user = !echo $JUPYTERHUB_USER\n", "env = !echo $CONDA_PREFIX\n", "if env[0] == '':\n", " env[0] = 'Python 3 (base)'\n", "if env[0] != '/home/jovyan/.local/envs/rtc_analysis':\n", " display(Markdown(f'WARNING:'))\n", " display(Markdown(f'This notebook should be run using the \"rtc_analysis\" conda environment.'))\n", " display(Markdown(f'It is currently using the \"{env[0].split(\"/\")[-1]}\" environment.'))\n", " display(Markdown(f'Select the \"rtc_analysis\" from the \"Change Kernel\" submenu of the \"Kernel\" menu.'))\n", " display(Markdown(f'If the \"rtc_analysis\" environment is not present, use Create_OSL_Conda_Environments.ipynb to create it.'))\n", " display(Markdown(f'Note that you must restart your server after creating a new environment before it is usable by notebooks.'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 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", "- [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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "from pathlib import Path\n", "\n", "import pandas as pd # for DatetimeIndex\n", "from osgeo import gdal # for GetRasterBand, Open, ReadAsArray\n", "gdal.UseExceptions()\n", "import numpy as np #for log10, mean, percentile, power\n", "\n", "%matplotlib inline\n", "import matplotlib.pylab as plb # for add_patch, add_subplot, figure, hist, imshow, set_title, xaxis,_label, text \n", "import matplotlib.pyplot as plt # for add_subplot, axis, figure, imshow, legend, plot, set_axis_off, set_data,\n", " # set_title, set_xlabel, set_ylabel, set_ylim, subplots, title, twinx\n", "import matplotlib.patches as patches # for Rectangle\n", "import matplotlib.animation as an # for FuncAnimation\n", "from matplotlib import rc \n", "\n", "import opensarlab_lib as asfn\n", "asfn.jupytertheme_matplotlib_format()\n", "\n", "from IPython.display import HTML" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 1. Load Data Stack\n", "\n", " \n", "\n", "This notebook will be using a 70-image deep L-band SAR data stack over Nepal for a first experience with time series processing. The L-band data were acquired by the ALOS PALSAR sensor and are available to us through the services of the [Alaska Satellite Facility](https://www.asf.alaska.edu/). \n", "\n", "Nepal is an interesting site for this analysis due to the significant seasonality of precipitation that is characteristic for this region. Nepal is said to have five seasons: spring, summer, monsoon, autumn and winter. Precipitation is low in the winter (November - March) and peaks dramatically in the summer, with top rain rates in July, August, and September (see figure to the right). As SAR is sensitive to changes in soil moisture, these weather patterns have a noticeable impact on the Radar Cross Section ($\\sigma$) time series information. \n", "\n", "We will analyze the variation of $\\sigma$ values over time and will interpret them in the context of rainfall rates in the imaged area. \n", "\n", "Before we get started, let's first **create a working directory for this analysis and change into it:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "path = Path(\"/home/jovyan/notebooks/SAR_Training/English/Master/data_time_series_example\")\n", "\n", "if not path.exists():\n", " path.mkdir()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will **retrieve the relevant data** from an [Amazon Web Service (AWS)](https://aws.amazon.com/) cloud storage bucket **using the following command:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s3_path = 's3://asf-jupyter-data-west/time_series.zip'\n", "time_series_path = path/Path(s3_path).name\n", "\n", "!aws --region=us-west-2 --no-sign-request s3 cp $s3_path $time_series_path" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's **unzip the file (overwriting previous extractions) and clean up after ourselves:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if time_series_path.exists():\n", " asfn.asf_unzip(str(path), str(time_series_path))\n", " time_series_path.unlink()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following lines set path variables needed for data processing. This step is not necessary but it saves a lot of extra typing later. **Define variables for the main data directory as well as for the files containing data and image information:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "datadirectory = path/'time_series/S32644X696260Y3052060sS1-EBD'\n", "datefile = datadirectory/'S32644X696260Y3052060sS1_D_vv_0092_mtfil.dates'\n", "imagefile = datadirectory/'S32644X696260Y3052060sS1_D_vv_0092_mtfil.vrt'\n", "imagefile_cross = datadirectory/'S32644X696260Y3052060sS1_D_vh_0092_mtfil.vrt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 2. Switch to the Data Directory:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!ls *.vrt #Uncomment this line to see a List of the files " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 3. Assess Image Acquisition Dates\n", "\n", "Before we start analyzing the available image data, we want to examine the content of our data stack. **First, we read the image acquisition dates for all files in the time series and create a *pandas* date index.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if datefile.exists():\n", " with open(str(datefile), 'r') as f:\n", " dates = f.readlines()\n", " tindex = pd.DatetimeIndex(dates)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the date index, we **make and print a lookup table for band numbers and dates:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if imagefile.exists():\n", " j = 1\n", " print('Bands and dates for', imagefile)\n", " for i in tindex:\n", " print(\"{:4d} {}\".format(j, i.date()),end=' ')\n", " j += 1\n", " if j%5 == 1: print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 4. Explore the Available Image Data \n", "\n", "To **open an image file using the gdal.Open() function.** This returns a variable (img) that can be used for further interactions with the file:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if imagefile.exists():\n", " img = gdal.Open(str(imagefile))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To **explore the image (number of bands, pixels, lines),** you can use several functions associated with the image object (img) created in the last code cell:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(img.RasterCount) # Number of Bands\n", "print(img.RasterXSize) # Number of Pixels\n", "print(img.RasterYSize) # Number of Lines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.1 Reading Data from an Image Band\n", "\n", "**To access any band in the image**, use GDAL's *GetRasterBand(x)* function. Replace the band_num value with the number of the band you wish to access." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "band_num = 70 \n", "band = img.GetRasterBand(band_num)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once a band is seleted, several functions associated with the band are available for further processing, e.g., *band.ReadAsArray(xoff=0,yoff=0,xsize=None,ysize=None)*\n", "\n", "**Let's read the entire raster layer for the band:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "raster = band.ReadAsArray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.2 Extracting Subsets from a Larger Image Frame\n", "\n", "Because of the potentially large data volume when dealing with time series data stacks, it may be prudent to read only a subset of data. \n", "\n", "Using GDAL's *ReadAsArray()* function, subsets can be requested by defining pixel offsets and subset size:\n", "\n", "**img.ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None)**\n", "\n", "- *xoff, yoff* are the offsets from the upper left corner in pixel/line coordinates. \n", "- *xsize, ysize* specify the size of the subset in x-direction (left to right) and y-direction (top to bottom).\n", "\n", "For example, we can **read only a subset of 5x5 pixels with an offset of 5 pixels and 20 lines:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "raster_sub = band.ReadAsArray(5, 20, 50, 50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is a two dimensional numpy array in the datatpye the data were stored in. **We can inspect these data in python by typing the array name on the commandline**:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "raster_sub" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3 Displaying Bands in the Time Series of SAR Data\n", "\n", " From the lookup table we know that bands 20 and 27 in the Nepal data stack are from mid February and late August. **Let's take look at these images**. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "raster_1 = img.GetRasterBand(20).ReadAsArray()\n", "raster_2 = img.GetRasterBand(27).ReadAsArray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**4.3.1 Write a Plotting Function**\n", "\n", "Matplotlib's plotting functions allow for powerful options to display imagery. We are following some standard approaches for setting up figures.\n", "First we are looking at a **raster band** and it's associated **histogram**.\n", "\n", "Our function, *show_image()* takes several parameters:\n", " \n", "- raster = a numpy two dimensional array \n", "- tindex = a panda index array for dates\n", "- bandnbr = the band number the corresponds to the raster \n", "- vmin = minimim value to display \n", "- vmax = maximum value to display\n", "- output_filename = name of output file, if saving the plot\n", "\n", "Note: By default, data will be linearly stretched between vmin and vmax.\n", "\n", "*We won't use this function in this notebook but it is a useful utility method, which can be copied and pasted for use in other analyses*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def show_image_histogram(raster, tindex, band_nbr, vmin=None, vmax=None, output_filename=None):\n", " assert 'plb' in globals(), 'Error: matplotlib.pylab must be imported as \"plb\"'\n", " assert 'plt' in globals(), 'Error: matplotlib.pyplot must be imported as \"plt\"' \n", " \n", " fig = plb.figure(figsize=(16, 8))\n", " ax1 = fig.add_subplot(121)\n", " ax2 = fig.add_subplot(122)\n", " \n", " # plot image\n", " ax1.imshow(raster, cmap='gray', vmin=vmin, vmax=vmax)\n", " ax1.set_title('Image Band {} {}'.format(band_nbr, tindex[band_nbr-1].date()))\n", " vmin = np.percentile(raster, 2) if vmin==None else vmin\n", " vmax = np.percentile(raster, 98) if vmax==None else vmax\n", " ax1.xaxis.set_label_text('Linear stretch Min={} Max={}'.format(vmin, vmax))\n", " \n", " #plot histogram\n", " h = ax2.hist(raster.flatten(), bins=200, range=(0, 10000))\n", " ax2.xaxis.set_label_text('Amplitude (Uncalibrated DN Values)')\n", " ax2.set_title('Histogram Band {} {}'.format(band_nbr, tindex[band_nbr-1].date()))\n", " \n", " if output_filename:\n", " plt.savefig(f'{output_filename}', dpi=300, transparent='true')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**We won't be calling our new function elsewhere in this notebook, so test it now:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# change the last argument to decide where you want to store image. \n", "show_image_histogram(raster_1, tindex, 20, vmin=2000, vmax=10000, output_filename=path/f'show_time_series.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 5. SAR Time Series Visualization, Animation, and Analysis\n", "\n", "This section introduces you to the handling and analysis of SAR time series stacks. A focus will be put on time series visualization, which allow us to inspect time series in more depth. Note that html animations are not exported into the pdf file, but will display interactively.\n", "\n", "### 5.1 Reading the SAR Time Series Subset\n", "\n", "Let's read an image subset (offset 400, 400 / size 600, 600) of the entire time series data stack. The data are linearly scaled amplitudes represented as unsigned 16 bit integer.\n", "\n", "We use the GDAL *ReadAsArray(xoff,yoff,xsize,ysize)* function where *xoff* is the offset in pixels from upper left; *yoff* is the offset in lines from upper left; *xsize* is the number of pixels and *ysize* is the number of lines of the subset.\n", "\n", "If *ReadAsArray()* is called without any parameters, the entire image data stack is read. \n", "\n", "Let's first **define a subset and make sure it is in the right geographic location.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Open the image and read the first raster band\n", "band = img.GetRasterBand(1)\n", "\n", "# Define the subset\n", "subset = (400, 400, 600, 600)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to **extract this subset from all slices of the data stack.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot one band together with the outline of the selected subset to verify its geographic location.\n", "raster = band.ReadAsArray()\n", "vmin = np.percentile(raster.flatten(), 5)\n", "vmax = np.percentile(raster.flatten(), 95)\n", "fig = plb.figure(figsize=(10, 10))\n", "ax = fig.add_subplot(111)\n", "ax.imshow(raster, cmap='gray', vmin=vmin, vmax=vmax)\n", "# plot the subset as rectangle\n", "_ = ax.add_patch(patches.Rectangle((subset[0], subset[1]), subset[2], subset[3], fill=False, edgecolor='red'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "raster0 = band.ReadAsArray(*subset)\n", "bandnbr = 0 # Needed for updates\n", "rasterstack = img.ReadAsArray(*subset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Close img, as it is no longer needed in the notebook:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "img = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.2 Calibration and Data Conversion between dB and Power Scales\n", "\n", "Focused SAR image data natively come in uncalibrated digital numbers (DN) and need to be calibrated to 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", "Let's **apply the calibration constant for our data and export it in *dB* scale:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "caldB = 20*np.log10(rasterstack) - 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": [ "### 5.3 Create a Time Series Animation\n", "\n", "First, **Create a directory in which to store our plots and move into it:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "product_path = path/'plots_and_animations'\n", "if not product_path.exists():\n", " product_path.mkdir()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to **create a time series animation** from the calibrated SAR data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture \n", "fig = plt.figure(figsize=(10, 10))\n", "ax = fig.add_subplot(111)\n", "ax.axis('off')\n", "vmin = np.percentile(caldB.flatten(), 5)\n", "vmax = np.percentile(caldB.flatten(), 95)\n", "r0dB = 20*np.log10(raster0) - 83\n", "im = ax.imshow(r0dB,cmap='gray', vmin=vmin, vmax=vmax)\n", "ax.set_title(\"{}\".format(tindex[0].date()))\n", "\n", "def animate(i):\n", " ax.set_title(\"{}\".format(tindex[i].date()))\n", " im.set_data(caldB[i])\n", "\n", "# Interval is given in milliseconds\n", "ani = an.FuncAnimation(fig, animate, frames=rasterstack.shape[0], interval=400)" ] }, { "cell_type": "markdown", "metadata": {}, "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": [ "**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": {}, "source": [ "### 5.3 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. Compute means using the data in **power scale** ($\\gamma^o_{pwr}$) .\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": [ "rs_means_pwr = np.mean(calPwr,axis=(1, 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Convert the resulting mean value time-series to dB scale for visualization and check that we got the means over time:**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rs_means_dB = 10.*np.log10(rs_means_pwr)\n", "rs_means_pwr.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Plot and save the time series of means (time_series_means.png):**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# 3. Now let's plot the time series of means\n", "fig = plt.figure(figsize=(16, 4))\n", "ax1 = fig.add_subplot(111)\n", "ax1.plot(tindex, rs_means_pwr)\n", "ax1.set_xlabel('Date')\n", "ax1.set_ylabel(r'$\\overline{\\gamma^o}$ [power]')\n", "\n", "\n", "ax2 = ax1.twinx()\n", "ax2.plot(tindex, rs_means_dB, color='red')\n", "ax2.set_ylabel(r'$\\overline{\\gamma^o}$ [dB]')\n", "fig.legend(['power', 'dB'], loc=1)\n", "plt.title(r'Time series profile of average band backscatter $\\gamma^o$ ')\n", "plt.savefig(f'{product_path}/time_series_means', dpi=72, transparent='true')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.4 Create Two-Panel Figure with Animated Global Mean $\\mu_{\\gamma^0_{dB}}$\n", "\n", "We use a few Matplotlib functions to **create a side-by-side animation of the dB-scaled imagery and the respective global means $\\mu_{\\gamma^0_{dB}}$.**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 4), gridspec_kw={'width_ratios':[1, 3]})\n", "\n", "vmin = np.percentile(rasterstack.flatten(), 5)\n", "vmax = np.percentile(rasterstack.flatten(), 95)\n", "im = ax1.imshow(raster0, cmap='gray', vmin=vmin, vmax=vmax)\n", "ax1.set_title(\"{}\".format(tindex[0].date()))\n", "ax1.set_axis_off()\n", "\n", "ax2.axis([tindex[0].date(), tindex[-1].date(), rs_means_dB.min(), rs_means_dB.max()])\n", "ax2.set_ylabel('$\\overline{\\gamma^o}$ [dB]')\n", "ax2.set_xlabel('Date')\n", "ax2.set_ylim((-10, -5))\n", "l, = ax2.plot([], [])\n", "\n", "def animate(i):\n", " ax1.set_title(\"{}\".format(tindex[i].date()))\n", " im.set_data(rasterstack[i])\n", " ax2.set_title(\"{}\".format(tindex[i].date()))\n", " l.set_data(tindex[:(i+1)], rs_means_dB[:(i+1)])\n", "\n", "# Interval is given in milliseconds\n", "ani = an.FuncAnimation(fig, animate, frames=rasterstack.shape[0], interval=400)" ] }, { "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": [ "**Save the animated time-series and histogram (animation_histogram.gif):**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ani.save(f'{product_path}/animation_histogram.gif', writer='pillow', fps=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Time_Series_Example.ipynb - Version 1.5.3 - February 2024*\n", "\n", "*Version Changes*\n", "\n", "- *Use raw strings for LateX in matplotlib labels*" ] } ], "metadata": { "kernelspec": { "display_name": "rtc_analysis", "language": "python", "name": "conda-env-.local-rtc_analysis-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.12.1" } }, "nbformat": 4, "nbformat_minor": 4 }