diff --git a/docs/examples/streaks.ipynb b/docs/examples/streaks.ipynb index 9742703..41e118f 100644 --- a/docs/examples/streaks.ipynb +++ b/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()" ] }, { diff --git a/src/xsarsea/utils.py b/src/xsarsea/utils.py index bd5f191..fa7de65 100644 --- a/src/xsarsea/utils.py +++ b/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 diff --git a/src/xsarsea/windspeed/cmod7.py b/src/xsarsea/windspeed/cmod7.py index 1d20398..52006ec 100644 --- a/src/xsarsea/windspeed/cmod7.py +++ b/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): diff --git a/src/xsarsea/windspeed/gmfs.py b/src/xsarsea/windspeed/gmfs.py index 9e688ae..54fe933 100644 --- a/src/xsarsea/windspeed/gmfs.py +++ b/src/xsarsea/windspeed/gmfs.py @@ -1,5 +1,4 @@ import numpy as np -import warnings from ..utils import timing from .utils import logger from functools import lru_cache @@ -7,7 +6,6 @@ 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 [ diff --git a/src/xsarsea/windspeed/gmfs_impl.py b/src/xsarsea/windspeed/gmfs_impl.py index 5b2379b..a650186 100644 --- a/src/xsarsea/windspeed/gmfs_impl.py +++ b/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 diff --git a/src/xsarsea/windspeed/models.py b/src/xsarsea/windspeed/models.py index 5ce6c8c..20dab2d 100644 --- a/src/xsarsea/windspeed/models.py +++ b/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) diff --git a/src/xsarsea/windspeed/utils.py b/src/xsarsea/windspeed/utils.py index 44ee28a..bdd1e76 100644 --- a/src/xsarsea/windspeed/utils.py +++ b/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 diff --git a/src/xsarsea/xsarsea.py b/src/xsarsea/xsarsea.py index 652b7ae..f94e401 100644 --- a/src/xsarsea/xsarsea.py +++ b/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)