dsig + sigma0_detrend review

docs/examples/streaks.ipynb

@@ -51,11 +51,11 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "f370d8c1-c1a8-46a9-bb75-c9386068156a",
+   "id": "29afbe30",
    "metadata": {},
    "outputs": [],
    "source": [
-    "xsarsea.windspeed.gmfs.GmfModel.activate_gmfs_impl()"
+    "xsarsea.windspeed.available_models()"
    ]
   },
   {

src/xsarsea/utils.py

@@ -3,7 +3,6 @@
 import fsspec
 import aiohttp
 import zipfile
-from pathlib import Path
 import yaml
 import time
 import logging
@@ -40,16 +39,16 @@ def _load_config():
     -------
     dict
     """
-    user_config_file = Path('~/.xsarsea/config.yml').expanduser()
+    user_config_file = os.path.expanduser('~/.xsarsea/config.yml')
     default_config_file = files('xsarsea').joinpath('config.yml')
 
-    if user_config_file.exists():
+    if os.path.exists(user_config_file):
         config_file = user_config_file
     else:
         config_file = default_config_file
 
     config = yaml.load(
-        config_file.open(),
+        open(config_file, 'r'),
         Loader=yaml.FullLoader)
     return config
 

src/xsarsea/windspeed/cmod7.py

@@ -1,8 +1,8 @@
 import os
 import xarray as xr
-from .utils import logger
+from xsarsea.windspeed.utils import logger
 import numpy as np
-from .models import LutModel
+from xsarsea.windspeed.models import LutModel
 
 
 class Cmod7Model(LutModel):

src/xsarsea/windspeed/gmfs.py

@@ -1,13 +1,11 @@
 import numpy as np
-import warnings
 from ..utils import timing
 from .utils import logger
 from functools import lru_cache
 from numba import njit, vectorize, guvectorize, float64, float32
 import xarray as xr
 import dask.array as da
 from .models import Model
-import time
 
 
 class GmfModel(Model):
@@ -83,6 +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
@@ -260,8 +259,9 @@ def __call__(self, inc, wspd, phi=None, broadcast=False, numba=True):
         # if all scalar, will return scalar
         all_scalar = all(np.isscalar(v)
                          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")
+        # 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 [

src/xsarsea/windspeed/gmfs_impl.py

@@ -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.], inc_range=[17., 50.], pol='VH', units='linear', defer=True)
 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.], inc_range=[17., 50.], pol='VH', units='linear', defer=True)
 def gmf_s1_v2(incidence, speed, phi=None):
     """
     Sentinel-1 VH GMF : relation between sigma0, incidence and windspeed. 
@@ -332,7 +332,6 @@ def gmf_rcm_noaa(incidence, speed, phi=None):
                        0.12713502524515713, 4.2806865431046752])
 
     # Z1
-
     a0_Z1 = Z1_p[0]
     b0_Z1 = Z1_p[1]
     b1_Z1 = Z1_p[2]
@@ -359,38 +358,3 @@ def gmf_rcm_noaa(incidence, speed, phi=None):
     sigmoid_Z2 = 1 / (1 + np.exp(-c2*(speed-c3)))
     sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
     return sig_Final
-
-
-def get_PR_Mouche1(inc_angle, wind_dir, wind_speed=None):
-    """
-    Get the polarization ratio VV over HH using the Mouche model.
-    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
-    """
-    # getRatioVVoverHH
-    theta = inc_angle
-    phi = wind_dir
-
-    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
-
-    pr = C0_theta + C1_theta * \
-        np.cos(np.deg2rad(phi)) + C2_theta*np.cos(2*np.deg2rad(phi))
-
-    return pr

src/xsarsea/windspeed/models.py

@@ -35,8 +35,8 @@ def __init__(self, name, **kwargs):
         self.resolution = kwargs.pop('resolution', None)
 
         if not hasattr(self, 'inc_range'):
-            self.inc_range = [17., 50.]
-
+            # self.inc_range = [17., 50.]
+            self.inc_range = [16., 66.]
         # steps for generated luts
         self.inc_step_lr = kwargs.pop(
             'inc_step_lr', 1.0)

src/xsarsea/windspeed/utils.py

@@ -1,10 +1,8 @@
-import os
-from pathlib import Path
 import warnings
 import numpy as np
 import yaml
 from importlib_resources import files
-
+import os
 import logging
 logger = logging.getLogger('xsarsea.windspeed')
 # logger.addHandler(logging.NullHandler())
@@ -37,7 +35,7 @@ def get_dsig(name, inc, sigma0_cr, nesz_cr):
         def sigmoid(x, c0, c1, d0, d1):
             sig = d0 + d1 / (1 + np.exp(-c0*(x-c1)))
             return sig
-        
+
         poptsig = np.array([1.57952257, 25.61843791,  1.46852088,  1.4058646])
         c = sigmoid(inc, *poptsig)
         return (1 / np.sqrt(b*(sigma0_cr / nesz_cr)**c))
@@ -48,11 +46,11 @@ def sigmoid(x, c0, c1, d0, d1):
         return (1 / np.sqrt(b*(sigma0_cr / nesz_cr)**c))
 
     elif name == "sarwing_lut_cmodms1ahw":
-        return 1.25 / (sigma0_cr / nesz_cr) ** 4.
+        return (1.25 / (sigma0_cr / nesz_cr)) ** 4.
 
     else:
         raise ValueError(
-            "dsig names different than 'gmf_s1_v2' or 'gmf_rs2_v2' or 'gmf_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' are not handled. You can compute your own dsig_cr.")
 
 
 def nesz_flattening(noise, inc):
@@ -135,18 +133,18 @@ def _load_config_luts(config_path):
     dict
     """
 
-    user_config_file = Path(config_path)
+    user_config_file = open(config_path, 'r')
     default_config_file = files('xsarsea').joinpath("windspeed").joinpath(
         'config_luts_default_direct_01_01_10.yml')
 
-    if user_config_file.exists():
+    if os.exists(user_config_file.exists):
         config_file = user_config_file
     else:
         # logger.info(f"Using default config file {default_config_file}")
         # config_file = default_config_file
         raise FileNotFoundError(f"Config file {user_config_file} not found")
     config = yaml.load(
-        config_file.open(),
+        open(config_file),
         Loader=yaml.FullLoader)
 
     return config

src/xsarsea/xsarsea.py

@@ -26,7 +26,15 @@ def sigma0_detrend(sigma0, inc_angle, wind_speed_gmf=np.array([10.]), wind_dir_g
     """
     model = get_model(model)
 
+    if wind_speed_gmf.ndim > 1 or wind_dir_gmf.ndim > 1:
+        raise ValueError("wind_speed_gmf and wind_dir_gmf must be 0D or 1D")
+    for var in [wind_speed_gmf, wind_dir_gmf]:
+        if var.ndim == 1 and var.size > 1:
+            raise ValueError(
+                "wind_speed_gmf and wind_dir_gmf size must be 1 or 0")
+
     # get model for one line (all incidences)
+    """
     try:
         # if using dask, model is unpicklable. The workaround is to use map_blocks
         sigma0_gmf_sample = inc_angle.isel(line=0).map_blocks(
@@ -37,8 +45,15 @@ def sigma0_detrend(sigma0, inc_angle, wind_speed_gmf=np.array([10.]), wind_dir_g
     except AttributeError:
         # this should be the standard way
         # see https://github.com/dask/distributed/issues/3450#issuecomment-585255484
-        sigma0_gmf_sample = model(inc_angle.isel(
-            line=0), wind_speed_gmf, wind_dir_gmf, broadcast=True)
+    """
+    sigma0_gmf_sample = model(inc_angle.isel(
+        line=0), wind_speed_gmf, wind_dir_gmf, broadcast=True)
+
+    for var in ["wspd", "phi", "incidence"]:
+        if var in sigma0_gmf_sample.dims:
+            sigma0_gmf_sample = sigma0_gmf_sample.squeeze(var)
+        if var in sigma0_gmf_sample.coords:
+            sigma0_gmf_sample = sigma0_gmf_sample.drop_vars(var)
 
     gmf_ratio_sample = sigma0_gmf_sample / np.nanmean(sigma0_gmf_sample)
     detrended = sigma0 / gmf_ratio_sample.broadcast_like(sigma0)