This script is a tutorial for looking at temporal change using Python and Copernicus Sentinel-2 data in NetCDF/CF format from satellittdata.no. The tutorial covers:
The "scientific" motivation is how to monitor temporal change in a confined region from computing a classical false color vegetation RGB composite and Normalized Difference Vegetation Index. Please note that this tutorial is made as show-case and thus the scientific value of the output should not be considered as optimal.
NOTE: to run this code, you don't need to download the Sentinel data prior to running the script since the data is provided in the script directly by means of data streaming using OPeNDAP. However, some python packages needs to be installed. To create a conda environment with all the necessary packages, use the following command:
conda create -n s2_temp_change python=3.7 netCDF4=1.4.0 matplotlib numpy scikit-image jupyter In your terminal, activate the environment and run the jupyter-notebook:
conda activate s2_temp_change jupyter-notebook
"""
===========================================================
Name: monitor_change_sentinel_2
Author(s): Trygve Halsne 19.09.2018 (dd.mm.YYYY)
Institution: Norwegian Meteorological Institute
Credits: Norwegian Space Agency, Copernicus EU, ESA
===========================================================
"""
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import datetime as dt
import numpy as np
from skimage import exposure
%matplotlib inline
def linearStretch(image, pLow, pHigh):
"""Operator for linear stretching of histogram."""
tmp_image = []
for i in range(image.shape[np.argmin(image.shape)]):
band = image[:,:,i]
iMin, iMax = np.percentile(band[~np.isnan(band)], (pLow, 100-pHigh))
band_rescaled = exposure.rescale_intensity(band, in_range=(iMin, iMax))
tmp_image.append(band_rescaled)
img_rescale = np.dstack(tmp_image)
return img_rescale
def get_uint8_image( scaled_image, vmin, vmax):
''' Scale image from float (or any) input array to uint8
Parameters
----------
image : 2D matrix
vmin : float - minimum value
vmax : float - maximum value
Returns
-------
2D matrix
'''
# redistribute into range [0,255]
uint8Image = 255 * ((scaled_image - vmin) / float((vmax - vmin)))
uint8Image[uint8Image < 0] = 0
uint8Image[uint8Image > 255] = 255
uint8Image[~np.isfinite(uint8Image)] = 0
return uint8Image.astype('uint8')
def save_plot(output_name, data):
''' Stores plot on local machine '''
fig = plt.figure(frameon=False)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax.imshow(data)
ax.set_aspect('auto')
ax.text(0.025, 0.95, date,
verticalalignment='center', fontsize=20,
transform=ax.transAxes, color='white',
bbox={'facecolor':'black','alpha':0.5, 'edgecolor':'gray', 'boxstyle':'round'})
fig.savefig(output_name)
plt.show()
plt.close('all')
# OPeNDAP urls for reading data
urls = ['http://nbstds.met.no/thredds/dodsC/NBS/S2A/2015/08/15/S2A_OPER_PRD_MSIL1C_PDMC_20160705T144923_R051_V20150815T105041_20150815T105041/S2A_OPER_MSI_L1C_TL_EPA__20160705T124729_A000763_T32VPM_N02_04.nc',
'http://nbstds.met.no/thredds/dodsC/NBS/S2A/2015/08/22/S2A_OPER_PRD_MSIL1C_PDMC_20160706T181324_R008_V20150822T104035_20150822T104035/S2A_OPER_MSI_L1C_TL_EPA__20160706T123006_A000863_T32VPM_N02_04.nc']
for i,url in enumerate(urls):
t0 = dt.datetime.now()
ncin = Dataset(url)
x0 = 1100; x1 = 2100; y0 = 2000; y1 = 3000 # Area covering Gardermoen
b3_data = ncin.variables['B3'][0,y0:y1,x0:x1]
b4_data = ncin.variables['B4'][0,y0:y1,x0:x1]
b8_data = ncin.variables['B8'][0,y0:y1,x0:x1]
print("Read all variables in ", dt.datetime.now() - t0)
datetime = ncin.getncattr('start_time')
date = datetime.split('T')[0]
rgb_false_color_veg = linearStretch(np.dstack([b8_data,b4_data,b3_data]),0.1,0.1)
NDVI = 0.0 + (b8_data.astype(np.float32) - b4_data.astype(np.float32))/(b8_data.astype(np.float32) + b4_data.astype(np.float32))*1.0
imagery = {'fcv': rgb_false_color_veg,'ndvi': NDVI}
# Saving plots to file
for composite,img in imagery.items():
output_name = str(composite + '_' + str(i) + '.png')
save_plot(output_name, img)
ncin.close()
Read all variables in 0:00:11.399908
Read all variables in 0:00:10.495714