Download this notebook from github.

---

Create HH LUT

This notebook will show how we create HH LUTs

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

# optional debug messages
import logging
logging.basicConfig()
logging.getLogger('xsarsea.windspeed').setLevel(logging.DEBUG) #or .setLevel(logging.INFO)

Definition Zhang & Mouche

inc_angle = np.arange(17,50,0.5)
ar = [1.3794, -3.19e-2, 1.4e-3]
br = [-0.1711, 2.6e-3]


def get_pol_ratio_zhangA(inc_angle,wind_speed):
    # do not force the ratio to be greater than 1
    ars2 = np.polynomial.polynomial.polyval(inc_angle, ar)
    brs2 = np.polynomial.polynomial.polyval(inc_angle, br)

    pol_ratio = ars2 * (wind_speed ** brs2)
    return pol_ratio

def get_pol_ratio_zhangB(inc_angle,wind_speed):
    # do force the ratio to be greater than 1
    ars2 = np.polynomial.polynomial.polyval(inc_angle, ar)
    brs2 = np.polynomial.polynomial.polyval(inc_angle, br)

    pol_ratio = ars2 * (wind_speed ** brs2)
    # we force the ratio to be greater than 1
    pol_ratio = np.where(pol_ratio < 1, 1, pol_ratio)
    return pol_ratio


def get_pol_ratio_mouche(inc_angle, wind_dir, wind_speed=None):

    theta=inc_angle
    phi=wind_dir
    # Alexis Mouche, D. Hauser,
    # V. Kudryavtsev and JF. Daloze,
    # "Multi polarization ocean radar
    # cross-section from ENVISAT ASAR
    # observations, airborne polarimetric
    # radar measurements and empirical or
    # semi-empirical models", ESA
    # ERS/ENVISAT Symposium, Salzburg,
    # September 2004
    A0 = 0.00650704
    B0 = 0.128983
    C0 = 0.992839
    Api2 = 0.00782194
    Bpi2 = 0.121405
    Cpi2 = 0.992839
    Api = 0.00598416
    Bpi = 0.140952
    Cpi = 0.992885

    P0_theta = A0*np.exp(B0*theta)+C0
    Ppi2_theta = Api2*np.exp(Bpi2*theta)+Cpi2
    Ppi_theta = Api*np.exp(Bpi*theta)+Cpi

    C0_theta = (P0_theta+Ppi_theta+2*Ppi2_theta)/4
    C1_theta = (P0_theta-Ppi_theta)/2
    C2_theta = (P0_theta+Ppi_theta-2*Ppi2_theta)/4

    polrat = C0_theta + C1_theta*np.cos(np.deg2rad(phi)) +  C2_theta*np.cos(2*np.deg2rad(phi))

    return polrat

Plot

plt.figure(figsize=(10,5))

for wspd in [3,7,10,15,20]:
    plt.plot(inc_angle, 10*np.log10(get_pol_ratio_zhangB(inc_angle,wspd)), label=f'Wspd = {wspd} m/s')

plt.legend(loc="lower right")
plt.title(f'ZhangB Polarization Ratio vs Incidence Angle')
plt.xlabel('Incidence Angle [deg]')
plt.ylabel('Polarization Ratio [dB]')
plt.xlim([17,50])
plt.grid()


plt.figure(figsize=(10,5))

for phi in [30, 60, 90, 120, 150, 180]:
    plt.plot(inc_angle, 10*np.log10(get_pol_ratio_mouche(inc_angle,phi)), label=f'Phi = {phi} deg')

plt.legend(loc="upper left")
plt.title(f'Mouche1 Polarization Ratio vs Incidence Angle')
plt.xlabel('Incidence Angle [deg]')
plt.ylabel('Polarization Ratio [dB]')
plt.xlim([17,50])
plt.grid()

create HH LUT

def create_gmfHH(fct, vv_name, pr_name, mod, res):
    pol_ratio_data = xr.apply_ufunc(
        fct,
        mod.incidence, mod.wspd,
        vectorize=True
    )
    pol_ratio_data_extended = pol_ratio_data.expand_dims({'phi': mod.phi}).broadcast_like(mod)

    mod_hh = (mod/pol_ratio_data_extended)
    mod_hh.name = 'sigma0_model'
    mod_hh = mod_hh.to_dataset()
    mod_hh.attrs['units'] = 'linear'
    mod_hh.attrs['construction'] = f'{vv_name} / {pr_name}'
    mod_hh.attrs['description'] = f'Backscatter coefficient in HH polarization build from {vv_name.upper()} model and {pr_name[0].upper() + pr_name[1:]} Polarization Ratio model'
    mod_hh.attrs['resolution'] = f'{res}'
    mod_hh.attrs['pol'] = 'HH'


    if vv_name == "cmod7":

        mod_hh.attrs['inc_range'] = np.array([16.,66.])
        mod_hh.attrs['wspd_range'] = np.array([0.2,50.])
        mod_hh.attrs['phi_range'] = np.array([0., 180.])

    elif vv_name == "cmod5n":

        mod_hh.attrs['inc_range'] = np.array([17.,50.])
        mod_hh.attrs['wspd_range'] = np.array([0.2, 50.])
        mod_hh.attrs['phi_range'] = np.array([0., 180.])

    else :
        raise ValueError("vv_name must be cmod7 or cmod5n")

    mod_hh.attrs['model'] = f'{vv_name}_R{res}_hh_{pr_name}'


    mod_hh.attrs['wspd_step'] = np.round(
        np.unique(np.diff(mod_hh.wspd)), decimals=2)[0]
    mod_hh.attrs['inc_step'] = np.round(
        np.unique(np.diff(mod_hh.incidence)), decimals=2)[0]
    mod_hh.attrs['phi_step'] = np.round(
            np.unique(np.diff(mod_hh.phi)), decimals=2)[0]
    fname = f'./nc_lut_gmf_{vv_name}_R{res}_hh_{pr_name}.nc'
    mod_hh.to_netcdf(fname, mode="w")
    print("model saved at ", fname)
    mod_hh.close()

vv_name = "cmod5n"
#create_gmfHH(get_pol_ratio_zhangB, vv_name, "zhangB", xsarsea.windspeed.get_model(f"gmf_{vv_name}").to_lut(**{'resolution':'high'}), 'high')