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 satellittdata.no. 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 satellittdata.no here (Note: you have to be logged in): https://satellittdata.no/en/metsis/map/wms?dataset=S1B_EW_GRDM_1SDH_20190827T052357_20190827T052457_017767_0216F4_E6FA%2CS2A_MSIL1C_20190827T123701_N0208_R095_T39XVL_20190827T143509&solr_core=nbs-l1&calling_results_page=basket

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 (dd.mm.YYYY)
Institution:   Norwegian Meteorological Institute, 2017

Credits: Norwegian Space Agency, Copernicus, ESA
=======================================================
"""

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import cartopy.crs 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 = 'http://nbstds.met.no/thredds/dodsC/NBS/S1B/2019/08/27/EW/S1B_EW_GRDM_1SDH_20190827T052357_20190827T052457_017767_0216F4_E6FA.nc'
s2_url = 'http://nbstds.met.no/thredds/dodsC/NBS/S2A/2019/08/27/S2A_MSIL1C_20190827T123701_N0208_R095_T39XVL_20190827T143509.nc'

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')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.show()

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');plt.show()

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_denoised[s0_denoised<=0]=10e-7
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())
plt.show()

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,
                                  target_extent)
s1_regridded = pr.kd_tree.resample_nearest(source_sw_def,
                       s0_db_rescaled, target_ad,
                       reduce_data=False, radius_of_influence=50,
                       fill_value=None)

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/__init__.py:30: 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]:
plt.figure();plt.imshow(s2_b3,cmap='gray',
                        extent=[target_extent[0],target_extent[2],target_extent[3],target_extent[1]])
#plt.xticks(np.arange(target_extent[0],target_extent[2],10))
plt.figure();plt.imshow(np.flipud(s1_regridded),cmap='gray',
                        extent=[target_proj_x_coords[0],target_proj_x_coords[-1],
                                target_proj_y_coords[0],target_proj_y_coords[-1]])

plt.show()

Plotting The Sea Ice Masks

In [11]:
plt.figure();plt.imshow(s2_ice_sheets,cmap='gray')
plt.figure();plt.imshow(np.flipud(s1_ice_sheets),cmap='gray')
                        
plt.show()
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