diff --git a/README.rst b/README.rst index 0e897d5..f677960 100644 --- a/README.rst +++ b/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/ \ No newline at end of file diff --git a/docs/examples/create_hh_lut.ipynb b/docs/examples/create_hh_lut.ipynb new file mode 100644 index 0000000..5406c34 --- /dev/null +++ b/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 +} diff --git a/docs/examples/gmfs_and_luts.ipynb b/docs/examples/gmfs_and_luts.ipynb index 808ed09..dc89f26 100644 --- a/docs/examples/gmfs_and_luts.ipynb +++ b/docs/examples/gmfs_and_luts.ipynb @@ -24,6 +24,7 @@ "source": [ "import xsar\n", "from xsarsea import windspeed\n", + "import xsarsea\n", "import xarray as xr\n", "import numpy as np\n", "import holoviews as hv\n", @@ -50,7 +51,8 @@ "source": [ "## Available models\n", "\n", - "Available models can be retrieved with [xsarsea.windspeed.available_models](../basic_api.rst#xsarsea.windspeed.available_models)" + "Available models can be retrieved with [xsarsea.windspeed.available_models](../basic_api.rst#xsarsea.windspeed.available_models).\n", + "By default, analytical models are already in the dataframe. " ] }, { @@ -68,19 +70,42 @@ "id": "f4999b03", "metadata": {}, "source": [ - "## Available models\n", + "## Add all models\n", "\n", - "Adding gmfs_impl (analytical models) with [xsarsea.windspeed.gmfs.GmfModel.activate_gmfs_impl](../basic_api.rst#xsarsea.windspeed.gmfs.GmfModel.activate)" + "replace paths by your own path containing all luts\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "0bdf23bd", - "metadata": {}, + "id": "583fccc9-3b9a-404c-9823-0f4c99aa0d8b", + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "windspeed.gmfs.GmfModel.activate_gmfs_impl()" + "nc_luts_path = xsarsea.get_test_file('nc_luts_reduce')\n", + "path_cmod7 = xsarsea.get_test_file(\"cmod7_and_python_script\")\n", + "windspeed.register_luts(nc_luts_path, path_cmod7)\n", + "windspeed.available_models()" + ] + }, + { + "cell_type": "markdown", + "id": "b1a91ca0-6246-48fa-b6d0-10fe4815298f", + "metadata": {}, + "source": [ + "## Add models by type \n" + ] + }, + { + "cell_type": "markdown", + "id": "b79464cf", + "metadata": {}, + "source": [ + "### Adding analytical models\n", + "\n", + "they are already available in the dataframe by default " ] }, { @@ -100,17 +125,17 @@ "metadata": {}, "outputs": [], "source": [ - "nc_luts_path = xsar.get_test_file('nc_luts_reduce')\n", + "nc_luts_path = xsarsea.get_test_file('nc_luts_reduce')\n", "windspeed.register_nc_luts(nc_luts_path)\n", "windspeed.available_models()" ] }, { "cell_type": "markdown", - "id": "4118ad0e-11a7-4d6e-970e-50a371af44f3", + "id": "856cf5ce", "metadata": {}, "source": [ - "## Addind CMOD7\n" + "### Adding cmod7" ] }, { @@ -127,32 +152,6 @@ " print(e)" ] }, - { - "cell_type": "markdown", - "id": "2a4321c1-e185-47a6-bb0d-e4f27dcea819", - "metadata": {}, - "source": [ - "### Adding Sarwing models (LUT)\n", - "\n", - "Sarwing model are not available by default, because they needs to be loaded from external file with [xsarsea.windspeed.register_sarwing_luts](../basic_api.rst#xsarsea.windspeed.register_sarwing_luts)\n", - "\n", - "A basic subset of sarwing lut can be retrieved with `xsar.get_test_file('sarwing_luts_subset')`. \n", - "\n", - "To get full sarwing luts, download them from https://cyclobs.ifremer.fr/static/sarwing_datarmor/xsardata/sarwing_luts, or use path `/home/datawork-cersat-public/cache/project/sarwing/GMFS/v1.6` (ifremer only)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "85e9e0fd-6b5f-4fd7-a74e-5af3a56bcec6", - "metadata": {}, - "outputs": [], - "source": [ - "sarwing_luts_subset_path = xsar.get_test_file('sarwing_luts_subset')\n", - "windspeed.register_sarwing_luts(sarwing_luts_subset_path)\n", - "windspeed.available_models()" - ] - }, { "cell_type": "markdown", "id": "e6b10527-aa1f-427c-b7a3-37c1b1441c28", @@ -254,7 +253,7 @@ "id": "9df945ab", "metadata": {}, "source": [ - "It won't have the same impact on SarwingLutModel or [NcLutModel](../basic_api.rst#xsarsea.windspeed.models.NcLutModel) or [cmod7Model](../basic_api.rst#xsarsea.windspeed.gmfs.cmod7Model).\n", + "It won't have the same impact on [NcLutModel](../basic_api.rst#xsarsea.windspeed.models.NcLutModel) or [cmod7Model](../basic_api.rst#xsarsea.windspeed.gmfs.cmod7Model).\n", "\n", "Indeed, these are saved at at desired format with a certain resolution. \n", "\n", @@ -269,7 +268,7 @@ "outputs": [], "source": [ "# here we specify the exact same params than the saved LUT have\n", - "windspeed.get_model('nc_lut_sarwing_lut_cmodms1ahw').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.1, 'resolution' : 'high'})" + "windspeed.get_model('nc_lut_cmodms1ahw').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.1, 'resolution' : 'high'})" ] }, { @@ -280,7 +279,7 @@ "outputs": [], "source": [ "# here we specify different params than the saved LUT have : interpolation is made\n", - "windspeed.get_model('nc_lut_sarwing_lut_cmodms1ahw').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.2, 'resolution' : 'high'})" + "windspeed.get_model('nc_lut_cmodms1ahw').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.2, 'resolution' : 'high'})" ] }, { @@ -288,7 +287,7 @@ "id": "2dd55b26", "metadata": {}, "source": [ - "## In practice\n", + "## In practice for wind retrieval\n", "[xsarsea.windspeed.invert_from_model](../basic_api.rst#xsarsea.windspeed.invert_from_model) can be called with **kwargs.\n", "\n", "We can use kwargs to force the use of high resolution Luts. \n", @@ -351,6 +350,15 @@ " return sig" ] }, + { + "cell_type": "markdown", + "id": "301d9d62", + "metadata": {}, + "source": [ + "Note : \n", + " if defer = True, you will have to use again xsarsea.windspeed.gmfs.GmfModel.activate_gmfs_impl() to activate the new model" + ] + }, { "cell_type": "code", "execution_count": null, @@ -411,8 +419,8 @@ "# loading low resolution & interpolating \n", "windspeed.get_model('nc_lut_gmf_cmod7_Rlow_hh_mouche1').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.1, 'resolution' : 'high'})\n", "# or \n", - "# loading high resolution\n", - "windspeed.get_model('nc_lut_gmf_cmod7_Rhigh_hh_mouche1').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.1, 'resolution' : 'high'})" + "# loading high resolution [not computed in the doc cause nc_lut_gmf_cmod7_Rhigh_hh_mouche1 is too big]\n", + "#windspeed.get_model('nc_lut_gmf_cmod7_Rhigh_hh_mouche1').to_lut(**{'wspd_step' : 0.1, 'phi_step' : 1, 'inc_step' : 0.1, 'resolution' : 'high'})" ] }, { @@ -478,7 +486,7 @@ "metadata": {}, "outputs": [], "source": [ - "model_compare([ 'gmf_dummy', 'nc_lut_sarwing_lut_cmodms1ahw'] )" + "model_compare([ 'gmf_dummy', 'nc_lut_cmodms1ahw'] )" ] }, { diff --git a/docs/examples/streaks.ipynb b/docs/examples/streaks.ipynb index 62cdd0b..1baafe4 100644 --- a/docs/examples/streaks.ipynb +++ b/docs/examples/streaks.ipynb @@ -361,14 +361,6 @@ "\n", "gridspace.opts(plot_size=(200,200)) + (hv_img * hv_vf).opts(legend_position='right', frame_height=500)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a762fdcb-ed24-4917-8004-d6ffee51ed25", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/docs/examples/windspeed_inversion.ipynb b/docs/examples/windspeed_inversion.ipynb deleted file mode 100644 index f3c63b6..0000000 --- a/docs/examples/windspeed_inversion.ipynb +++ /dev/null @@ -1,566 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "e7588254-89c4-4ad0-94d4-df629f2eb94c", - "metadata": {}, - "source": [ - "# Wind speed invertion\n", - "\n", - "This notebook example show how to invert wind speed from sigma0, using models from LUT or GMF, using [xsarsea.windspeed.invert_from_model](../basic_api.rst#xsarsea.windspeed.invert_from_model)" - ] - }, - { - "cell_type": "markdown", - "id": "983cebef-95b2-48f1-9708-aa9ad762b787", - "metadata": {}, - "source": [ - "> .. warning::\n", - " **Use of ancillary wind**\n", - "\n", - "> We suggest to go from `ancillary_wind = -np.conj(sarwing_ds.owi_ancillary_wind)` to `ancillary_wind = sarwing_ds.owi_ancillary_wind` ; then it won't match sarwing results anymore\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7136c9bf-eeba-42bc-b255-08cdc5710752", - "metadata": {}, - "outputs": [], - "source": [ - "import xsarsea\n", - "from xsarsea import windspeed\n", - "import xarray as xr\n", - "import numpy as np\n", - "import holoviews as hv\n", - "hv.extension('bokeh')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "77f0e14b-9f80-4fc2-89e7-19113b6977d3", - "metadata": {}, - "outputs": [], - "source": [ - "# optional debug messages\n", - "import logging\n", - "logging.basicConfig()\n", - "logging.getLogger('xsarsea.windspeed').setLevel(logging.DEBUG) # or .setLevel(logging.INFO)" - ] - }, - { - "cell_type": "markdown", - "id": "19887be3-890a-4071-bb96-07814f9bd474", - "metadata": {}, - "source": [ - "## Read sarwing owi file\n", - "\n", - "download more from https://cyclobs.ifremer.fr/static/sarwing_datarmor/processings/c39e79a/default/reports/shoc_dailyupdate_names/report.html (in the 'safe' column, download files named like `*-owi-xx-*.nc`)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "77630d90-573a-428a-9029-1c3272289879", - "metadata": {}, - "outputs": [], - "source": [ - "#sarwing_owi_file = xsarsea.get_test_file('s1a-iw-owi-xx-20210909t130650-20210909t130715-039605-04AE83.nc')\n", - "sarwing_owi_file = \"/home/vincelhx/Documents/autoentreprise/IFREMER/s1a-iw-owi-xx-20210909t130650-20210909t130715-039605-04AE83.nc\"\n", - "sarwing_ds = xsarsea.read_sarwing_owi(sarwing_owi_file)\n", - "sarwing_ds" - ] - }, - { - "cell_type": "markdown", - "id": "4ced9944-f9ba-4fa8-aa95-4e14c49fe0f8", - "metadata": {}, - "source": [ - "## Get ancillary wind\n", - "\n", - "Ecmwf wind is stored in owi file in *geographical* (deg N/S) convention. `xsarsea.windspeed` need it relative to *sample* (ie antenna), as a complex number.\n", - "\n", - "We use [xsarsea.dir_meteo_to_sample](../basic_api.rst#xsarsea.dir_meteo_to_sample) to convert `sarwing_ds.owiEcmwfWindDirection` (deg) to radians, relative to sample, using `sarwing_ds.owiHeading`\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a1e03453-5cc6-4395-a998-a9a31224e634", - "metadata": {}, - "outputs": [], - "source": [ - "owi_ecmwf_wind = sarwing_ds.owiEcmwfWindSpeed * np.exp(1j* xsarsea.dir_meteo_to_sample(sarwing_ds.owiEcmwfWindDirection, sarwing_ds.owiHeading))\n", - "sarwing_ds = xr.merge([\n", - " sarwing_ds,\n", - " owi_ecmwf_wind.to_dataset(name='owi_ancillary_wind'),\n", - "])" - ] - }, - { - "cell_type": "markdown", - "id": "a88e3200-dd0d-44e0-9760-a9bd9b2a9581", - "metadata": {}, - "source": [ - "## Ancillary wind control plot\n", - "\n", - "### Directions in default convention\n", - "\n", - "We check that ancillary wind is correct, by plotting the wind speed, the wind direction, and the wind direction components.\n", - "\n", - "The *sample* component increase with *sample* axis, and *line* component increase with *line*.\n", - "\n", - "Note that we have to use `kdims=['sample','line']` to properly display the vectorfield, because memory order `['sample','line']` is geographycally transposed " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3cb363c1-f250-4156-ad1c-f257dbcbf768", - "metadata": {}, - "outputs": [], - "source": [ - "sub_sarwing_ds = sarwing_ds.isel(line=slice(None, None, 10), sample=slice(None, None, 10))\n", - "\n", - "vectorfield = hv.VectorField(\n", - " (\n", - " sub_sarwing_ds.sample, sub_sarwing_ds.line,\n", - " xsarsea.dir_meteo_to_sample(sub_sarwing_ds.owiEcmwfWindDirection, sub_sarwing_ds.owiHeading),\n", - " sub_sarwing_ds.owiEcmwfWindSpeed\n", - " )\n", - ")\n", - "\n", - "owi_ecmwf_wind_sample = np.real(owi_ecmwf_wind)\n", - "owi_ecmwf_wind_line = np.imag(owi_ecmwf_wind)\n", - "(\n", - " hv.Image(sarwing_ds.owiEcmwfWindSpeed, kdims=['sample','line']).opts(title='speed and dir', clim=(0,50), cmap='jet') * vectorfield\n", - " + hv.Image(owi_ecmwf_wind_sample, kdims=['sample','line']).opts(title='sample component',cmap='bwr') * vectorfield\n", - " + hv.Image(owi_ecmwf_wind_line, kdims=['sample','line']).opts(title='line component',cmap='bwr') * vectorfield\n", - ").opts(title='Direction component in standard convention')" - ] - }, - { - "cell_type": "markdown", - "id": "41d3e477-a0fc-43cb-b356-927339873f8f", - "metadata": {}, - "source": [ - "### Directions in gmf or lut convention\n", - "\n", - "In gmf or lut, the sample wind direction component is **negative** if the wind vector is in the same direction as increasing sample, and line component is **positive** if if the wind vector is in the same direction as increasing line.\n", - "\n", - "Switching to one convention to another can be done with `-np.conj(complex_dir)`\n", - "\n", - "Note that in `xsarsea.windspeed`, all directions are in **gmf or lut convention**" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6d7ae283-21cd-4e72-9466-7ca0a21279fa", - "metadata": {}, - "outputs": [], - "source": [ - "(\n", - " hv.Image(sarwing_ds.owiEcmwfWindSpeed, kdims=['sample','line']).opts(title='speed and dir', clim=(0,50), cmap='jet') * vectorfield\n", - " + hv.Image(-owi_ecmwf_wind_sample, kdims=['sample','line']).opts(title='sample component',cmap='bwr') * vectorfield\n", - " + hv.Image(owi_ecmwf_wind_line, kdims=['sample','line']).opts(title='line component',cmap='bwr') * vectorfield\n", - ").opts(title='Direction component in gmf or lut convention')" - ] - }, - { - "cell_type": "markdown", - "id": "7e98f14a-90e1-4cd4-9d8b-1faf7fa11ca8", - "metadata": {}, - "source": [ - "### luts resolution params :\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8b74ab2f-4413-4fa7-ad30-9df5d73bba97", - "metadata": {}, - "outputs": [], - "source": [ - "kwargs = {\"wspd_step\" : 0.1, \"inc_step\" : 0.1, \"phi_step\" : 1.0, \"resolution\": \"high\"}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cd1d2343-1c8b-4ca3-922d-77589e4c1200", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed.gmfs.GmfModel.activate_gmfs_impl()" - ] - }, - { - "cell_type": "markdown", - "id": "b4128cac-1f6c-4db0-a3ed-4908fa3ec7ee", - "metadata": {}, - "source": [ - "## Copol inversion\n", - "\n", - "Copol wind inversion, using model and ancillary wind" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "499ad90a-2d90-4a16-a016-cb3d55c7d23d", - "metadata": {}, - "outputs": [], - "source": [ - "wind_co = windspeed.invert_from_model(\n", - " sarwing_ds.owiIncidenceAngle,\n", - " sarwing_ds.owiNrcs,\n", - " ancillary_wind = -np.conj(sarwing_ds.owi_ancillary_wind),\n", - " model='gmf_cmod5n',\n", - " **kwargs)\n", - "windspeed_co = np.abs(wind_co)" - ] - }, - { - "cell_type": "markdown", - "id": "089780fa-aeab-4b28-b39c-55f5c221e205", - "metadata": {}, - "source": [ - "## differences with sarwing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "607e8ed3-41ce-49af-8e5c-be6a128933de", - "metadata": {}, - "outputs": [], - "source": [ - "xsarsea.windspeed.available_models()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "81eaaaa8-8b07-4669-8b2e-2afd0c310e35", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed_diff =np.abs(windspeed_co) - sarwing_ds.owiWindSpeed_Tab_copol\n", - "(\n", - " (\n", - " hv.Image(np.abs(windspeed_co), kdims=['sample','line']).opts(title='xsarsea') \n", - " + hv.Image(sarwing_ds.owiWindSpeed_Tab_copol, kdims=['sample','line']).opts(title='sarwing' )\n", - " ).opts(hv.opts.Image(cmap='jet', clim=(0,50))) \n", - " + hv.Image(windspeed_diff, kdims=['sample','line']).opts(cmap='jet', clim=(-0.005,0.005), title='xsarsea-sarwing\\nmean=%.4f std=%.4f' % (np.mean(windspeed_diff), np.std(windspeed_diff)))\n", - ").opts(hv.opts.Image(colorbar=True, tools=['hover']))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "42e17d14-9c82-49d2-af81-b4077bdaeb07", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"max diff:\",np.nanmax(np.abs(windspeed_diff)),\"mean diff:\", np.mean(windspeed_diff).values,\"std : \", np.std(windspeed_diff).values)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4322b65b", - "metadata": {}, - "outputs": [], - "source": [ - "def normalize_angle(diff):\n", - " return (diff + 180) % 360 - 180\n", - "\n", - "sarwing_ds[\"winddir_co\"] = (90 - (np.angle(-np.conj(wind_co),deg=True)) + sarwing_ds.owiHeading)%360\n", - "sub_sarwing_ds = sarwing_ds.isel(line=slice(None, None, 10), sample=slice(None, None, 10))\n", - "\n", - "angle_diff = normalize_angle(sarwing_ds[\"winddir_co\"] - sarwing_ds.owiWindDirection_Tab_copol)\n", - "\n", - "(\n", - " (hv.Image(sarwing_ds[\"winddir_co\"]).opts(title='xsarsea') + \n", - " hv.Image(sarwing_ds.owiWindDirection_Tab_copol).opts(title='sarwing' )).opts(hv.opts.Image(cmap='jet', clim=(0,360))) + \n", - " hv.Image(normalize_angle(angle_diff)).opts(cmap='jet', clim=(-5,5), title='xsarsea-sarwing\\nmean=%.4f std=%.4f' % (np.mean(angle_diff), np.std(angle_diff)))\n", - ").opts(hv.opts.Image(colorbar=True, tools=['hover']))" - ] - }, - { - "cell_type": "markdown", - "id": "20a4a4cc-d209-4062-961f-0a6420c0fba4", - "metadata": {}, - "source": [ - "## Crosspol inversion\n", - "\n", - "Sarwing crosspol inversion is done with model `cmodms1ahw`, that is only available from lut.\n", - "\n", - "So we use [xsarsea.windspeed.register_nc_luts](../basic_api.rst#xsarsea.windspeed.register_all_nc_luts) to register this lut.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3a772aa1-48e5-4654-aeb7-d655d51f7c71", - "metadata": {}, - "outputs": [], - "source": [ - "nc_luts_path = xsarsea.get_test_file('nc_luts_reduce')\n", - "windspeed.register_nc_luts(nc_luts_path)\n", - "nc_luts_path" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d8daa945-0b8c-4c65-8459-0ff2acb1b38d", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed.available_models()" - ] - }, - { - "cell_type": "markdown", - "id": "0cb5a5ef-2af6-4a38-9830-6f2152f9703e", - "metadata": {}, - "source": [ - "To match sarwing outputs, we will flatten `sarwing_ds.owiNesz_cross` with [xsarsea.windspeed.nesz_flattening](../basic_api.rst#xsarsea.windspeed.nesz_flattening) , and compute `dsig_cr`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4d5fe397-05a8-460b-afc3-5ddbebbf80c3", - "metadata": {}, - "outputs": [], - "source": [ - "# nesz cross flattening, and dsig_cr\n", - "nesz_cross_flat = windspeed.nesz_flattening(sarwing_ds.owiNesz_cross, sarwing_ds.owiIncidenceAngle)\n", - "dsig_cr = (1.25 / (sarwing_ds.owiNrcs_cross / nesz_cross_flat )) ** 4.\n", - "\n", - "windspeed_cr = windspeed.invert_from_model(\n", - " sarwing_ds.owiIncidenceAngle,\n", - " sarwing_ds.owiNrcs_cross,\n", - " dsig_cr=dsig_cr,\n", - " model='nc_lut_sarwing_lut_cmodms1ahw',\n", - " **kwargs)\n", - "windspeed_cr = np.abs(windspeed_cr)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "21eba4e7-0d3a-4c6b-a161-e5e92aa59d6d", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed_diff = windspeed_cr - sarwing_ds.owiWindSpeed_Tab_crosspol\n", - "(\n", - " (hv.Image(windspeed_cr).opts(title='xsarsea') + hv.Image(sarwing_ds.owiWindSpeed_Tab_crosspol).opts(title='sarwing' )).opts(hv.opts.Image(cmap='jet', clim=(0,50))) + \n", - " hv.Image(windspeed_diff).opts(cmap='jet', clim=(-0.2,0.2), title='xsarsea-sarwing')\n", - ").opts(hv.opts.Image(colorbar=True, tools=['hover']))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "592300ca-c41d-416e-aca8-e246ec497167", - "metadata": {}, - "outputs": [], - "source": [ - "np.nanmax(np.abs(windspeed_diff)),np.mean(windspeed_diff),np.std(windspeed_diff)" - ] - }, - { - "cell_type": "markdown", - "id": "353cb230-406e-466b-a85b-af69530c155f", - "metadata": {}, - "source": [ - "## Dualpol inversion\n", - "\n", - "Dualpol give better wind inversion results.\n", - "\n", - "Dualpol inversion also returns copol wind." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ab1f5b6f-1be7-4122-9be8-5af0991ed960", - "metadata": {}, - "outputs": [], - "source": [ - "#sarwing_luts_subset_path = xsarsea.get_test_file('sarwing_luts_subset')\n", - "sarwing_luts_subset_path = \"/home/vincelhx/Documents/autoentreprise/IFREMER/sarwing_luts_subset\"\n", - "windspeed.register_sarwing_luts(sarwing_luts_subset_path)\n" - ] - }, - { - "cell_type": "markdown", - "id": "b0854d7f", - "metadata": {}, - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "af0841b0-b5a4-42d8-994d-360cecd32797", - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "wind_co, wind_dual = windspeed.invert_from_model(\n", - " sarwing_ds.owiIncidenceAngle,\n", - " sarwing_ds.owiNrcs,\n", - " sarwing_ds.owiNrcs_cross,\n", - " ancillary_wind=-np.conj(sarwing_ds.owi_ancillary_wind),\n", - " dsig_cr = dsig_cr,\n", - " model=('cmod5n','cmodms1ahw'),\n", - " **kwargs)\n", - "windspeed_co = np.abs(wind_co)\n", - "windspeed_dual = np.abs(wind_dual)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f500e621-ffd9-41c9-bb9b-68f7aba5cf79", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "id": "29379eaa-c340-4f52-9586-dafecb0cfa19", - "metadata": {}, - "source": [ - "## differences with sarwing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0ecead16-b097-44a5-93e4-3a6dc175d0a5", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed_diff = windspeed_dual - sarwing_ds.owiWindSpeed_Tab_dualpol_2steps\n", - "(\n", - " (hv.Image(windspeed_dual).opts(title='xsarsea') + \n", - " hv.Image(sarwing_ds.owiWindSpeed_Tab_dualpol_2steps).opts(title='reference' )).opts(hv.opts.Image(cmap='jet', clim=(0,50))) + \n", - " hv.Image(windspeed_diff).opts(cmap='jet', clim=(-0.5,0.5), title='xsarsea-reference\\nmean=%.4f std=%.4f' % (np.mean(windspeed_diff), np.std(windspeed_diff)))\n", - ").opts(hv.opts.Image(colorbar=True, tools=['hover']))" - ] - }, - { - "cell_type": "markdown", - "id": "d79a3983-5c50-4ba4-a4df-61cb6b9ef955", - "metadata": {}, - "source": [ - "## direction output" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d21b329b-0218-4b81-8994-0e9317fe11a7", - "metadata": {}, - "outputs": [], - "source": [ - "sarwing_ds[\"winddir_dual\"] = (90 - (np.angle(-np.conj(wind_dual),deg=True)) + sarwing_ds.owiHeading)%360\n", - "sub_sarwing_ds = sarwing_ds.isel(line=slice(None, None, 10), sample=slice(None, None, 10))\n", - "\n", - "vectorfield = hv.VectorField(\n", - " (\n", - " sub_sarwing_ds.sample, sub_sarwing_ds.line,\n", - " xsarsea.dir_meteo_to_sample(sarwing_ds[\"winddir_dual\"] ,sub_sarwing_ds.owiHeading),\n", - " np.abs(wind_dual).isel(line=slice(None, None, 10), sample=slice(None, None, 10))\n", - " )\n", - ")\n", - "\n", - "hv.Image(windspeed_dual, kdims=['sample','line']).opts(title='speed and dir', clim=(0,50), cmap='jet') * vectorfield\n" - ] - }, - { - "cell_type": "markdown", - "id": "d1db2cf4", - "metadata": {}, - "source": [ - "## sarwing direction output" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6dc527c2-1f19-447e-9bd3-7e43acaf1a66", - "metadata": {}, - "outputs": [], - "source": [ - "sub_sarwing_ds = sarwing_ds.isel(line=slice(None, None, 10), sample=slice(None, None, 10))\n", - "\n", - "vectorfield = hv.VectorField(\n", - " (\n", - " sub_sarwing_ds.sample, sub_sarwing_ds.line,\n", - " xsarsea.dir_meteo_to_sample(sub_sarwing_ds.owiWindDirection ,sub_sarwing_ds.owiHeading),\n", - " sub_sarwing_ds.owiWindSpeed\n", - " )\n", - ")\n", - "\n", - "hv.Image(sarwing_ds.owiWindSpeed, kdims=['sample','line']).opts(title='speed and dir', clim=(0,50), cmap='jet') * vectorfield\n" - ] - }, - { - "cell_type": "markdown", - "id": "e4016131-62f2-4557-a57f-795741f9cad0", - "metadata": {}, - "source": [ - "## direction difference" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0a48cad7-f888-4d13-8df7-a9191690d7bb", - "metadata": {}, - "outputs": [], - "source": [ - "def normalize_angle(diff):\n", - " return (diff + 180) % 360 - 180\n", - " \n", - "angle_diff = normalize_angle(sarwing_ds[\"winddir_dual\"] - sarwing_ds.owiWindDirection_Tab_dualpol_2steps)\n", - "\n", - "(\n", - " (hv.Image(sarwing_ds[\"winddir_dual\"]).opts(title='xsarsea') + \n", - " hv.Image(sarwing_ds.owiWindDirection_Tab_dualpol_2steps).opts(title='reference' )).opts(hv.opts.Image(cmap='jet', clim=(0,360))) + \n", - " hv.Image(normalize_angle(angle_diff)).opts(cmap='jet', clim=(-5,5), title='xsarsea-reference\\nmean=%.4f std=%.4f' % (np.mean(angle_diff), np.std(angle_diff)))\n", - ").opts(hv.opts.Image(colorbar=True, tools=['hover']))" - ] - } - ], - "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": 5 -} diff --git a/docs/examples/windspeed_retrieval_L1.ipynb b/docs/examples/windspeed_retrieval_L1.ipynb index 2a514b6..f404136 100644 --- a/docs/examples/windspeed_retrieval_L1.ipynb +++ b/docs/examples/windspeed_retrieval_L1.ipynb @@ -10,12 +10,13 @@ }, { "cell_type": "markdown", - "id": "6ddb8b3b-5e67-4dc2-b36b-e234510a80ff", + "id": "f81f9e26-4e02-4f9e-a2dc-7993f5e4be6c", "metadata": {}, "source": [ "> .. warning::\n", " **Use of ancillary wind**\n", - "> On this notebook, we changed from `ancillary_wind = -np.conj(sarwing_ds.owi_ancillary_wind)` to `sarwing_ds.owi_ancillary_wind`` ; then it won't match sarwing results anymore" + "> \n", + "> On this notebook, we changed from `ancillary_wind = -np.conj(sarwing_ds.owi_ancillary_wind)` to `sarwing_ds.owi_ancillary_wind` ; then it won't match sarwing results anymore" ] }, { @@ -35,17 +36,6 @@ "import os,sys,re, cv2" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "69879cf5-1185-47e4-9407-060292b2b44e", - "metadata": {}, - "outputs": [], - "source": [ - "sarwing_luts_subset_path = xsarsea.utils.get_test_file('sarwing_luts_subset')\n", - "windspeed.register_sarwing_luts(sarwing_luts_subset_path)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -61,20 +51,10 @@ }, { "cell_type": "markdown", - "id": "fa29e9e0-ed13-4448-bdbe-c61a6b2f33e3", - "metadata": { - "tags": [] - }, - "source": [ - "## Requirements for inversion" - ] - }, - { - "cell_type": "markdown", - "id": "214f46fe-6169-4a76-a68e-b7c8a205b5fa", + "id": "c3067241-0202-4f11-99f8-7251eb461c78", "metadata": {}, "source": [ - "xsar is required" + "## Requirements for inversion" ] }, { @@ -105,7 +85,7 @@ "outputs": [], "source": [ "safe_path = xsarsea.get_test_file(\"S1A_EW_GRDM_1SDV_20230217T002336_20230217T002412_047268_05AC30_Z005.SAFE\")\n", - "s1meta = xsar.Sentinel1Meta(safe_path)\n" + "s1meta = xsar.Sentinel1Meta(safe_path)" ] }, { @@ -135,15 +115,15 @@ "metadata": {}, "source": [ "land mask:\n", - "not applied yet " + "not applied in the doc" ] }, { "cell_type": "markdown", - "id": "95eb760f-20fc-4fa6-b556-c5804a7e06d4", + "id": "91da3760-5a71-4af7-a673-3e0e9192dcc3", "metadata": {}, "source": [ - "getting associated ancillary data (ecmwf)" + "## getting ancillary data" ] }, { @@ -221,8 +201,9 @@ "tags": [] }, "source": [ - "creation of variables of interest for inversion \n", + "## variables of interest\n", "\n", + "creation of variables of interest for inversion \n", "here we could add a land/ice mask." ] }, @@ -264,20 +245,18 @@ }, { "cell_type": "markdown", - "id": "00eb8706-f7cc-419a-97c9-5898178926bc", - "metadata": { - "tags": [] - }, + "id": "b82cf244-7e9b-481e-bdfb-2599fa8d92ca", + "metadata": {}, "source": [ - "## Inversion" + "## Inversion\n" ] }, { "cell_type": "markdown", - "id": "dc83dd2b-44ff-46ab-b360-a77c52ed2ce3", + "id": "5936708c-c060-429f-89bc-718b812e8ca1", "metadata": {}, "source": [ - "### inversion parameters" + "### inversion parameters\n" ] }, { @@ -288,7 +267,7 @@ "outputs": [], "source": [ "apply_flattening = True\n", - "GMF_VV_NAME = \"cmod5n\"\n", + "GMF_VV_NAME = \"gmf_cmod5n\"\n", "GMF_VH_NAME = \"gmf_s1_v2\"" ] }, @@ -339,17 +318,7 @@ "id": "20801efb-5870-47d9-8fb1-bb3e499f9a34", "metadata": {}, "source": [ - "### get windspeed & direction in dfferent polarizations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "99949fd6-9950-4be8-b9c1-3b906f60b1b3", - "metadata": {}, - "outputs": [], - "source": [ - "windspeed.gmfs.GmfModel.activate_gmfs_impl()" + "### retrieve windspeed & direction in dfferent polarizations" ] }, { @@ -431,6 +400,18 @@ " windspeed_cross=(['line', 'sample'], windspeed_cr))" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d7dff75-3c71-4fb7-a7a9-da823fc55df9", + "metadata": {}, + "outputs": [], + "source": [ + "windspeed_co = dataset.windspeed_co.compute() # compute data if needed\n", + "windspeed_cross = dataset.windspeed_cross\n", + "windspeed_dual = dataset.windspeed_dual.compute()" + ] + }, { "cell_type": "markdown", "id": "7cecc22f-f6be-43ad-bb72-0e243e3b15a9", @@ -446,9 +427,39 @@ "metadata": {}, "outputs": [], "source": [ - "hv.Image(dataset.windspeed_co.compute(), label='wind speed co-pol').opts(cmap='jet',colorbar=True,clim=(0,80),height=400, width=425) + \\\n", - "hv.Image(dataset.windspeed_cross, label='wind speed cr-pol').opts(cmap='jet',colorbar=True,clim=(0,80),height=400, width=425) + \\\n", - "hv.Image(dataset.windspeed_dual.compute(), label='wind speed dual-pol').opts(cmap='jet',colorbar=True,clim=(0,80),height=400, width=425)" + "import matplotlib.pyplot as plt\n", + "import cartopy.crs as ccrs\n", + "import cartopy.feature as cfeature\n", + "\n", + "projection = ccrs.PlateCarree()\n", + "fig, axs = plt.subplots(1, 3, figsize=(18, 6), subplot_kw={'projection': projection})\n", + "\n", + "# Wind Speed Co-Pol\n", + "img0 = axs[0].pcolormesh(dataset.longitude, dataset.latitude, windspeed_co, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())\n", + "axs[0].add_feature(cfeature.COASTLINE)\n", + "axs[0].add_feature(cfeature.BORDERS, linestyle=':')\n", + "axs[0].set_title('Wind Speed Co-Pol')\n", + "gl0 = axs[0].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')\n", + "gl0.top_labels = gl0.right_labels = False\n", + "\n", + "# Wind Speed Cr-Pol\n", + "img1 = axs[1].pcolormesh(dataset.longitude, dataset.latitude, windspeed_cross, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())\n", + "axs[1].add_feature(cfeature.COASTLINE)\n", + "axs[1].add_feature(cfeature.BORDERS, linestyle=':')\n", + "axs[1].set_title('Wind Speed Cr-Pol')\n", + "gl1 = axs[1].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')\n", + "gl1.top_labels = gl1.right_labels = False\n", + "\n", + "# Wind Speed Dual-Pol\n", + "img2 = axs[2].pcolormesh(dataset.longitude, dataset.latitude, windspeed_dual, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())\n", + "axs[2].add_feature(cfeature.COASTLINE)\n", + "axs[2].add_feature(cfeature.BORDERS, linestyle=':')\n", + "axs[2].set_title('Wind Speed Dual-Pol')\n", + "gl2 = axs[2].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')\n", + "gl2.top_labels = gl2.right_labels = False\n", + "\n", + "fig.colorbar(img0, ax=axs, orientation='horizontal', fraction=0.02, pad=0.1, aspect=40)\n", + "plt.show()\n" ] }, { @@ -499,12 +510,20 @@ "hv.Image(dataset.ancillary_wind_speed, kdims=['sample','line']).opts(title='ECMWF speed and dir', clim=(0,50), cmap='jet') * vectorfield\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "df234faf-a32b-459a-b90b-ffc4e3d254a1", + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "id": "93e46998-9c0e-4e55-ac62-b13a7393f907", "metadata": {}, "source": [ - "### save as a level-2 netcdf" + "## saving" ] }, { @@ -573,40 +592,8 @@ "metadata": {}, "outputs": [], "source": [ - "#ds_1000.to_netcdf(\"my_L2_product\")" + "#ds_1000.to_netcdf(\"my_L2_product.nc\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0ef458c7-ddcf-4ab9-98f2-cbe2157377c3", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3ba99aa5-ad0c-48ca-afbb-ec454c4720e4", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d99d681b-fc05-4dde-93d5-2790c2f0bceb", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c6d8da25-2b67-4398-8f54-fd3be7b376b2", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/docs/examples/xsarsea.ipynb b/docs/examples/xsarsea.ipynb index 88b3df8..dd205fb 100644 --- a/docs/examples/xsarsea.ipynb +++ b/docs/examples/xsarsea.ipynb @@ -22,9 +22,7 @@ "outputs": [], "source": [ "import xsar\n", - "import xsarsea\n", - "import os\n", - "import datetime" + "import xsarsea" ] }, { @@ -68,7 +66,6 @@ "source": [ "# get test file. You can replace with an path to other SAFE\n", "filename = xsar.get_test_file('S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_018268_01EB76_Z010.SAFE')\n", - "#filename = xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z010.SAFE')\n", "filename" ] }, @@ -80,18 +77,28 @@ "source": [ "# open the dataset with xsar\n", "sar_ds = xsar.open_dataset(filename, resolution='1000m')\n", - "sar_ds[['sigma0','incidence']]\n" + "sar_ds[['longitude','latitude','sigma0','incidence']]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Sigma0 detrending\n", + "## Sigma0 detrending with xsarsea\n", "\n", "Sigma0 detrending is done by [xsarsea.sigma0_detrend](../basic_api.rst#xsarsea.sigma0_detrend) function\n", "\n", - "As the resulting xarray dataset have the same coordinates as the original sigma0, we can add a `sigma0_detrend` variable to the dataset." + "As the resulting xarray dataset have the same coordinates as the original sigma0, we can add a `sigma0_detrend` variable to the dataset.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default, analytical models are available. We load all luts : \n", + "replace paths by your own path containing all luts\n", + "\n", + "\n" ] }, { @@ -100,7 +107,17 @@ "metadata": {}, "outputs": [], "source": [ - "xsarsea.windspeed.gmfs.GmfModel.activate_gmfs_impl()" + "nc_luts_path = xsarsea.get_test_file('nc_luts_reduce')\n", + "path_cmod7 = xsarsea.get_test_file(\"cmod7_and_python_script\")\n", + "xsarsea.windspeed.register_luts(nc_luts_path, path_cmod7)\n", + "xsarsea.windspeed.available_models()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### default detrend with gmf_cmod5n" ] }, { @@ -119,11 +136,86 @@ "metadata": {}, "outputs": [], "source": [ - "(\n", - " hv.Image(sar_ds.sigma0.sel(pol='VV')).opts(title=\"original sigma0\") \\\n", - " + hv.Image(sar_ds.sigma0_detrend.isel(pol=0)).opts(title=\"detrended sigma0\")\n", - ").opts(hv.opts.Image(cmap='gray', clim=(0,0.4)))" + "original_sigma0_plot = hv.Image(sar_ds.sigma0.sel(pol='VV')).opts(\n", + " title=\"Original Sigma0\",\n", + " cmap='gray',\n", + " clim=(0, 0.4),\n", + " colorbar=True,\n", + " width=400,\n", + " height=400,\n", + " tools=['hover']\n", + ")\n", + "\n", + "detrended_sigma0_plot = hv.Image(sar_ds.sigma0_detrend.isel(pol=0)).opts(\n", + " title=\"Detrended Sigma0\",\n", + " cmap='gray',\n", + " clim=(0, 0.4),\n", + " colorbar=True,\n", + " width=400,\n", + " height=400,\n", + " tools=['hover']\n", + ")\n", + "\n", + "(original_sigma0_plot + detrended_sigma0_plot).cols(2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### custom detrend\n", + "\n", + "very small differences" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "detrend_vv_gmf_cmod7 = xsarsea.sigma0_detrend(sar_ds.sigma0.sel(pol='VV'), sar_ds.incidence, model='gmf_cmod7')\n", + "detrend_vv_gmf_cmod5 = xsarsea.sigma0_detrend(sar_ds.sigma0.sel(pol='VV'), sar_ds.incidence, model='gmf_cmod5')\n", + "detrend_vv_gmf_cmod5n = xsarsea.sigma0_detrend(sar_ds.sigma0.sel(pol='VV'), sar_ds.incidence, model='gmf_cmod5n')\n", + "## HH GMF \n", + "detrend_vv_gmf_cmod5n_zhangA = xsarsea.sigma0_detrend(sar_ds.sigma0.sel(pol='VV'), sar_ds.incidence, model='cmod5n_Rlow_hh_zhangA')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "import numpy as np\n", + "\n", + "fig, axs = plt.subplots(1, 4, figsize=(20, 5)) \n", + "\n", + "datasets = [detrend_vv_gmf_cmod7, detrend_vv_gmf_cmod5, detrend_vv_gmf_cmod5n, detrend_vv_gmf_cmod5n_zhangA]\n", + "titles = ['GMF CMOD7', 'GMF CMOD5', 'GMF CMOD5N', 'GMF CMOD5N ZhangA']\n", + "\n", + "for ax, data, title in zip(axs, datasets, titles):\n", + " img = ax.pcolormesh(sar_ds.longitude, sar_ds.latitude, data, cmap='gray', vmin=0, vmax=0.4) \n", + " ax.set_title(title)\n", + " ax.set_xlabel(\"longitude\")\n", + " ax.set_ylabel(\"latitude\")\n", + " divider = make_axes_locatable(ax)\n", + " cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", + " plt.colorbar(img, cax=cax)\n", + "\n", + "# Improve spacing between plots\n", + "plt.tight_layout()\n", + "plt.show()\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/docs/index.rst b/docs/index.rst index 2eec9a9..ab73989 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -29,12 +29,12 @@ Those examples show how to: * :doc:`examples/streaks` -* :doc:`examples/windspeed_inversion` - * :doc:`examples/windspeed_retrieval_L1` * :doc:`examples/gmfs_and_luts` +* :doc:`examples/create_hh_lut` + Help & Reference ................ @@ -66,8 +66,8 @@ Last documentation build: |today| examples/xsarsea examples/streaks examples/windspeed_retrieval_L1 - examples/windspeed_inversion examples/gmfs_and_luts + examples/create_hh_lut .. toctree:: :maxdepth: 1 diff --git a/src/xsarsea/__init__.py b/src/xsarsea/__init__.py index a78a1b0..0b98bb6 100644 --- a/src/xsarsea/__init__.py +++ b/src/xsarsea/__init__.py @@ -1,8 +1,8 @@ __all__ = ['sigma0_detrend', 'dir_meteo_to_sample', - 'dir_sample_to_meteo', 'dir_meteo_to_oceano', 'dir_oceano_to_meteo', 'dir_to_180', 'dir_to_360', 'get_test_file'] + 'dir_sample_to_meteo', 'dir_meteo_to_oceano', 'dir_oceano_to_meteo', 'dir_to_180', 'dir_to_360', 'get_test_file', 'read_sarwing_owi'] from .utils import get_test_file -from .xsarsea import dir_meteo_to_sample, dir_sample_to_meteo, sigma0_detrend, dir_meteo_to_oceano, dir_oceano_to_meteo, dir_to_180, dir_to_360, read_sarwing_owi +from .xsarsea import dir_meteo_to_sample, dir_sample_to_meteo, sigma0_detrend, dir_meteo_to_oceano, dir_oceano_to_meteo, dir_to_180, dir_to_360, read_sarwing_owi try: from importlib import metadata except ImportError: # for Python<3.8 diff --git a/src/xsarsea/gradients.py b/src/xsarsea/gradients.py index 0db0fc4..9544faa 100644 --- a/src/xsarsea/gradients.py +++ b/src/xsarsea/gradients.py @@ -105,7 +105,7 @@ def histogram(self): exclude_dims=set(self._window_dims.values()), output_core_dims=[['angles'], []], vectorize=True, - output_dtypes=[np.float, np.float] + output_dtypes=[np.float64, np.float64] ) grad_hist = grad_hist.rename( 'weight').assign_coords(angles=angles_bins) @@ -807,7 +807,7 @@ def gradient_histogram(g2, c, angles_bins): k_all = np.round((angle - angles_start) / angles_step) # output array - grads = np.zeros_like(angles_bins, dtype=np.float) + grads = np.zeros_like(angles_bins, dtype=np.float64) # filter nan abs_g2 = np.abs(g2) @@ -862,7 +862,7 @@ def circ_smooth(hist): input_core_dims=[["angles"], ["kernel_len"]], output_core_dims=[['angles']], vectorize=True, - output_dtypes=[np.float], + output_dtypes=[np.float64], ) # unwrap diff --git a/src/xsarsea/utils.py b/src/xsarsea/utils.py index fa7de65..393c372 100644 --- a/src/xsarsea/utils.py +++ b/src/xsarsea/utils.py @@ -78,27 +78,30 @@ def get_test_file(fname, iszip=True): base_url = 'https://cyclobs.ifremer.fr/static/sarwing_datarmor/xsardata' file_url = '%s/%s.zip' % (base_url, fname) - if not iszip: - import urllib - file_url = '%s/%s' % (base_url, fname) - warnings.warn("Downloading %s" % file_url) - urllib.request.urlretrieve( - file_url, os.path.join(config['data_dir'], fname)) + if not os.path.exists(os.path.join(res_path, fname)): - else: - if not os.path.exists(os.path.join(res_path, fname)): + if not iszip: + import urllib + file_url = '%s/%s' % (base_url, fname) warnings.warn("Downloading %s" % file_url) - with fsspec.open( - 'filecache::%s' % file_url, - https={'client_kwargs': { - 'timeout': aiohttp.ClientTimeout(total=3600)}}, - filecache={'cache_storage': os.path.join( - os.path.join(config['data_dir'], 'fsspec_cache'))} - ) as f: - - warnings.warn("Unzipping %s" % os.path.join(res_path, fname)) - with zipfile.ZipFile(f, 'r') as zip_ref: - zip_ref.extractall(res_path) + urllib.request.urlretrieve( + file_url, os.path.join(res_path, fname)) + + else: + if not os.path.exists(os.path.join(res_path, fname)): + warnings.warn("Downloading %s" % file_url) + with fsspec.open( + 'filecache::%s' % file_url, + https={'client_kwargs': { + 'timeout': aiohttp.ClientTimeout(total=3600)}}, + filecache={'cache_storage': os.path.join( + os.path.join(config['data_dir'], 'fsspec_cache'))} + ) as f: + + warnings.warn("Unzipping %s" % + os.path.join(res_path, fname)) + with zipfile.ZipFile(f, 'r') as zip_ref: + zip_ref.extractall(res_path) return os.path.join(res_path, fname) diff --git a/src/xsarsea/windspeed/__init__.py b/src/xsarsea/windspeed/__init__.py index 98f808d..2e71cd8 100644 --- a/src/xsarsea/windspeed/__init__.py +++ b/src/xsarsea/windspeed/__init__.py @@ -2,10 +2,11 @@ windspeed module, for retrieving wind speed from sigma0 and models. """ __all__ = ['invert_from_model', 'available_models', 'get_model', 'register_cmod7', - 'register_sarwing_luts', 'register_nc_luts', 'register_luts', 'nesz_flattening', 'GmfModel'] + 'register_pickle_luts', 'register_nc_luts', 'register_luts', 'nesz_flattening', 'GmfModel','Model'] + from .windspeed import invert_from_model from .models import available_models, get_model, register_nc_luts, register_luts -from .sarwing_luts import register_sarwing_luts +from .pickle_luts import register_pickle_luts from .cmod7 import register_cmod7 from .utils import nesz_flattening, get_dsig from .gmfs import GmfModel diff --git a/src/xsarsea/windspeed/gmfs.py b/src/xsarsea/windspeed/gmfs.py index 54fe933..f2e6203 100644 --- a/src/xsarsea/windspeed/gmfs.py +++ b/src/xsarsea/windspeed/gmfs.py @@ -81,7 +81,7 @@ def inner(func): cls._name_prefix, gmf_name)) wspd_range = kwargs.pop('wspd_range', None) - + if wspd_range is None: if len(set(pol)) == 1: # copol @@ -261,7 +261,6 @@ def __call__(self, inc, wspd, phi=None, broadcast=False, numba=True): for v in [inc, wspd, phi] if v is not None) # logger.debug("Initial input shapes, inc: %s, wspd: %s, phi: %s", # inc.shape, wspd.shape, phi.shape if phi is not None else "None") - # if all 1d, will return 2d or 3d shape('incidence', 'wspd', 'phi'), unless broadcast is True all_1d = all(hasattr(v, 'ndim') and v.ndim == 1 for v in [ diff --git a/src/xsarsea/windspeed/gmfs_impl.py b/src/xsarsea/windspeed/gmfs_impl.py index 96c8b76..1669b72 100644 --- a/src/xsarsea/windspeed/gmfs_impl.py +++ b/src/xsarsea/windspeed/gmfs_impl.py @@ -23,7 +23,7 @@ def gmf_cmod5_generic(neutral=False): 6.2437, 2.3893, 0.3249, 4.159, 1.693]) name = 'gmf_cmod5n' - @GmfModel.register(name, wspd_range=[0.2, 50.], pol='VV', units='linear', defer=True) + @GmfModel.register(name, wspd_range=[0.2, 50.], pol='VV', units='linear', defer=False) def gmf_cmod5(inc, wspd, phi): zpow = 1.6 thetm = 40. @@ -81,7 +81,7 @@ def gmf_cmod5(inc, wspd, phi): gmf_cmod5_generic(neutral=True) -@GmfModel.register(wspd_range=[0.2, 50.], pol='VV', units='linear', defer=True) +@GmfModel.register(wspd_range=[0.2, 50.], pol='VV', units='linear', defer=False) def gmf_cmodifr2(inc_angle, wind_speed, wind_dir): C = np.zeros(26) @@ -194,7 +194,7 @@ def gmf_cmodifr2(inc_angle, wind_speed, wind_dir): # return sig -@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=True) +@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=False) def gmf_rs2_v2(incidence, speed, phi=None): """ Radarsat-2 VH GMF : relation between sigma0, incidence and windspeed. @@ -249,7 +249,7 @@ def gmf_rs2_v2(incidence, speed, phi=None): return sig_Final -@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=True) +@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=False) def gmf_s1_v2(incidence, speed, phi=None): """ Sentinel-1 VH GMF : relation between sigma0, incidence and windspeed. @@ -304,7 +304,7 @@ def gmf_s1_v2(incidence, speed, phi=None): return sig_Final -@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=True) +@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear', defer=False) def gmf_rcm_noaa(incidence, speed, phi=None): """ Radarsat Consteallation Mission VH GMF : relation between sigma0, incidence and windspeed. diff --git a/src/xsarsea/windspeed/models.py b/src/xsarsea/windspeed/models.py index be41ea5..48197c9 100644 --- a/src/xsarsea/windspeed/models.py +++ b/src/xsarsea/windspeed/models.py @@ -417,13 +417,10 @@ def register_nc_luts(topdir, gmf_names=None): >>> xsarsea.windspeed.register_nc_luts(xsarsea.get_test_file('nc_luts_subset')) - register all sarwing lut from ifremer path + register all pickle lut from ifremer path - >>> xsarsea.windspeed.register_sarwing_luts('/home/datawork-cersat-public/cache/project/sarwing/xsardata/nc_luts') + >>> xsarsea.windspeed.register_pickle_luts('/home/datawork-cersat-public/cache/project/sarwing/xsardata/nc_luts') - Notes - _____ - Sarwing lut can be downloaded from https://cyclobs.ifremer.fr/static/sarwing_datarmor/xsardata/nc_luts See Also -------- @@ -528,14 +525,14 @@ def register_luts(topdir=None, topdir_cmod7=None): ---------- topdir: str top dir path to nc luts. - + topdir_cmod7: str top dir path to cmod7 luts. kwargs: dict kwargs to pass to register_nc_luts """ - + # register gmf luts import xsarsea.windspeed as windspeed windspeed.GmfModel.activate_gmfs_impl() @@ -547,3 +544,4 @@ def register_luts(topdir=None, topdir_cmod7=None): # register cmod7 if topdir_cmod7 != None: windspeed.register_cmod7(topdir_cmod7) + \ No newline at end of file diff --git a/src/xsarsea/windspeed/sarwing_luts.py b/src/xsarsea/windspeed/pickle_luts.py similarity index 88% rename from src/xsarsea/windspeed/sarwing_luts.py rename to src/xsarsea/windspeed/pickle_luts.py index d1692c8..1e17a04 100644 --- a/src/xsarsea/windspeed/sarwing_luts.py +++ b/src/xsarsea/windspeed/pickle_luts.py @@ -6,7 +6,7 @@ from .models import LutModel -class SarwingLutModel(LutModel): +class PickleLutModel(LutModel): _name_prefix = 'sarwing_lut__' _priority = 10 @@ -19,7 +19,7 @@ def _raw_lut(self, **kwargs): if not os.path.isdir(self.path): raise FileNotFoundError(self.path) - logger.info('load sarwing lut from %s' % self.path) + logger.info('load pickle lut from %s' % self.path) sigma0_db_path = os.path.join(self.path, 'sigma.npy') sigma0_db = np.ascontiguousarray(np.transpose(np.load(sigma0_db_path))) @@ -70,7 +70,7 @@ def _raw_lut(self, **kwargs): return da_sigma0_db.transpose(*final_dims) -def register_sarwing_luts(path): +def register_pickle_luts(path): """ Register LUTs from a specified path. The path can be a directory containing multiple LUTs or a single LUT file. @@ -85,11 +85,11 @@ def register_sarwing_luts(path): -------- Register a single LUT: - >>> xsarsea.windspeed.register_sarwing_luts(xsarsea.get_test_file('sarwing_luts_subset/GMF_cmodms1ahw')) + >>> xsarsea.windspeed.register_pickle_luts(xsarsea.get_test_file('sarwing_luts_subset/GMF_cmodms1ahw')) Register all LUTs from a directory: - >>> xsarsea.windspeed.register_sarwing_luts('/home/datawork-cersat-public/cache/project/sarwing/GMFS/v1.6') + >>> xsarsea.windspeed.register_pickle_luts('/home/datawork-cersat-public/cache/project/sarwing/GMFS/v1.6') Notes ----- @@ -102,8 +102,8 @@ def register_sarwing_luts(path): """ def register_lut(file_path): - sarwing_name = os.path.basename(file_path) - name = sarwing_name.replace('GMF_', SarwingLutModel._name_prefix) + name = os.path.basename(file_path) + name = name.replace('GMF_', PickleLutModel._name_prefix) # Guess available pols from filenames if os.path.exists(os.path.join(file_path, 'wind_speed_and_direction.pkl')): pol = 'VV' @@ -112,7 +112,7 @@ def register_lut(file_path): else: pol = None - sarwing_model = SarwingLutModel(name, file_path, pol=pol) + pickleLutmodel = PickleLutModel(name, file_path, pol=pol) last_folder_name = os.path.basename(os.path.normpath(path)) if last_folder_name.startswith('GMF_'): diff --git a/src/xsarsea/windspeed/utils.py b/src/xsarsea/windspeed/utils.py index bdd1e76..33cc72b 100644 --- a/src/xsarsea/windspeed/utils.py +++ b/src/xsarsea/windspeed/utils.py @@ -48,9 +48,12 @@ def sigmoid(x, c0, c1, d0, d1): elif name == "sarwing_lut_cmodms1ahw": return (1.25 / (sigma0_cr / nesz_cr)) ** 4. + elif name == "nc_lut_cmodms1ahw": + return (1.25 / (sigma0_cr / nesz_cr)) ** 4. + else: raise ValueError( - "dsig names different than 'gmf_s1_v2' or 'gmf_rs2_v2' or 'sarwing_lut_cmodms1ahw' are not handled. You can compute your own dsig_cr.") + "dsig names different than 'gmf_s1_v2' or 'gmf_rs2_v2' or 'sarwing_lut_cmodms1ahw' or 'nc_lut_cmodms1ahw' are not handled. You can compute your own dsig_cr.") def nesz_flattening(noise, inc): diff --git a/test/test_xsarsea.py b/test/test_xsarsea.py index dca1f30..dfd42ac 100644 --- a/test/test_xsarsea.py +++ b/test/test_xsarsea.py @@ -5,7 +5,7 @@ import dask.array as da -@windspeed.gmfs.GmfModel.register(inc_range=[17., 50.], wspd_range=[3., 80.], pol='VH', units='linear', defer=True) +@windspeed.gmfs.GmfModel.register(inc_range=[17., 50.], wspd_range=[3., 80.], pol='VH', units='linear', defer=False) def gmf_dummy(inc, wspd, phi=None): a0 = 0.00013106836021008122 a1 = -4.530598283705591e-06 @@ -91,7 +91,7 @@ def test_inversion(): ]) sarwing_luts_subset_path = xsarsea.get_test_file('sarwing_luts_subset') - windspeed.sarwing_luts.register_sarwing_luts(sarwing_luts_subset_path) + windspeed.pickle_luts.register_sarwing_luts(sarwing_luts_subset_path) nesz_cross_flat = windspeed.nesz_flattening( sarwing_ds.owiNesz_cross, sarwing_ds.owiIncidenceAngle)