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:
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:
In your terminal, activate the environment and run the jupyter-notebook:
"""
=======================================================
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
# 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)
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]]
# 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
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
# 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)
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()
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'")
# Extracting ice/water from Sentinel-1
s1_ice_sheets = s1_regridded<1
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()
plt.figure();plt.imshow(s2_ice_sheets,cmap='gray')
plt.figure();plt.imshow(np.flipud(s1_ice_sheets),cmap='gray')
plt.show()
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