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
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
"""
===============================================
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
# 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]
# 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()
# 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()
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()