#!/usr/bin/env python
# coding: utf-8

# # Looking at temporal change in Sentinel-2 data
# 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: 
# 
# - streaming of data using OPeNDAP (i.e no downloads involved)
# - subsetting of data by means of selecting a set of variables within a confined area of interest
# - creating an RGB in Python
# - computing NDVI in Python
# 
# 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

# In[14]:


""" 
===========================================================
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
get_ipython().run_line_magic('matplotlib', 'inline')


# ## Defining some functions

# In[18]:


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')


# ## Reading and storing images to disk

# In[19]:


# 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()