Observing Etna Eruption¶
This use case was prepared with support from European Space Agency.
1 Introduction¶
In early spring 2017, biggest Sicilian volcano - Etna - awakened and erupted several times. The Etna eruption started on January 23rd 2017 but intensified after February 27th and the activity continued till late April.
In this exercise we will demonstrate how such events can be observed using satellite images. We will utilize imagery acquired from satellites Sentinel 2 (s2) and Landsat 8 (l8). We will introduce the basic physical concepts that we need in order to understand the example but we will not go into theoretical details.
First, we need to set up two "technical" things:
The text and code you see here are written in Jupyter Notebook. You will need to "execute cells" to see satellite images. You can do that by clicking on the cell with code and then pressing "ctrl + Enter".
We will use Sentinel Hub to download the data. A provisional access is already provided for your convenience (i.e. a long number assigned to INSTANCE_ID variable) but it will support only limited number of requests to avoid abuse. If you want to play around and request more data, you will need to create your own account and replace the value of INSTANCE_ID with your own one. Check instructions to do so here.
Right below is a first cell that needs to be executed. Click on it and press: "ctrl + Enter".
from sentinelhub import SHConfig
INSTANCE_ID = 'cf8091e5-2837-463e-b45d-8dc0b121cad7'
if INSTANCE_ID:
config = SHConfig()
config.instance_id = INSTANCE_ID
else:
config = None
# Set options
%reload_ext autoreload
%autoreload 2
%matplotlib inline
A text "In [*]" before a the cell indicates that the cell is being executed. When the "*" turns into a number, e.g. "In [3]" it means that the cell was successfully executed and you can see the results of execution below the cell. Some cells do not return any results but just set up things that will be used latter.The execution of some cells might take a while, which is also indicated with an “hourglass icon” in the name of your browser's tab:
Back to Etna now. We will download satellite data from the Sentinel Hub using sentinelhub-py. To prepare a request for downloading, we follow the example in sentinelhub-py documentation, where one can find the descriptions of all parameters. We adjusted the values so that we will download all "true color" images of Etna acquired between February 27th 2017 and April 18th 2017. The coordinates of bounding box of the area of interest (AOI) can be read from any online map, which displays the coordinates of current position of mouse (and defines coordinate system to which they refer). We used EO Browser.
from sentinelhub import (WmsRequest, WcsRequest, MimeType, CRS, BBox,
DataSource, CustomUrlParam, geo_utils, DataCollection)
import datetime
import numpy as np
import requests as r
etna_bbox_wgs84 = [14.96, 37.70, 15.06, 37.78]
etna_bbox = BBox(bbox=etna_bbox_wgs84, crs=CRS.WGS84)
true_color_request_s2 = WmsRequest(data_folder='data',
data_collection=DataCollection.SENTINEL2_L1C,
layer='TRUE-COLOR-S2-L1C',
bbox=etna_bbox,
time=('2017-02-27','2017-04-18'),
width=512,
maxcc=0.8,
config=config,
custom_url_params={
CustomUrlParam.SHOWLOGO: False})
true_color_img_s2 = true_color_request_s2.get_data(save_data=True)
true_color_dates_s2 = true_color_request_s2.get_dates()
print('We downloaded %d images, each one acquired on a different date (i.e. \
time series). ' % len(true_color_img_s2))
print('Single image is of type {} and has a shape {}. First two dimensions \
indicate the size of downloaded image in pixels and the last is a number of downloaded\
bands.'.format(type(true_color_img_s2[-1]), true_color_img_s2[-1].shape))
We downloaded 11 images, each one acquired on a different date (i.e. time series). Single image is of type <class 'numpy.ndarray'> and has a shape (516, 512, 4). First two dimensions indicate the size of downloaded image in pixels and the last is a number of downloaded bands.
Important: If the execution of the previous cell is not finished (i.e. the "In [*]" still shows "*") you won’t be able to continue. You might need to wait a while.
To have a means to display the downloaded images appropriately, we will prepare a function "plot_image", which will add labels and scale bar to the image. Run the cell below to display one of the images (we selected the 4th image in the time series, since it shows smoke and snow. You could choose any other image by changing the value of "display_image" variable).
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt
def plot_image(image, bbox, factor=1):
"""
Utility function for plotting RGB images.
:param image: image to be plotted
:type image: numpy array
:param bbox: bounding box object, which will be used to set the
extent of the plot
:type bbox: sentinelhub.common.bbox
:param factor: factor
:type factor: int
"""
fig = plt.subplots(nrows=1, ncols=1, figsize=(9, 9))
bbox_polygon = bbox.get_polygon()
if np.issubdtype(image.dtype, np.floating):
plt.imshow(np.minimum(image * factor, 1), extent=[
bbox_polygon[0][0], bbox_polygon[2][0],
bbox_polygon[0][1], bbox_polygon[2][1]])
else:
plt.imshow(image, extent=[bbox_polygon[0][0],
bbox_polygon[2][0],
bbox_polygon[0][1],
bbox_polygon[2][1]])
#Add axes' lables
if bbox.crs == CRS.WGS84 :
plt.xlabel('longitude [°]')
plt.ylabel('latitude [°]')
else:
plt.xlabel('easting [m]')
plt.ylabel('northing [m]')
#Add scalebar
width = image.shape[1]
height = image.shape[2]
resx, resy = geo_utils.bbox_to_resolution(bbox,
width,
height)
scalebar = ScaleBar((width*resx)/(bbox_polygon[2][0]
-bbox_polygon[0][0]),
'm',
location='lower right',
box_alpha=0.9)
plt.gca().add_artist(scalebar)
display_image = 4
plot_image(true_color_img_s2[display_image], etna_bbox)
plt.title('Etna, true color image acquired with Sentinel 2 on '+ str(
true_color_dates_s2[display_image].date()));
print()