<jupyter_start><jupyter_text>Exploring SAR Time Series Data over Ecosystems and Deforestation Sites  <NAME>; University of Alaska Fairbanks & <NAME>, Earth Big Data, LLC  This notebook introduces you to the time series signatures over forested sites and sites affected by deforestation. The data analysis is done 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. This notebook covers the following data analysis concepts:- 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.- How to extract time series information for individual pixels of an image.- Typical time series signatures over forests and deforestation sites.  Important Notes about JupyterHub  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.<jupyter_code>%%javascript
var kernel = Jupyter.notebook.kernel;
var command = ["notebookUrl = ",
               "'", window.location, "'" ].join('')
kernel.execute(command)
from IPython.display import Markdown
from IPython.display import display

user = !echo $JUPYTERHUB_USER
env = !echo $CONDA_PREFIX
if env[0] == '':
    env[0] = 'Python 3 (base)'
if env[0] != '/home/jovyan/.local/envs/rtc_analysis':
    display(Markdown(f'<text style=color:red><strong>WARNING:</strong></text>'))
    display(Markdown(f'<text style=color:red>This notebook should be run using the "rtc_analysis" conda environment.</text>'))
    display(Markdown(f'<text style=color:red>It is currently using the "{env[0].split("/")[-1]}" environment.</text>'))
    display(Markdown(f'<text style=color:red>Select the "rtc_analysis" from the "Change Kernel" submenu of the "Kernel" menu.</text>'))
    display(Markdown(f'<text style=color:red>If the "rtc_analysis" environment is not present, use <a href="{notebookUrl.split("/user")[0]}/user/{user[0]}/notebooks/conda_environments/Create_OSL_Conda_Environments.ipynb"> Create_OSL_Conda_Environments.ipynb </a> to create it.</text>'))
    display(Markdown(f'<text style=color:red>Note that you must restart your server after creating a new environment before it is usable by notebooks.</text>'))<jupyter_output><empty_output><jupyter_text>0. Importing Relevant Python Packages  In this notebook we will use the following scientific libraries:     Pandas 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.      GDAL 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.     NumPy 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.      Matplotlib 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.<jupyter_code>%%capture
import os # for chdir, getcwd, path.basename, path.exists
from math import ceil

import pandas as pd # for DatetimeIndex
import numpy as np #for log10, mean, percentile, power
from osgeo import gdal # for GetRasterBand, Open, ReadAsArray

%matplotlib notebook
import matplotlib.patches as patches  # for Rectangle
import matplotlib.pyplot as plt # for add_subplot, axis, figure, imshow, legend, plot, set_axis_off, set_data,
                                # set_title, set_xlabel, set_ylabel, set_ylim, subplots, title, twinx
plt.rcParams.update({'font.size': 12})

import asf_notebook as asfn
asfn.jupytertheme_matplotlib_format()<jupyter_output><empty_output><jupyter_text>1. Load Data Stack    This notebook will be using a 78-image deep dual-polarization C-band SAR data stack over Madre de Dios in Peru to analyze time series signatures of vegetation covers, water bodies, and areas affected by deforestation. The C-band data were acquired by ESA's Sentinel-1 SAR sensor constellation and are available to you through the services of the Alaska Satellite Facility. The site in question is interesting as it has experienced extensive logging over the last 10 years (see image to the right; Monitoring of the Andean Amazon Project). Since the 1980s, people have been clearing forests in this area for farming, cattle ranching, logging, and (recently) gold mining. Creating RGB color composites is an easy way to visualize ongoing changes in the landscape.Before we get started, let's first create a working directory for this analysis and change into it:<jupyter_code>path = "/home/jovyan/notebooks/SAR_Training/English/Hazards/data_Ex2-4_S1-MadreDeDios"
asfn.new_directory(path)
os.chdir(path)
print(f"Current working directory: {os.getcwd()}")<jupyter_output><empty_output><jupyter_text>We will retrieve the relevant data from an Amazon Web Service (AWS) cloud storage bucket using the following command:<jupyter_code>time_series_path = 's3://asf-jupyter-data/MadreDeDios.zip'
time_series = os.path.basename(time_series_path)
!aws --region=us-east-1 --no-sign-request s3 cp $time_series_path $time_series<jupyter_output><empty_output><jupyter_text>Now, let's unzip the file (overwriting previous extractions) and clean up after ourselves:<jupyter_code>if asfn.path_exists(time_series):
    asfn.asf_unzip(os.getcwd(), time_series)
    os.remove(time_series)<jupyter_output><empty_output><jupyter_text>2. Define Data Directory and Path to VRT   Create a variable containing the VRT filename and the image acquisition dates:<jupyter_code>!gdalbuildvrt -separate raster_stack.vrt tiffs/*_VV.tiff
image_file_VV = "raster_stack.vrt"
!gdalbuildvrt -separate raster_stack_VH.vrt tiffs/*_VH.tiff
image_file_VH = "raster_stack_VH.vrt"<jupyter_output><empty_output><jupyter_text>Create an index of timedelta64 data with Pandas:<jupyter_code>!ls tiffs/*_VV.tiff | sort | cut -c 7-21 > raster_stack_VV.dates
datefile_VV = 'raster_stack_VV.dates'
dates_VV = open(datefile_VV).readlines()
tindex_VV = pd.DatetimeIndex(dates_VV)

!ls tiffs/*_VH.tiff | sort | cut -c 7-21 > raster_stack_VH.dates
datefile_VH = 'raster_stack_VH.dates'
dates_VH = open(datefile_VH).readlines()
tindex_VH = pd.DatetimeIndex(dates_VH)<jupyter_output><empty_output><jupyter_text>3. Assess Image Acquisition Dates    Before we start analyzing the available image data, we want to examine the content of our data stack. From the date index, we make and print a lookup table for band numbers and dates:<jupyter_code>stindex=[]
for i in [datefile_VV,datefile_VH]:
    sdates=open(i).readlines()
    stindex.append(pd.DatetimeIndex(sdates))
    j=1
    print('\nBands and dates for',i.strip('.dates'))
    for k in stindex[-1]:
        print("{:4d} {}".format(j, k.date()),end=' ')
        j+=1
        if j%5==1: print()<jupyter_output><empty_output><jupyter_text>4. Create Minimum Image to Identify Likely Areas of Deforestation    4.1 Load Time Series Stack  First, we load the raster stack into memory and calculate the minimum backscatter in the time series:<jupyter_code>img = gdal.Open(image_file_VV)
band = img.GetRasterBand(1)
raster0 = band.ReadAsArray()
band_number = 0 # Needed for updates
rasterstack_VV = img.ReadAsArray()<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>print(img.RasterCount) # Number of Bands
print(img.RasterXSize) # Number of Pixels
print(img.RasterYSize) # Number of Lines<jupyter_output><empty_output><jupyter_text>The following line calculates the minimum backscatter per pixel across the time series:<jupyter_code>db_mean = np.min(rasterstack_VV, axis=0)<jupyter_output><empty_output><jupyter_text>4.2 Visualize the Minimum Image and Select a Coordinate for a Time Series  Write a class to create an interactive plot from which we can select interesting image locations for a time series.<jupyter_code>class pixelPicker:
    def __init__(self, image, width, height):
        self.x = None
        self.y = None
        self.fig = plt.figure(figsize=(width, height))
        self.ax = self.fig.add_subplot(111, visible=False)
        self.rect = patches.Rectangle(
            (0.0, 0.0), width, height, 
            fill=False, clip_on=False, visible=False)
        self.rect_patch = self.ax.add_patch(self.rect)
        self.cid = self.rect_patch.figure.canvas.mpl_connect('button_press_event', 
                                                             self)
        self.image = image
        self.plot = self.gray_plot(self.image, fig=self.fig, return_ax=True)
        self.plot.set_title('Select a Point of Interest')
        
        
    def gray_plot(self, image, vmin=None, vmax=None, fig=None, return_ax=False):
        '''
        Plots an image in grayscale.
        Parameters:
        - image: 2D array of raster values
        - vmin: Minimum value for colormap
        - vmax: Maximum value for colormap
        - return_ax: Option to return plot axis
        '''
        if vmin is None:
            vmin = np.nanpercentile(self.image, 1)
        if vmax is None:
            vmax = np.nanpercentile(self.image, 99)
        #if fig is None:
        #   my_fig = plt.figure() 
        ax = fig.add_axes([0.1,0.1,0.8,0.8])
        ax.imshow(image, cmap=plt.cm.gist_gray, vmin=vmin, vmax=vmax)
        if return_ax:
            return(ax)
        
    
    def __call__(self, event):
        print('click', event)
        self.x = event.xdata
        self.y = event.ydata
        for pnt in self.plot.get_lines():
            pnt.remove()
        plt.plot(self.x, self.y, 'ro')<jupyter_output><empty_output><jupyter_text>Now we are ready to plot the minimum image. Click a point interest for which you want to analyze radar brightness over time:<jupyter_code># Large plot of multi-temporal average of VV values to inspect pixel values
fig_xsize = 7.5
fig_ysize = 7.5
my_plot = pixelPicker(db_mean, fig_xsize, fig_ysize)<jupyter_output><empty_output><jupyter_text>Save the selected coordinates:<jupyter_code>sarloc = (ceil(my_plot.x), ceil(my_plot.y))
print(sarloc)<jupyter_output><empty_output><jupyter_text>5. Plot SAR Brightness Time Series at Point Locations    5.1 SAR Brightness Time Series at Point Locations   We will pick a pixel location identified in the SAR image above and plot the time series for this identified point. By focusing on image locations undergoing deforestation, we should see the changes in the radar cross section related to the deforestation event.    First, for processing of the imagery in this notebook we generate a list of image handles and retrieve projection and georeferencing information.<jupyter_code>imagelist=[image_file_VV, image_file_VH]
geotrans=[]
proj=[]
img_handle=[]
xsize=[]
ysize=[]
bands=[]
for i in imagelist:
    img_handle.append(gdal.Open(i))
    geotrans.append(img_handle[-1].GetGeoTransform())
    proj.append(img_handle[-1].GetProjection())
    xsize.append(img_handle[-1].RasterXSize)
    ysize.append(img_handle[-1].RasterYSize)
    bands.append(img_handle[-1].RasterCount)
# for i in proj:
#     print(i)
# for i in geotrans:
#     print(i)
# for i in zip(['C-VV','C-VH','NDVI','B3','B4','B5'],bands,ysize,xsize):
#     print(i)<jupyter_output><empty_output><jupyter_text>Now, let's pick a 5x5 image area around a center pixel defined in variable sarloc...<jupyter_code>ref_x=geotrans[0][0]+sarloc[0]*geotrans[0][1]
ref_y=geotrans[0][3]+sarloc[1]*geotrans[0][5]
print('UTM Coordinates      ',ref_x, ref_y)
print('SAR pixel/line       ',sarloc[0], sarloc[1])
subset_sentinel=(sarloc[0], sarloc[1], 5, 5)<jupyter_output><empty_output><jupyter_text>... and extract the time series for this small area around the selected center pixel:<jupyter_code>s_ts=[]
for idx in (0, 1):
    means=[]
    for i in range(bands[idx]):
        rs=img_handle[idx].GetRasterBand(i+1).ReadAsArray(*subset_sentinel)
        rs_means_pwr = np.mean(rs)
        rs_means_dB = 10.*np.log10(rs_means_pwr)
        means.append(rs_means_dB)
    s_ts.append(pd.Series(means,index=stindex[idx]))
        
means = []<jupyter_output><empty_output><jupyter_text>Plot the extracted time series for VV and VH polarizations:<jupyter_code>%matplotlib inline
fig, ax = plt.subplots(1, 1, figsize=(8, 4))

s_ts[0].plot(ax=ax, color='red', label='C-VV')#,xlim=(min(min(stindex),min(stindex[0])),
                                             #     max(max(stindex),max(stindex[0]))))
s_ts[1].plot(ax=ax, color='blue', label='C-VH')
ax.set_xlabel('Date')
ax.set_ylabel('Sentinel-1 $\gamma^o$ [dB]')

ax.set_title('Sentinel-1 Backscatter')
plt.grid()
_ = ax.legend(loc='best')
_ = fig.suptitle('Time Series Profiles of Sentinel-1 SAR Backscatter')
figname = f"RCSTimeSeries-{ref_x:.0f}_{ref_y:.0f}.png"
plt.savefig(figname, dpi=300, transparent='true')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0941517163441164, 'fraction_numerical': 0.009182087865517729, 'mean_word_length': 4.54602428515472, 'pattern_counts': {'":': 0, '<': 91, '<?xml version=': 0, '>': 93, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '42884', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4520, 'n_tokens_neox': 4327, 'n_words': 1420, 'path': 'SAR_Training/English/Hazards/Exercise3B-ExploreSARTimeSeriesDeforestation.ipynb', 'repo_name': 'danielbeny/opensarlab-notebooks', 'size': '14176'}