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:
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.
import url_widget as url_w
notebookUrl = url_w.URLWidget()
display(notebookUrl)
from IPython.display import Markdown
from IPython.display import display
notebookUrl = notebookUrl.value
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>'))
In this notebook we will use the following scientific libraries:
%%capture
from pathlib import Path
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
gdal.UseExceptions()
%matplotlib widget
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 opensarlab_lib as asfn
asfn.jupytertheme_matplotlib_format()
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:
path = Path.cwd()/"data_Ex2-4_S1-MadreDeDios"
if not path.exists():
path.mkdir()
We will retrieve the relevant data from an Amazon Web Service (AWS) cloud storage bucket using the following command:
time_series_path = 's3://asf-jupyter-data-west/MadreDeDios.zip'
time_series = Path(time_series_path).name
!aws --region=us-west-2 --no-sign-request s3 cp $time_series_path $time_series
Now, let's unzip the file (overwriting previous extractions) and clean up after ourselves:
if Path(time_series).exists():
asfn.asf_unzip(str(path), time_series)
Path(time_series).unlink()
Create a variable containing the VRT filename and the image acquisition dates:
!gdalbuildvrt -separate {path}/raster_stack.vrt {path}/tiffs/*_VV.tiff
image_file_VV = path/"raster_stack.vrt"
!gdalbuildvrt -separate {path}/raster_stack_VH.vrt {path}/tiffs/*_VH.tiff
image_file_VH = path/"raster_stack_VH.vrt"
Create an index of timedelta64 data with Pandas:
!ls {path}/tiffs/*_VV.tiff | sort | sed 's/[^0-9]*//g' | cut -c 4-11 > {path}/raster_stack_VV.dates
datefile_VV = path/'raster_stack_VV.dates'
dates_VV = open(str(datefile_VV)).readlines()
tindex_VV = pd.DatetimeIndex(dates_VV)
!ls {path}/tiffs/*_VH.tiff | sed 's/[^0-9]*//g' | cut -c 4-11 > {path}/raster_stack_VH.dates
datefile_VH = path/'raster_stack_VH.dates'
dates_VH = open(str(datefile_VH)).readlines()
tindex_VH = pd.DatetimeIndex(dates_VH)
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:
stindex=[]
for i in [str(datefile_VV),str(datefile_VH)]:
with open(i) as f:
sdates = f.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()
img = gdal.Open(str(image_file_VV))
band = img.GetRasterBand(1)
raster0 = band.ReadAsArray()
band_number = 0 # Needed for updates
rasterstack_VV = img.ReadAsArray()
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:
print(img.RasterCount) # Number of Bands
print(img.RasterXSize) # Number of Pixels
print(img.RasterYSize) # Number of Lines
The following line calculates the minimum backscatter per pixel across the time series:
db_mean = np.min(rasterstack_VV, axis=0)
Write a class to create an interactive plot from which we can select interesting image locations for a time series.
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')
Now we are ready to plot the minimum image. Click a point interest for which you want to analyze radar brightness over time:
# 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)
Save the selected coordinates
sarloc = (ceil(my_plot.x), ceil(my_plot.y))
print(sarloc)
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.
imagelist = [str(image_file_VV), str(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)
Now, let's pick a 5x5 image area around a center pixel defined in variable sarloc...
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)
... and extract the time series for this small area around the selected center pixel:
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 = []
Create a directory in which to store our plots and animations:
product_path = path/'plots_and_animations'
if not product_path.exists():
product_path.mkdir()
Plot the extracted time series for VV and VH polarizations:
fig, ax = plt.subplots(1, 1, figsize=(8, 4))
s_ts[0].plot(ax=ax, color='red', label='C-VV')
s_ts[1].plot(ax=ax, color='blue', label='C-VH')
ax.set_xlabel('Date')
ax.set_ylabel(r'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(product_path/f"{figname}", dpi=300, transparent='true')
Explore this data set some more by picking different point coordinates to explore. Use the time series animation together with the minimum plot to identify interesting areas and explore the radar brightness history. Discuss with your colleagues what you find.
Exercise3-ExploreSARTimeSeriesDeforestation - Version 1.5.1 - February 2024
Version Changes: