In [1]:
""" A test script for monitoring vegetation change over a predefined area of
    interest by means of Sentinel-2 data from the Norwegian National Ground
    Segment (NBS) read through OPeNDAP.

    Made as a show case using Python 2.7

===========================================================
Name:          monitor_change_sentinel_2
Author(s):     Trygve Halsne 19.09.2018 (dd.mm.YYYY)
Modifications:
Copyright:     (c) Norwegian Meteorological Institute, 2018
===========================================================
"""


from netCDF4 import Dataset
import matplotlib.pyplot as plt
import datetime as dt
import numpy as np
from skimage import exposure

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.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.iteritems():
        output_name = str(composite + '_' + str(i) + '.png')
        save_plot(output_name, img)
    ncin.close()
Read all variables in  0:00:09.292439
Read all variables in  0:00:10.468528