Aligning Sentinel 1 and 2 data in Python

This script is a tutorial for aligning, in terms of collocating, two products in different projections to the same projection. This is carried out on the NetCDF/CF versions of the Sentinel-1 GRD and Sentinel-2 L1C from In addition, this tutorial covers:

  • streaming of data using OPeNDAP (i.e no downloads involved)
  • subsetting of data by means of OPeNDAP and slicing
  • basic feature extraction using e.g. the skimage and scipy modules

The "scientific" motivation is measuring the temporal change in sea ice concentration over a particular area in a multi sensor context. Please note that this tutorial is made as show-case and thus the scientific value of the output should not be considered as optimal.

You can inspect the products using the visualization in here (Note: you have to be logged in):

To create a conda environment with all the necessary packages, use the following command:

  • conda create -n s1_s2_sic python=3.7 netCDF4=1.4.0 matplotlib cartopy numpy scikit-image scipy pyproj jupyter pyresample -c conda-forge

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

  • conda activate s1_s2_sic
  • jupyter-notebook
In [1]:
Author(s):     Trygve Halsne 27.10.2019 (
Institution:   Norwegian Meteorological Institute, 2017

Credits: Norwegian Space Agency, Copernicus, ESA

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import as ccrs
import numpy as np
%matplotlib inline
In [2]:
# Creating a function to check that the products covers the same AOI
def create_extent(lon,lat,nx,ny):
    """ A basic function extracting the boundary of a 2D raster layer"""
    swath_left=np.array([(lon[row,0],lat[row,0]) for row in range(ny)]).T
    swath_right=np.array([(lon[row,-1],lat[row,-1]) for row in range(ny)]).T
    swath_front=np.array([(lon[0,col],lat[0,col]) for col in range(nx)]).T
    swath_rear=np.array([(lon[-1,col],lat[-1,col]) for col in range(nx)]).T

    return np.concatenate((swath_left,swath_right,swath_front,swath_front,swath_rear),axis=1)

Reading input data:

  • s1: Sentinel-1
  • s2: Sentinel-2
In [3]:
s1_url = ''
s2_url = ''

ncin_s1 = Dataset(s1_url)
ncin_s2 = Dataset(s2_url)

# Subsetting parameters for extracting overlayed data
s1_indices = [5400,6400,3000,4000] #miny, maxy, minx, maxx
s2_indices = [0,1500,0,2000] #miny, maxy, minx, maxx

# Reading data of interest

s1_hv = ncin_s1['Amplitude_HV'][0,s1_indices[0]:s1_indices[1], s1_indices[2]:s1_indices[3]]
s1_hv_noise = ncin_s1['noiseCorrectionMatrix_HV'][0,s1_indices[0]:s1_indices[1], s1_indices[2]:s1_indices[3]]
s1_hv_s0 = ncin_s1['sigmaNought_HV'][0,s1_indices[0]:s1_indices[1], s1_indices[2]:s1_indices[3]]
s1_lon = ncin_s1['lon'][s1_indices[0]:s1_indices[1], s1_indices[2]:s1_indices[3]]
s1_lat = ncin_s1['lat'][s1_indices[0]:s1_indices[1], s1_indices[2]:s1_indices[3]]

s2_b3 = ncin_s2['B3'][0,s2_indices[0]:s2_indices[1], s2_indices[2]:s2_indices[3]]
s2_lon = ncin_s2['lon'][s2_indices[0]:s2_indices[1], s2_indices[2]:s2_indices[3]]
s2_lat = ncin_s2['lat'][s2_indices[0]:s2_indices[1], s2_indices[2]:s2_indices[3]]
In [4]:
# Plot lat/lon extent of subsetted products
ny,nx = s2_b3.shape
s2_extent = create_extent(s2_lon,s2_lat,nx,ny)
ny,nx = s1_lon.shape
s1_extent = create_extent(s1_lon,s1_lat,nx,ny)

plt.scatter(s1_extent[0,:],s1_extent[1,:], c='b')
plt.scatter(s2_extent[0,:],s2_extent[1,:], c='r')

Non-sophisticated sea ice masking of Sentinel-2 data

In [5]:
# Non-sophisticated sea ice masking of Sentinel-2 data
from skimage import io, filters, measure
from scipy import ndimage
import matplotlib.pyplot as plt

#s2_ice_sheets = ndimage.binary_fill_holes(s2_b3 < threshold_s2)
s2_ice_sheets = ndimage.binary_fill_holes(s2_b3 < 2500)
#plt.imshow(s2_ice_sheets, cmap='gray');

s2_ice_labels = measure.label(s2_ice_sheets)
print("Number of sea ice features: {}".format(s2_ice_labels.max()))
Number of sea ice features: 16299

Sea ice masking of Sentinel-1 data using 30Db as value

In [6]:
# Sentinel-1 sea ice masking
from skimage import exposure

s0_denoised = (abs(s1_hv)**2-s1_hv_noise)/s1_hv_s0**2
s0_speckle = ndimage.median_filter(s0_denoised,5)
s0_db = 10*np.log10(s0_speckle)
s0_db_rescaled = exposure.rescale_intensity(s0_db, in_range=(-30,-4), out_range=(0, 255)).astype(np.uint8)

Plotting the data on a map

In [7]:
fig, ax = plt.subplots(nrows=1,ncols=2,subplot_kw={'projection': ccrs.AzimuthalEquidistant()})
ax[0].pcolormesh(s1_lon,s1_lat,s0_db_rescaled, transform=ccrs.PlateCarree())
ax[1].pcolormesh(s2_lon[::4,::4],s2_lat[::4,::4],s2_b3[::4,::4], transform=ccrs.PlateCarree())

Aligning the data

In [8]:
import pyresample as pr
import pyproj

source_sw_def = pr.geometry.SwathDefinition(lons=s1_lon, lats=s1_lat)
target_proj4 = pyproj.Proj(ncin_s2['UTM_projection'].proj4_string)
target_nx = 2000
target_ny = 1500
target_extent = (ncin_s2['x'][s2_indices[2]].data ,ncin_s2['y'][s2_indices[0]].data,
                 ncin_s2['x'][s2_indices[3]].data,ncin_s2['y'][s2_indices[1]].data )

target_ad = pr.utils.get_area_def('target', 'target', 'target',
                                  target_proj4.srs, target_nx, target_ny,
s1_regridded = pr.kd_tree.resample_nearest(source_sw_def,
                       s0_db_rescaled, target_ad,
                       reduce_data=False, radius_of_influence=50,

target_proj_x_coords = target_ad.projection_x_coords
target_proj_y_coords = target_ad.projection_y_coords
/home/trygveh/anaconda3/envs/test_multi_2/lib/python3.7/site-packages/pyresample/utils/ UserWarning: 'get_area_def' has moved, import it with 'from pyresample import get_area_def'
  warnings.warn("'get_area_def' has moved, import it with 'from pyresample import get_area_def'")
In [9]:
# Extracting ice/water from Sentinel-1
s1_ice_sheets = s1_regridded<1

Plotting the aligned data

In [10]:

Plotting The Sea Ice Masks

In [11]:
In [12]:
print("Sentinel-1 sea ice concentration: {}".format(1-np.sum(s1_ice_sheets)/(s1_ice_sheets).size))
print("Sentinel-2 sea ice concentration: {}".format(1-np.sum(s2_ice_sheets)/s2_ice_sheets.size))

print("Product time difference: {} hours".format((ncin_s2['time'][0]-ncin_s1['time'][0])/(60*60)))
Sentinel-1 sea ice concentration: 0.8205586666666667
Sentinel-2 sea ice concentration: 0.5890176666666667
Product time difference: 7.217777777777778 hours