As part of the ILB8 session of the OTC23 online training course, this Jupyter notebook will use OVL portal to retrieve inputs data

author: dr.fab, OceanDataLab
date: 2022-10-24

Computation of SAR image spectra¶

Sentinel-1 sea surface roughness data are modulated by long ocean waves

When zoomed to the maximum level, these modulations are clearly visible in the SAR images and can be used to retrieve informations on long wave properties

This notebook will teach you how to:

Import data from OVL
Compute spectra in cartesian and log-polar coordinates

Sentinel1 sea surface roughness data will be extracted from the OVL portal as demonstrated in our “OVL extraction How To” YouTube video.

import os
import numpy
import numpy as np
from numpy.fft import fft2, fftshift, ifft2, fftfreq
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.signal import medfilt

1. Define mathematical functions to compute spectra¶

def looks2xspec(im1,im2,periodo_size):
    imshape = numpy.array(im1.shape, dtype='int32')

    ###########################################################################
    # Set periodograms/looks/specs sizes and positions
    ###########################################################################

    aziwindow = numpy.hanning(periodo_size+2)[1:-1]
    ranwindow = numpy.hanning(periodo_size+2)[1:-1]
    window = numpy.sqrt(numpy.outer(aziwindow, ranwindow))
    ###########################################################################
    # Compute periodograms and compute co/cross spectra
    ###########################################################################

    count = numpy.ceil(imshape[0]/periodo_size).astype('int32')
    perpos = (numpy.floor(numpy.linspace(0, imshape[0]-periodo_size, num=count)+0.5).astype('int32'),
              numpy.floor(numpy.linspace(0, imshape[1]-periodo_size, num=count)+0.5).astype('int32'))
    specshape = numpy.array((periodo_size, periodo_size), dtype='int32')
    specs = numpy.zeros(specshape, dtype='complex64')

    for appos in iter(perpos[0]):
        for rppos in iter(perpos[1]):
            sub1 = im1[appos:appos+periodo_size, rppos:rppos+periodo_size]
            sub1 = sub1 - numpy.mean(sub1)
            per1 = fftshift(fft2(sub1*window))/periodo_size
            sub2 = im2[appos:appos+periodo_size, rppos:rppos+periodo_size]
            sub2 = sub2 - numpy.mean(sub2)
            per2 = fftshift(fft2(sub2*window))/periodo_size
            specs += per1 * numpy.conj(per2)
            
    specs[int(periodo_size/2-2):int(periodo_size/2+3), int(periodo_size/2-2):int(periodo_size/2+3)] = 0
    specs[int(periodo_size/2-3):int(periodo_size/2+4), int(periodo_size/2)] = 0
    specs[int(periodo_size/2), int(periodo_size/2-3):int(periodo_size/2+4)] = 0

    return specs

2. Extract data from OVL portal¶

Choose the downloaded .npy file and set the image diagonal length in km (as displayed in the OVL portal)¶

you may change the periodio_size according to the spectral resolution needed (64, 128, 256, 512, ...)¶

sarname = '/home/drfab/Téléchargements/3857_Sentinel-1A_SAR_roughness-s1a-iw-grd-vv-20221016t182149-20221016t182214-045471-057009-001-CA83.npy'
diagonal_km = 25.56

sigo = numpy.flipud(np.load(sarname))

fig = plt.figure(figsize=(16,16))
plt.imshow(numpy.flipud(sigo),interpolation='bicubic',cmap='gray')

3. Compute and plot spectra in cartesian coordinates¶

ground_spacing = diagonal_km *1000 / np.sqrt(sigo.shape[1]**2+sigo.shape[0]**2)

periodo_size = 256

spec = looks2xspec(sigo, sigo,periodo_size)
kran = (numpy.arange(periodo_size)-periodo_size/2.)/periodo_size*numpy.pi/ground_spacing*2
kazi = (numpy.arange(periodo_size)-periodo_size/2.)/periodo_size*numpy.pi/ground_spacing*2
fig = plt.figure(figsize = (20,16))
cut = int(periodo_size/3)
plt.imshow(numpy.flipud(numpy.abs(spec[cut:-cut+1,cut:-cut-1])*1e-5), extent=[kran[0+cut], kran[-1-cut],
            kazi[-1-cut], kazi[1+cut]], aspect='auto',cmap='jet')
kcircle=[2*np.pi/800, 2*np.pi/400, 2*np.pi/200, 2*np.pi/100]
theta = 2. * np.pi * np.linspace(0, 1, num=361)
for kcir in kcircle:
    plt.plot(kcir*np.cos(theta), kcir*np.sin(theta), ':w')
    kcirstr = '%im' % (np.round(2*np.pi/kcir))
    plt.text(0, -(kcir+.005), kcirstr, ha='center', va='top',fontsize=None,color='white')

plt.xlabel('wavenumber [rad/m]')
plt.ylabel('wavenumber [rad/m]')
plt.title('Absolute value')
plt.colorbar()

4. Convert and plot spectra in Log polar coordinates¶

from scipy.interpolate import RectBivariateSpline
def cart2logpol(spec, kazi, kran, nphi=72, nk=60, kmin=2.*np.pi/1200,
                kmax=2.*np.pi/30, heading=None):
    """
    """
    # Define log-polar grid
    phi = np.linspace(0., 360., num=nphi, endpoint=False)
    dphi = phi[1] - phi[0]
    alpha = (kmax / kmin) ** (1. / (nk - 1.))
    k = kmin * alpha ** np.arange(nk)
    dk = (np.sqrt(alpha) - 1. / np.sqrt(alpha)) * k
    # Bilinear interpolation
    intfunc = RectBivariateSpline(kazi, kran, spec, kx=1, ky=1)
    if heading is None:
        # phi clockwise, phi=0 means up azimuth
        angincart = np.deg2rad(90. - phi)
    else:
        # (heading is expected clockwise from north)
        # phi clockwise, phi=0 means north
        angincart = np.deg2rad(90. - (phi - heading))
    intkazi = k[np.newaxis, :] * np.sin(angincart[:, np.newaxis])
    intkran = k[np.newaxis, :] * np.cos(angincart[:, np.newaxis])
    polspec = np.zeros((nphi, nk), dtype=spec.dtype)
    indint = np.where((intkazi >= kazi[0]) & (intkazi <= kazi[-1]) & \
                      (intkran >= kran[0]) & (intkran <= kran[-1]))
    polspec[indint] = intfunc(intkazi[indint], intkran[indint], grid=False)
    #intspec = intfunc(intkazi, intkran, grid=False)
    # Energy conservation
    dkazi = kazi[1] - kazi[0]
    dkran = kran[1] - kran[0]
    kcart = np.sqrt(kazi[:, np.newaxis] ** 2. + kran[np.newaxis, :] ** 2.)
    indcart = np.where((kcart >= kmin) & (kcart <= kmax))
    enecart = 4. * np.sqrt(np.sum(np.abs(spec[indcart])) * dkazi * dkran)
    areapol = k * dk * np.deg2rad(dphi)
    enepol = 4. * np.sqrt(np.sum(np.abs(polspec) * areapol[np.newaxis, :]))
    #print enecart, enepol, (enecart / enepol) ** 2.
    polspec *= (enecart / enepol) ** 2.
    return polspec, phi, dphi, k, dk

polspec, phi, dphi, ks, dk = cart2logpol(numpy.abs(spec)*1e-5, kazi, kran,nphi=72*2)

s_phigeo = phi
_dphigeo = dphi
_pl = np.concatenate(([s_phigeo[0] - _dphigeo / 2.], s_phigeo[:-1] + _dphigeo / 2., [s_phigeo[-1] + _dphigeo / 2.]))
_pl = np.deg2rad(90. - _pl)
lmax = 100
indk = np.where((ks <= 2. * np.pi / lmax) & (ks >= 2. * np.pi / 800))[0]
_kl = np.concatenate(([0], ks[indk] + dk[indk] / 2.))
_xl = _kl[np.newaxis, :] * np.cos(_pl[:, np.newaxis])
_yl = _kl[np.newaxis, :] * np.sin(_pl[:, np.newaxis])
klims = [-2. * np.pi / lmax, 2. * np.pi / lmax]
klims2 = [-2. * np.pi / lmax *0.7, 2. * np.pi / lmax*0.7]

fig = plt.figure(figsize=(20,16))
plt.subplots_adjust(left=0.15, right=0.975, bottom=0.1, top=0.95, wspace=0.145)

plt.pcolormesh(_xl, _yl, medfilt(polspec[:, indk],3), cmap=plt.get_cmap('jet'))#, vmin=vmin, vmax=vmax)
plt.plot([0, 0], klims, 'silver',linewidth=0.6)
plt.plot(klims, [0, 0], 'silver',linewidth=0.6)
plt.plot(klims2, klims2, 'silver',linewidth=0.6)
plt.plot(klims2[::-1], klims2, 'silver',linewidth=0.6)
# plt.pcolormesh(_xl, _yl, polspec[:, indk], cmap=plt.get_cmap('jet'))#, vmin=vmin, vmax=vmax)
cirx = np.cos(np.linspace(0, 2. * np.pi, num=100))
ciry = np.sin(np.linspace(0, 2. * np.pi, num=100))
for wl in [400, 200, 100]:
    plt.plot(2. * np.pi / wl * cirx, 2. * np.pi / wl * ciry,color='silver',linewidth=0.6 )
    plt.text(0, 2. * np.pi / wl,'{:d}m'.format(wl),color='k',horizontalalignment='center')
plt.text(0, 2. * np.pi / 95,'North',color='k',horizontalalignment='center')
plt.xlabel('Eastward wavenumber (rad/m)')
plt.ylabel('Northward wavenumber (rad/m)')
#plt.axis('off')
plt.gca().set_aspect('equal')
plt.colorbar()