Merge pull request #79 from vincelhx/master

updating documentation

README.rst

@@ -46,4 +46,4 @@ ________________
 Documentation and examples
 --------------------------
 
-https://cyclobs.ifremer.fr/static/sarwing_datarmor/xsarsea
+https://cerweb.ifremer.fr/datarmor/doc_sphinx/xsarsea/

docs/examples/create_hh_lut.ipynb

@@ -0,0 +1,234 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Create HH LUT\n",
+    "\n",
+    "This notebook will show how we create HH LUTs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import xarray as xr\n",
+    "\n",
+    "# optional debug messages\n",
+    "import logging\n",
+    "logging.basicConfig()\n",
+    "logging.getLogger('xsarsea.windspeed').setLevel(logging.DEBUG) #or .setLevel(logging.INFO)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Definition Zhang & Mouche "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "inc_angle = np.arange(17,50,0.5)\n",
+    "ar = [1.3794, -3.19e-2, 1.4e-3]\n",
+    "br = [-0.1711, 2.6e-3]\n",
+    "    \n",
+    "\n",
+    "def get_pol_ratio_zhangA(inc_angle,wind_speed):\n",
+    "    # do not force the ratio to be greater than 1\n",
+    "    ars2 = np.polynomial.polynomial.polyval(inc_angle, ar)\n",
+    "    brs2 = np.polynomial.polynomial.polyval(inc_angle, br)\n",
+    "\n",
+    "    pol_ratio = ars2 * (wind_speed ** brs2)\n",
+    "    return pol_ratio\n",
+    "\n",
+    "def get_pol_ratio_zhangB(inc_angle,wind_speed):\n",
+    "    # do force the ratio to be greater than 1\n",
+    "    ars2 = np.polynomial.polynomial.polyval(inc_angle, ar)\n",
+    "    brs2 = np.polynomial.polynomial.polyval(inc_angle, br)\n",
+    "\n",
+    "    pol_ratio = ars2 * (wind_speed ** brs2)\n",
+    "    # we force the ratio to be greater than 1\n",
+    "    pol_ratio = np.where(pol_ratio < 1, 1, pol_ratio)\n",
+    "    return pol_ratio\n",
+    "\n",
+    "\n",
+    "def get_pol_ratio_mouche(inc_angle, wind_dir, wind_speed=None):\n",
+    "\n",
+    "    theta=inc_angle\n",
+    "    phi=wind_dir\n",
+    "    # Alexis Mouche, D. Hauser,\n",
+    "    # V. Kudryavtsev and JF. Daloze,\n",
+    "    # \"Multi polarization ocean radar\n",
+    "    # cross-section from ENVISAT ASAR\n",
+    "    # observations, airborne polarimetric\n",
+    "    # radar measurements and empirical or\n",
+    "    # semi-empirical models\", ESA\n",
+    "    # ERS/ENVISAT Symposium, Salzburg,\n",
+    "    # September 2004\n",
+    "    A0 = 0.00650704\n",
+    "    B0 = 0.128983\n",
+    "    C0 = 0.992839\n",
+    "    Api2 = 0.00782194\n",
+    "    Bpi2 = 0.121405\n",
+    "    Cpi2 = 0.992839\n",
+    "    Api = 0.00598416\n",
+    "    Bpi = 0.140952\n",
+    "    Cpi = 0.992885\n",
+    "\n",
+    "    P0_theta = A0*np.exp(B0*theta)+C0\n",
+    "    Ppi2_theta = Api2*np.exp(Bpi2*theta)+Cpi2\n",
+    "    Ppi_theta = Api*np.exp(Bpi*theta)+Cpi\n",
+    "\n",
+    "    C0_theta = (P0_theta+Ppi_theta+2*Ppi2_theta)/4\n",
+    "    C1_theta = (P0_theta-Ppi_theta)/2\n",
+    "    C2_theta = (P0_theta+Ppi_theta-2*Ppi2_theta)/4\n",
+    "\n",
+    "    polrat = C0_theta + C1_theta*np.cos(np.deg2rad(phi)) +  C2_theta*np.cos(2*np.deg2rad(phi))\n",
+    "\n",
+    "    return polrat\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure(figsize=(10,5))\n",
+    "\n",
+    "for wspd in [3,7,10,15,20]:\n",
+    "    plt.plot(inc_angle, 10*np.log10(get_pol_ratio_zhangB(inc_angle,wspd)), label=f'Wspd = {wspd} m/s')\n",
+    "    \n",
+    "plt.legend(loc=\"lower right\")\n",
+    "plt.title(f'ZhangB Polarization Ratio vs Incidence Angle')\n",
+    "plt.xlabel('Incidence Angle [deg]')\n",
+    "plt.ylabel('Polarization Ratio [dB]')\n",
+    "plt.xlim([17,50])\n",
+    "plt.grid()\n",
+    "\n",
+    "\n",
+    "plt.figure(figsize=(10,5))\n",
+    "\n",
+    "for phi in [30, 60, 90, 120, 150, 180]:\n",
+    "    plt.plot(inc_angle, 10*np.log10(get_pol_ratio_mouche(inc_angle,phi)), label=f'Phi = {phi} deg')\n",
+    "    \n",
+    "plt.legend(loc=\"upper left\")\n",
+    "plt.title(f'Mouche1 Polarization Ratio vs Incidence Angle')\n",
+    "plt.xlabel('Incidence Angle [deg]')\n",
+    "plt.ylabel('Polarization Ratio [dB]')\n",
+    "plt.xlim([17,50])\n",
+    "plt.grid()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## create HH LUT "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def create_gmfHH(fct, vv_name, pr_name, mod, res):\n",
+    "    pol_ratio_data = xr.apply_ufunc(\n",
+    "        fct,\n",
+    "        mod.incidence, mod.wspd,\n",
+    "        vectorize=True\n",
+    "    )\n",
+    "    pol_ratio_data_extended = pol_ratio_data.expand_dims({'phi': mod.phi}).broadcast_like(mod)\n",
+    "\n",
+    "    mod_hh = (mod/pol_ratio_data_extended)\n",
+    "    mod_hh.name = 'sigma0_model'\n",
+    "    mod_hh = mod_hh.to_dataset()\n",
+    "    mod_hh.attrs['units'] = 'linear'\n",
+    "    mod_hh.attrs['construction'] = f'{vv_name} / {pr_name}'\n",
+    "    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'\n",
+    "    mod_hh.attrs['resolution'] = f'{res}'\n",
+    "    mod_hh.attrs['pol'] = 'HH'\n",
+    "    \n",
+    "    \n",
+    "    if vv_name == \"cmod7\":\n",
+    "        \n",
+    "        mod_hh.attrs['inc_range'] = np.array([16.,66.])\n",
+    "        mod_hh.attrs['wspd_range'] = np.array([0.2,50.])\n",
+    "        mod_hh.attrs['phi_range'] = np.array([0., 180.])\n",
+    "    \n",
+    "    elif vv_name == \"cmod5n\":\n",
+    "        \n",
+    "        mod_hh.attrs['inc_range'] = np.array([17.,50.])\n",
+    "        mod_hh.attrs['wspd_range'] = np.array([0.2, 50.])\n",
+    "        mod_hh.attrs['phi_range'] = np.array([0., 180.])\n",
+    "    \n",
+    "    else : \n",
+    "        raise ValueError(\"vv_name must be cmod7 or cmod5n\")\n",
+    "    \n",
+    "    mod_hh.attrs['model'] = f'{vv_name}_R{res}_hh_{pr_name}'\n",
+    "\n",
+    "\n",
+    "    mod_hh.attrs['wspd_step'] = np.round(\n",
+    "        np.unique(np.diff(mod_hh.wspd)), decimals=2)[0]\n",
+    "    mod_hh.attrs['inc_step'] = np.round(\n",
+    "        np.unique(np.diff(mod_hh.incidence)), decimals=2)[0]\n",
+    "    mod_hh.attrs['phi_step'] = np.round(\n",
+    "            np.unique(np.diff(mod_hh.phi)), decimals=2)[0]\n",
+    "    fname = f'./nc_lut_gmf_{vv_name}_R{res}_hh_{pr_name}.nc'\n",
+    "    mod_hh.to_netcdf(fname, mode=\"w\")\n",
+    "    print(\"model saved at \", fname)\n",
+    "    mod_hh.close()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vv_name = \"cmod5n\"\n",
+    "#create_gmfHH(get_pol_ratio_zhangB, vv_name, \"zhangB\", xsarsea.windspeed.get_model(f\"gmf_{vv_name}\").to_lut(**{'resolution':'high'}), 'high')\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}