Dectect and extract the Folgefonna glacier in Sentinel-2 data

This script is a tutorial for carrying out image segmentation in Python on Copernicus Sentinel-2 data in NetCDF/CF format from satellittdata.no. This tutorials covers the topic for a user beeing interested in a specific feature in a specific part of the Sentinel-2 tile. In addition, 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
  • basic band math on optical imagery.
  • basic feature extraction by means of thresholding

The "scientific" motivation is how to automatically extract the extent of the Folgefonna glacier in optical imagery which can easily be extended to a time-series analysis.

NOTE: to run this code, you don't need to download the Sentinel data since it 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 folgefonna_ex python=3.7 netCDF4=1.4.0 matplotlib numpy scikit-image scipy geopandas jupyter fiona shapely cartopy

In your terminal, activate the environment and run the jupyter-notebook:

conda activate folgefonna_ex jupyter-notebook

In [1]:
""" 
===============================================
Name:          detect_folgefonna.ipynb
Author(s):     Trygve Halsne 26.10.2017 (dd.mm.YYYY)
Institution:   Norwegian MeteMeteorological Institute
Credits:       Norwegian Space Agency, Copernicus EU, ESA
===============================================
"""
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import sys
from scipy.spatial import distance
from skimage import measure
import geopandas as gpd
from fiona.crs import from_epsg
from shapely import geometry
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import cartopy.crs as ccrs
In [2]:
# Read product using OPeNDAP
url = "http://nbstds.met.no/thredds/dodsC/NBS/S2A/2015/08/18/S2A_OPER_PRD_MSIL1C_PDMC_20160707T072908_R094_V20150818T110049_20150818T110049/S2A_OPER_MSI_L1C_TL_EPA__20160707T045940_A000806_T32VLM_N02_04.nc"
ncin = Dataset(url, 'r')

# Decide area of interest
x_start = 4000
y_start = 3000
x_size = 2000
y_size = 3000

# Read image bands from OPeNDAP
b3_data = ncin.variables['B3'][0,y_start:y_start+y_size, x_start: x_start+x_size]
b11_data = ncin.variables['B11'][0,y_start:y_start+y_size, x_start: x_start+x_size]
lats = ncin.variables['lat'][y_start:y_start+y_size, x_start: x_start+x_size]
lons = ncin.variables['lon'][y_start:y_start+y_size, x_start: x_start+x_size]

Doing the bandmaths

In [3]:
# Computing the Normalized Difference Snow Index
NDSI = 0.0 + (b11_data.astype(np.float32) - b3_data.astype(np.float32))/(b11_data.astype(np.float32) + b3_data.astype(np.float32))*1.0

# Generate your desired figures/plots
fig2 = plt.figure(figsize=[15,15])  # a new figure window
ax2 = fig2.add_subplot(1, 1, 1)  # specify (nrows, ncols, axnum)
ax2.imshow(NDSI,cmap='gray')
ax2.set_title('NDSI: Normalised Difference Snow Index',fontsize=22)
plt.show()

Creating Polygon

In [4]:
# Find contours at a constant value of -.90
contours = measure.find_contours(NDSI, -.90)

contours_filtered = []
for c in contours:
    if c.shape[0]>2000:
        contours_filtered.append(c)


# Create an empty geopandas GeoDataFrame
newdata = gpd.GeoDataFrame()
# Create a new column called 'geometry' to the GeoDataFrame
newdata['geometry'] = None
# Set the GeoDataFrame's coordinate system to WGS84
newdata.crs = from_epsg(4326)


# Find contours with a certain distance
for n, contour in enumerate(contours_filtered):
    dists = distance.cdist(contour[::4], contour[::4], 'euclidean') # Every fourth element to save memory..
    
    if dists.max() > 600:
        coords = []
        for c in contour:
            lat = lats[int(c[0]),int(c[1])]
            lon = lons[int(c[0]),int(c[1])]
            coords.append([lon,lat])

        poly = geometry.Polygon([[c[0], c[1]] for c in coords])
        newdata.loc[n, 'geometry'] = poly

        newdata.loc[n,'idx_name'] = 'contour_' + str(n)
        break
    
# Write the data into a Shapefile
newdata.to_file("contour_folgefonna_ipy.shp")

#plt.show()
        

Plotting the polygon on a map

In [5]:
fig = plt.figure(figsize=[12,12])
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())

ax.set_extent([5, 7, 59, 61], ccrs.Geodetic())
shape_feature = ShapelyFeature(Reader('contour_folgefonna_ipy.shp').geometries(),
                               ccrs.PlateCarree(), edgecolor='black')
ax.coastlines(resolution='50m')
ax.add_feature(shape_feature, facecolor='red')
plt.show()