diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f74a9ae..0713437 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,69 +15,69 @@ jobs: steps: # Checkout the code - - name: Checkout code - uses: actions/checkout@v4 + - name: Checkout code + uses: actions/checkout@v4 - - uses: mamba-org/setup-micromamba@v1 - with: - micromamba-version: "1.5.9-1" # any version from https://github.com/mamba-org/micromamba-releases - channels: tcevaer, conda-forge - init-shell: bash - post-cleanup: "all" + - uses: mamba-org/setup-micromamba@v1 + with: + micromamba-version: "1.5.9-1" # any version from https://github.com/mamba-org/micromamba-releases + channels: tcevaer, conda-forge + init-shell: bash + post-cleanup: "all" - - name: Configure Conda channel priority to disabled - run: | - conda config --set channel_priority disabled + - name: Configure Conda channel priority to disabled + run: | + conda config --set channel_priority disabled - - name: Create environment and install tools - run: micromamba create -n grdwind_env pytest conda-build boa python=3.10 -y + - name: Create environment and install tools + run: micromamba create -n grdwind_env pytest conda-build boa python=3.10 -y - - name: Build package - run: | - cd recipe - eval "$(micromamba shell hook --shell bash)" - micromamba activate grdwind_env - conda mambabuild . + - name: Build package + run: | + cd recipe + eval "$(micromamba shell hook --shell bash)" + micromamba activate grdwind_env + conda mambabuild . - # Install the built package into the environment - - name: Install the built package - run: | - eval "$(micromamba shell hook --shell bash)" - micromamba activate grdwind_env - conda install --use-local grdwindinversion -y + # Install the built package into the environment + - name: Install the built package + run: | + eval "$(micromamba shell hook --shell bash)" + micromamba activate grdwind_env + conda install --use-local grdwindinversion -y - # Cache the test data if previously downloaded (up to 10 GB limit for the cache) - # WARNING : modify the key if the data is modified !! - - name: Cache test data - uses: actions/cache@v4 - id: cache - with: - path: ./test_data - key: test-data-v3 - restore-keys: test-data-v3 + # Cache the test data if previously downloaded (up to 10 GB limit for the cache) + # WARNING : modify the key if the data is modified !! + - name: Cache test data + uses: actions/cache@v4 + id: cache + with: + path: ./test_data + key: test-data-v3 + restore-keys: test-data-v3 - # Download test data if not already cached - - name: Download test data - if: steps.cache.outputs.cache-hit != 'true' # Only download if cache miss - run: | + # Download test data if not already cached + - name: Download test data + if: steps.cache.outputs.cache-hit != 'true' # Only download if cache miss + run: | mkdir -p ./test_data/ wget https://cloud.ifremer.fr/index.php/s/ExLQ2TnYAqozPWE/download -O /tmp/ecmwf.zip unzip /tmp/ecmwf.zip -d ./test_data/ wget https://cloud.ifremer.fr/index.php/s/kRgdOOPsjoZieZR/download -O /tmp/l1.zip unzip /tmp/l1.zip -d ./test_data/ - timeout-minutes: 200 # Adjust depending on the size of your data + timeout-minutes: 200 # Adjust depending on the size of your data - # Set up xsar configuration - - name: Setup xsar configuration - run: | + # Set up xsar configuration + - name: Setup xsar configuration + run: | mkdir -p ~/.xsar echo "data_dir: /tmp" > ~/.xsar/config.yaml echo "auxiliary_dir: ./test_data/auxiliary" >> ~/.xsar/config.yaml echo "path_dataframe_aux: ./test_data/auxiliary/active_aux.csv" >> ~/.xsar/config.yaml - # Set up grdwindinversion configuration - - name: Setup grdwindinversion configuration - run: | + # Set up grdwindinversion configuration + - name: Setup grdwindinversion configuration + run: | mkdir -p ~/.grdwindinversion echo "'ecmwf_0100_1h': ./test_data/ECMWF/forecast/hourly/0100deg/netcdf_light/%Y/%j/ECMWF_FORECAST_0100_%Y%m%d%H%M_10U_10V.nc" > ~/.grdwindinversion/data_config.yaml echo "'ecmwf_0125_1h': ./test_data/ECMWF/0.125deg/1h/forecasts/%Y/%j/ecmwf_%Y%m%d%H%M.nc" >> ~/.grdwindinversion/data_config.yaml @@ -85,9 +85,9 @@ jobs: #echo "'lut_cmod7_path': './test_data/GMFS/v1.6/GMF_cmod7_official/cmod7_and_python_script'" >> ~/.grdwindinversion/data_config.yaml #echo "'lut_ms1ahw_path': './test_data/GMFS/v1.6/GMF_cmodms1ahw'" >> ~/.grdwindinversion/data_config.yaml - # Run the tests - - name: Run tests - run: | + # Run the tests + - name: Run tests + run: | eval "$(micromamba shell hook --shell bash)" micromamba activate grdwind_env pytest diff --git a/grdwindinversion/config_prod.yaml b/grdwindinversion/config_prod.yaml index 4e34246..5550fa2 100644 --- a/grdwindinversion/config_prod.yaml +++ b/grdwindinversion/config_prod.yaml @@ -1,5 +1,7 @@ no_subdir: True winddir_convention: "meteorological" +add_gradientsfeatures: False +add_nrcs_model: False S1A: GMF_VV_NAME: "gmf_cmod5n" GMF_VH_NAME: "gmf_s1_v2" diff --git a/grdwindinversion/config_prod_recal.yaml b/grdwindinversion/config_prod_recal.yaml index 98fedfb..8100d83 100644 --- a/grdwindinversion/config_prod_recal.yaml +++ b/grdwindinversion/config_prod_recal.yaml @@ -6,9 +6,9 @@ S1A: apply_flattening: True recalibration: True ancillary: "ecmwf" - inc_step: 0.1 - wspd_step: 0.1 - phi_step: 1.0 + inc_step: 0.3 + wspd_step: 0.3 + phi_step: 2.0 resolution: "high" S1B: GMF_VV_NAME: "gmf_cmod5n" diff --git a/grdwindinversion/config_prod_recal_streaks_nrcsmod.yaml b/grdwindinversion/config_prod_recal_streaks_nrcsmod.yaml new file mode 100644 index 0000000..9f5cdb8 --- /dev/null +++ b/grdwindinversion/config_prod_recal_streaks_nrcsmod.yaml @@ -0,0 +1,48 @@ +no_subdir: True +winddir_convention: "meteorological" +add_gradientsfeatures: True +add_nrcs_model: True +S1A: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: True + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +S1B: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: True + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RS2: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rs2_v2" + dsig_VH_NAME: "gmf_rs2_v2" + apply_flattening: False + recalibration: True + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RCM: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rcm_noaa" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: True + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" diff --git a/grdwindinversion/config_prod_streaks.yaml b/grdwindinversion/config_prod_streaks.yaml new file mode 100644 index 0000000..93715be --- /dev/null +++ b/grdwindinversion/config_prod_streaks.yaml @@ -0,0 +1,48 @@ +no_subdir: True +winddir_convention: "meteorological" +add_gradientsfeatures: True +add_nrcs_model: False +S1A: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +S1B: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RS2: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rs2_v2" + dsig_VH_NAME: "gmf_rs2_v2" + apply_flattening: False + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RCM: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rcm_noaa" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" diff --git a/grdwindinversion/config_prod_streaks_nrcsmod.yaml b/grdwindinversion/config_prod_streaks_nrcsmod.yaml new file mode 100644 index 0000000..dd06fea --- /dev/null +++ b/grdwindinversion/config_prod_streaks_nrcsmod.yaml @@ -0,0 +1,48 @@ +no_subdir: True +winddir_convention: "meteorological" +add_gradientsfeatures: True +add_nrcs_model: True +S1A: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +S1B: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_s1_v2" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RS2: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rs2_v2" + dsig_VH_NAME: "gmf_rs2_v2" + apply_flattening: False + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" +RCM: + GMF_VV_NAME: "gmf_cmod5n" + GMF_VH_NAME: "gmf_rcm_noaa" + dsig_VH_NAME: "gmf_s1_v2" + apply_flattening: True + recalibration: False + ancillary: "ecmwf" + inc_step: 0.1 + wspd_step: 0.1 + phi_step: 1.0 + resolution: "high" diff --git a/grdwindinversion/inversion.py b/grdwindinversion/inversion.py index 3e294e6..2006a43 100644 --- a/grdwindinversion/inversion.py +++ b/grdwindinversion/inversion.py @@ -1,3 +1,4 @@ +import tempfile import traceback import xsar @@ -15,14 +16,12 @@ import re import string import os -from grdwindinversion.streaks import get_streaks -from grdwindinversion.utils import check_incidence_range, get_pol_ratio_name +from grdwindinversion.utils import check_incidence_range, get_pol_ratio_name, timing from grdwindinversion.load_config import getConf # optional debug messages import logging -logging.basicConfig() -logging.getLogger('xsarsea.windspeed').setLevel( - logging.INFO) # or .setLevel(logging.INFO) +logger = logging.getLogger('grdwindinversion.inversion') +logger.addHandler(logging.NullHandler()) def getSensorMetaDataset(filename): @@ -96,12 +95,10 @@ def getOutputName2(input_file, outdir, sensor, meta, subdir=True): new_format = f"{MISSIONID.lower()}--owi-xx-{meta_start_date.lower()}-{meta_stop_date.lower()}-_____-_____.nc" elif sensor == 'RCM': regex = re.compile( - "([A-Z0-9]+)_OK([0-9]+)_PK([0-9]+)_(.*?)_(.*?)_(.*?)_(.*?)_(.*?)_(.*?)_(.*?)") - template = string.Template( - "${MISSIONID}_OK${DATA1}_PK${DATA2}_${DATA3}_${BEAM}_${DATE}_${TIME}_${POLARIZATION1}_${POLARIZATION2}_${PRODUCT}") + r"(RCM[0-9])_OK([0-9]+)_PK([0-9]+)_([0-9]+)_([A-Z]+)_(\d{8})_(\d{6})_([A-Z]{2}(?:_[A-Z]{2})?)_([A-Z]+)$") match = regex.match(basename_match) - MISSIONID, DATA1, DATA2, DATA3, BEAM_MODE, DATE, TIME, POLARIZATION1, POLARIZATION2, LAST = match.groups() - new_format = f"{MISSIONID.lower()}-{BEAM_MODE.lower()}-owi-xx-{meta_start_date.lower()}-{meta_stop_date.lower()}-_____-_____.nc" + MISSIONID, DATA1, DATA2, DATA3, BEAM, DATE, TIME, POLARIZATION, PRODUCT = match.groups() + new_format = f"{MISSIONID.lower()}-{BEAM.lower()}-owi-xx-{meta_start_date.lower()}-{meta_stop_date.lower()}-_____-_____.nc" else: raise ValueError( "sensor must be S1A|S1B|RS2|RCM, got sensor %s" % sensor) @@ -209,6 +206,7 @@ def getAncillary(meta, ancillary_name='ecmwf'): ancillary_name) +@timing(logger=logger.debug) def inverse(dual_pol, inc, sigma0, sigma0_dual, ancillary_wind, dsig_cr, model_co, model_cross, **kwargs): """ Invert sigma0 to retrieve wind using model (lut or gmf). @@ -282,7 +280,8 @@ def inverse(dual_pol, inc, sigma0, sigma0_dual, ancillary_wind, dsig_cr, model_c return wind_co, None, None -def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flattening): +@timing(logger=logger.debug) +def makeL2asOwi(xr_dataset, config): """ Rename xr_dataset variables and attributes to match naming convention. @@ -290,12 +289,8 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte ---------- xr_dataset: xarray.Dataset dataset to rename - dual_pol: bool - True if dualpol, False if singlepol - copol: str - copolarization name - crosspol: str - crosspolarization name + config: dict + configuration dict Returns ------- @@ -320,13 +315,27 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte 'winddir_co': 'owiWindDirection_co', 'ancillary_wind_speed': 'owiAncillaryWindSpeed', 'ancillary_wind_direction': 'owiAncillaryWindDirection', - 'sigma0_detrend': 'owiNrcs_detrend' + 'sigma0_detrend': 'owiNrcs_detrend', }) if "offboresight" in xr_dataset: xr_dataset = xr_dataset.rename( {"offboresight": "owiOffboresightAngle"}) + if config["add_nrcs_model"]: + xr_dataset = xr_dataset.rename( + {"ancillary_nrcs": "owiAncillaryNrcs"}) + xr_dataset.owiAncillaryNrcs.attrs["units"] = "m^2 / m^2" + xr_dataset.owiAncillaryNrcs.attrs[ + "long_name"] = f"Ancillary Normalized Radar Cross Section - simulated from {config['l2_params']['copol_gmf']} & ancillary wind" + + if config["l2_params"]["dual_pol"]: + xr_dataset = xr_dataset.rename( + {"ancillary_nrcs_cross": "owiAncillaryNrcs_cross"}) + xr_dataset.owiAncillaryNrcs_cross.attrs["units"] = "m^2 / m^2" + xr_dataset.owiAncillaryNrcs_cross.attrs[ + "long_name"] = f"Ancillary Normalized Radar Cross Section - simulated from {config['l2_params']['crosspol_gmf']} & ancillary wind" + xr_dataset.owiLon.attrs["units"] = "degrees_east" xr_dataset.owiLon.attrs["long_name"] = "Longitude at wind cell center" xr_dataset.owiLon.attrs["standard_name"] = "longitude" @@ -343,23 +352,25 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte xr_dataset.owiElevationAngle.attrs["long_name"] = "Elevation angle at wind cell center" xr_dataset.owiElevationAngle.attrs["standard_name"] = "elevation" - xr_dataset['owiNrcs'] = xr_dataset['sigma0_ocean'].sel(pol=copol) + xr_dataset['owiNrcs'] = xr_dataset['sigma0_ocean'].sel( + pol=config["l2_params"]["copol"]) xr_dataset.owiNrcs.attrs = xr_dataset.sigma0_ocean.attrs xr_dataset.owiNrcs.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs.attrs['long_name'] = 'Normalized Radar Cross Section' xr_dataset.owiNrcs.attrs['definition'] = 'owiNrcs_no_noise_correction - owiNesz' - xr_dataset['owiMask_Nrcs'] = xr_dataset['sigma0_mask'].sel(pol=copol) + xr_dataset['owiMask_Nrcs'] = xr_dataset['sigma0_mask'].sel( + pol=config["l2_params"]["copol"]) xr_dataset.owiMask_Nrcs.attrs = xr_dataset.sigma0_mask.attrs # NESZ & DSIG xr_dataset = xr_dataset.assign( - owiNesz=(['line', 'sample'], xr_dataset.nesz.sel(pol=copol).values)) + owiNesz=(['line', 'sample'], xr_dataset.nesz.sel(pol=config["l2_params"]["copol"]).values)) xr_dataset.owiNesz.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNesz.attrs['long_name'] = 'Noise Equivalent SigmaNaught' xr_dataset['owiNrcs_no_noise_correction'] = xr_dataset['sigma0_ocean_raw'].sel( - pol=copol) + pol=config["l2_params"]["copol"]) xr_dataset.owiNrcs_no_noise_correction.attrs = xr_dataset.sigma0_ocean_raw.attrs xr_dataset.owiNrcs_no_noise_correction.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs_no_noise_correction.attrs[ @@ -378,7 +389,7 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte # sigma0_raw__corrected cross if "sigma0_raw__corrected" in xr_dataset: xr_dataset['owiNrcs_no_noise_correction_recalibrated'] = xr_dataset['sigma0_raw__corrected'].sel( - pol=copol) + pol=config["l2_params"]["copol"]) xr_dataset.owiNrcs_no_noise_correction_recalibrated.attrs = xr_dataset.sigma0_raw__corrected.attrs xr_dataset.owiNrcs_no_noise_correction_recalibrated.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs_no_noise_correction_recalibrated.attrs[ @@ -388,7 +399,7 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte xr_dataset.owiNrcs.attrs['definition'] = 'owiNrcs_no_noise_correction_recalibrated - owiNesz' - if dual_pol: + if config["l2_params"]["dual_pol"]: xr_dataset = xr_dataset.rename({ 'dsig_cross': 'owiDsig_cross', @@ -399,31 +410,31 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte 'sigma0_detrend_cross': 'owiNrcs_detrend_cross' }) - if apply_flattening: + if config["apply_flattening"]: xr_dataset = xr_dataset.rename({ 'nesz_cross_flattened': 'owiNesz_cross_flattened', }) # nrcs cross xr_dataset['owiNrcs_cross'] = xr_dataset['sigma0_ocean'].sel( - pol=crosspol) + pol=config["l2_params"]["crosspol"]) xr_dataset.owiNrcs_cross.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs_cross.attrs['long_name'] = 'Normalized Radar Cross Section' xr_dataset.owiNrcs_cross.attrs['definition'] = 'owiNrcs_cross_no_noise_correction - owiNesz_cross' xr_dataset['owiMask_Nrcs_cross'] = xr_dataset['sigma0_mask'].sel( - pol=crosspol) + pol=config["l2_params"]["crosspol"]) xr_dataset.owiMask_Nrcs_cross.attrs = xr_dataset.sigma0_mask.attrs # nesz cross xr_dataset = xr_dataset.assign(owiNesz_cross=( - ['line', 'sample'], xr_dataset.nesz.sel(pol=crosspol).values)) # no flattening + ['line', 'sample'], xr_dataset.nesz.sel(pol=config["l2_params"]["crosspol"]).values)) # no flattening xr_dataset.owiNesz_cross.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNesz_cross.attrs['long_name'] = 'Noise Equivalent SigmaNaught' xr_dataset['owiNrcs_cross_no_noise_correction'] = xr_dataset['sigma0_ocean_raw'].sel( - pol=crosspol) + pol=config["l2_params"]["crosspol"]) xr_dataset.owiNrcs_cross_no_noise_correction.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs_cross_no_noise_correction.attrs[ @@ -432,7 +443,7 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte #  sigma0_raw__corrected cross if "sigma0_raw__corrected" in xr_dataset: xr_dataset['owiNrcs_cross_no_noise_correction_recalibrated'] = xr_dataset['sigma0_raw__corrected'].sel( - pol=crosspol) + pol=config["l2_params"]["crosspol"]) xr_dataset.owiNrcs_cross_no_noise_correction_recalibrated.attrs = xr_dataset.sigma0_raw__corrected.attrs xr_dataset.owiNrcs_cross_no_noise_correction_recalibrated.attrs['units'] = 'm^2 / m^2' xr_dataset.owiNrcs_cross_no_noise_correction_recalibrated.attrs[ @@ -442,10 +453,19 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte xr_dataset.owiNrcs_cross.attrs['definition'] = 'owiNrcs_cross_no_noise_correction_recalibrated - owiNesz_cross' - if add_streaks: + if config["add_gradientsfeatures"]: xr_dataset = xr_dataset.rename({ - 'streaks_direction': 'owiStreaksDirection', + 'heterogeneity_mask': 'owiWindFilter' }) + else: + xr_dataset['owiWindFilter'] = xr.full_like(xr_dataset.owiNrcs, 0) + xr_dataset['owiWindFilter'].attrs['long_name'] = "Quality flag taking into account the local heterogeneity" + xr_dataset['owiWindFilter'].attrs['valid_range'] = np.array([0, 3]) + xr_dataset['owiWindFilter'].attrs['flag_values'] = np.array([ + 0, 1, 2, 3]) + xr_dataset['owiWindFilter'].attrs[ + 'flag_meanings'] = "homogeneous_NRCS, heterogeneous_from_co-polarization_NRCS, heterogeneous_from_cross-polarization_NRCS, heterogeneous_from_dual-polarization_NRCS" + xr_dataset['owiWindFilter'].attrs['comment'] = 'NOT COMPUTED YET' #  other variables @@ -458,15 +478,6 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte xr_dataset['owiWindQuality'].attrs['flag_meanings'] = "good medium low poor" xr_dataset['owiWindQuality'].attrs['comment'] = 'NOT COMPUTED YET' - xr_dataset['owiWindFilter'] = xr.full_like(xr_dataset.owiNrcs, 0) - xr_dataset['owiWindFilter'].attrs['long_name'] = "Quality flag taking into account the local heterogeneity" - xr_dataset['owiWindFilter'].attrs['valid_range'] = np.array([0, 3]) - xr_dataset['owiWindFilter'].attrs['flag_values'] = np.array([ - 0, 1, 2, 3]) - xr_dataset['owiWindFilter'].attrs[ - 'flag_meanings'] = "homogeneous_NRCS, heterogeneous_from_co-polarization_NRCS, heterogeneous_from_cross-polarization_NRCS, heterogeneous_from_dual-polarization_NRCS" - xr_dataset['owiWindFilter'].attrs['comment'] = 'NOT COMPUTED YET' - xr_dataset = xr_dataset.rename( {"line": "owiAzSize", "sample": "owiRaSize"}) @@ -476,8 +487,6 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte xr_dataset = xr_dataset.drop_vars(["sigma0_raw__corrected"]) xr_dataset = xr_dataset.drop_dims(['pol']) - xr_dataset.compute() - table_fillValue = { "owiWindQuality": -1, "owiHeading": 9999.99, @@ -505,7 +514,7 @@ def makeL2asOwi(xr_dataset, dual_pol, copol, crosspol, add_streaks, apply_flatte return xr_dataset, encoding -def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False, resolution='1000m'): +def preprocess(filename, outdir, config_path, overwrite=False, add_gradientsfeatures=False, resolution='1000m'): """ Main function to generate L2 product. @@ -549,6 +558,10 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False recalibration = config["recalibration"] meta = fct_meta(filename) + # si une des deux n'est pas VV VH HH HV on ne fait rien + if not all([pol in ["VV", "VH", "HH", "HV"] for pol in meta.pols.split(' ')]): + raise ValueError(f"Polarisation non gérée : meta.pols = {meta.pols}") + no_subdir_cfg = config_base.get("no_subdir", False) config["no_subdir"] = no_subdir_cfg @@ -560,6 +573,22 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False f'Using meteorological convention because "winddir_convention" was not found in config.') config["winddir_convention"] = winddir_convention + if "add_gradientsfeatures" in config_base: + add_gradientsfeatures = config_base["add_gradientsfeatures"] + else: + add_gradientsfeatures = False + logging.warning( + f'Not computing gradients by default') + config["add_gradientsfeatures"] = add_gradientsfeatures + + if "add_nrcs_model" in config_base: + add_nrcs_model = config_base["add_nrcs_model"] + else: + add_nrcs_model = False + logging.warning( + f'Not computing nrcs from model by default') + config["add_nrcs_model"] = add_nrcs_model + # creating a dictionnary of parameters config["l2_params"] = {} @@ -607,6 +636,9 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False logging.error(e) sys.exit(-1) + # load + xr_dataset = xr_dataset.load() + # defining dual_pol, and gmfs by channel if len(xr_dataset.pol.values) == 2: dual_pol = True @@ -707,7 +739,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False 90. - np.rad2deg(np.arctan2(xr_dataset.model_V10, xr_dataset.model_U10)) + 180) % 360 xr_dataset['ancillary_wind_direction'] = xr.where(xr_dataset['mask'], np.nan, - xr_dataset['ancillary_wind_direction'].compute()).transpose( + xr_dataset['ancillary_wind_direction']).transpose( *xr_dataset['ancillary_wind_direction'].dims) xr_dataset['ancillary_wind_direction'].attrs = {} xr_dataset['ancillary_wind_direction'].attrs['units'] = 'degrees_north' @@ -718,7 +750,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset['ancillary_wind_speed'] = np.sqrt( xr_dataset['model_U10']**2+xr_dataset['model_V10']**2) xr_dataset['ancillary_wind_speed'] = xr.where(xr_dataset['mask'], np.nan, - xr_dataset['ancillary_wind_speed'].compute()).transpose( + xr_dataset['ancillary_wind_speed']).transpose( *xr_dataset['ancillary_wind_speed'].dims) xr_dataset['ancillary_wind_speed'].attrs = {} xr_dataset['ancillary_wind_speed'].attrs['units'] = 'm s^-1' @@ -727,7 +759,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset['ancillary_wind_speed'].attrs['standart_name'] = 'wind_speed' xr_dataset['ancillary_wind'] = xr.where(xr_dataset['mask'], np.nan, - (xr_dataset.ancillary_wind_speed * np.exp(1j * xsarsea.dir_meteo_to_sample(xr_dataset.ancillary_wind_direction, xr_dataset.ground_heading))).compute()).transpose( + (xr_dataset.ancillary_wind_speed * np.exp(1j * xsarsea.dir_meteo_to_sample(xr_dataset.ancillary_wind_direction, xr_dataset.ground_heading)))).transpose( *xr_dataset['ancillary_wind_speed'].dims) xr_dataset.attrs['ancillary_source'] = xr_dataset['model_U10'].attrs['history'].split('decoded: ')[ @@ -736,7 +768,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False # nrcs processing xr_dataset['sigma0_ocean'] = xr.where(xr_dataset['mask'], np.nan, - xr_dataset['sigma0'].compute()).transpose(*xr_dataset['sigma0'].dims) + xr_dataset['sigma0']).transpose(*xr_dataset['sigma0'].dims) xr_dataset['sigma0_ocean'].attrs = xr_dataset['sigma0'].attrs #  we forced it to 1e-15 xr_dataset['sigma0_ocean'].attrs['comment'] = "clipped, no values <=0 ; 1e-15 instread" @@ -751,7 +783,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset['sigma0_ocean'] <= 0, 1e-15, xr_dataset['sigma0_ocean']) xr_dataset['sigma0_ocean_raw'] = xr.where(xr_dataset['mask'], np.nan, - xr_dataset['sigma0_raw'].compute()).transpose(*xr_dataset['sigma0_raw'].dims) + xr_dataset['sigma0_raw']).transpose(*xr_dataset['sigma0_raw'].dims) xr_dataset['sigma0_ocean_raw'].attrs = xr_dataset['sigma0_raw'].attrs @@ -765,7 +797,7 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset.sigma0.sel(pol=crosspol), xr_dataset.incidence, model=model_cross) if config["apply_flattening"]: xr_dataset = xr_dataset.assign(nesz_cross_flattened=( - ['line', 'sample'], windspeed.nesz_flattening(xr_dataset.nesz.sel(pol=crosspol), xr_dataset.incidence))) + ['line', 'sample'], windspeed.nesz_flattening(xr_dataset.nesz.sel(pol=crosspol), xr_dataset.incidence).data)) xr_dataset['nesz_cross_flattened'].attrs[ "comment"] = 'nesz has been flattened using windspeed.nesz_flattening' # dsig @@ -793,11 +825,58 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset.attrs["path_aux_cal_old"] = os.path.basename(os.path.dirname( os.path.dirname(xsar_dataset.datatree['recalibration'].attrs['path_aux_cal_old']))) - if add_streaks: + if add_nrcs_model: + # add timing + phi = np.abs( + np.rad2deg(xsarsea.dir_meteo_to_sample( + xr_dataset["ancillary_wind_direction"], xr_dataset["ground_heading"])) + ) + + varnames = ["ancillary_nrcs"] + gmf_names = [model_co] + if dual_pol: + varnames.append("ancillary_nrcs_cross") + gmf_names.append(model_cross) + + for idx, gmf_name in enumerate(gmf_names): + + @timing(logger=logger.info) + def apply_lut_to_dataset(): + lut = xsarsea.windspeed.get_model( + gmf_name).to_lut(unit="linear") + + def lut_selection(incidence, wspd, phi): + if "phi" in lut.coords: + return lut.sel( + incidence=incidence, wspd=wspd, phi=phi, method="nearest" + ) + else: + return lut.sel( + incidence=incidence, wspd=wspd, method="nearest" + ) + + xr_dataset[varnames[idx]] = xr.apply_ufunc( + lut_selection, + xr_dataset.incidence, + xr_dataset.ancillary_wind_speed, + phi, + vectorize=True, + dask="parallelized", + output_dtypes=[float], + ) + + apply_lut_to_dataset() + + if add_gradientsfeatures: + from grdwindinversion.streaks import GradientFeatures + xsar_dataset_100 = fct_dataset( meta, resolution='100m') + xr_dataset_100 = xsar_dataset_100.datatree['measurement'].to_dataset() xr_dataset_100 = xr_dataset_100.rename(map_model) + # load dataset + xr_dataset_100 = xr_dataset_100.load() # adding sigma0 detrend xr_dataset_100['sigma0_detrend'] = xsarsea.sigma0_detrend( @@ -819,15 +898,40 @@ def preprocess(filename, outdir, config_path, overwrite=False, add_streaks=False xr_dataset_100.land_mask.values = binary_dilation(xr_dataset_100['land_mask'].values.astype('uint8'), structure=np.ones((3, 3), np.uint8), iterations=3) xr_dataset_100['sigma0_detrend'] = xr.where( - xr_dataset_100['land_mask'], np.nan, xr_dataset_100['sigma0'].compute()).transpose(*xr_dataset_100['sigma0'].dims) - - xr_dataset['streaks_direction'] = get_streaks( - xr_dataset, xr_dataset_100) + xr_dataset_100['land_mask'], np.nan, xr_dataset_100['sigma0']).transpose(*xr_dataset_100['sigma0'].dims) + + xr_dataset_100['ancillary_wind'] = ( + xr_dataset_100.model_U10 + 1j * xr_dataset_100.model_V10) * np.exp(1j * np.deg2rad(xr_dataset_100.ground_heading)) + + downscales_factors = [1, 2, 4, 8] + # 4 and 8 must be in downscales_factors + assert all([x in downscales_factors for x in [4, 8]]) + + gradientFeatures = GradientFeatures( + xr_dataset=xr_dataset, + xr_dataset_100=xr_dataset_100, + windows_sizes=[1600, 3200], + downscales_factors=downscales_factors, + window_step=1 + ) + dataArraysHeterogeneity = gradientFeatures.get_heterogeneity_mask( + config) + xr_dataset = xr_dataset.merge(dataArraysHeterogeneity) + + # Create your streaks dataset + xr_dataset_streaks = xr.Dataset({ + 'dir_smooth': gradientFeatures.streaks_individual(), + 'dir_mean_smooth': gradientFeatures.streaks_mean_smooth(), + 'dir_smooth_mean': gradientFeatures.streaks_smooth_mean(), + }) + else: + xr_dataset_streaks = None - return xr_dataset, out_file, config + return xr_dataset, xr_dataset_streaks, out_file, config -def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add_streaks=False, resolution='1000m'): +@timing(logger=logger.info) +def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, resolution='1000m'): """ Main function to generate L2 product. @@ -854,8 +958,8 @@ def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add final dataset """ - xr_dataset, out_file, config = preprocess( - filename, outdir, config_path, overwrite, add_streaks, resolution) + xr_dataset, xr_dataset_streaks, out_file, config = preprocess( + filename, outdir, config_path, overwrite, resolution) model_co = config["l2_params"]["model_co"] model_cross = config["l2_params"]["model_cross"] @@ -899,8 +1003,6 @@ def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add model_co=model_co, model_cross=model_cross, ** kwargs) - wind_co.compute() - # windspeed_co xr_dataset['windspeed_co'] = np.abs(wind_co) xr_dataset["windspeed_co"].attrs["units"] = "m.s⁻1" @@ -917,7 +1019,6 @@ def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add # windspeed_dual / windspeed_cr / /winddir_dual / winddir_cr if dual_pol and wind_dual is not None: - wind_dual.compute() xr_dataset['windspeed_dual'] = np.abs(wind_dual) xr_dataset["windspeed_dual"].attrs["units"] = "m.s⁻1" xr_dataset["windspeed_dual"].attrs["long_name"] = "Wind speed inverted from model %s (%s) & %s (%s)" % ( @@ -951,8 +1052,9 @@ def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add "long_name"] = f"{ancillary_name} wind direction in oceanographic convention (clockwise, to), ex: 0°=to north, 90°=to east" xr_dataset, encoding = makeL2asOwi( - xr_dataset, dual_pol, copol, crosspol, add_streaks=add_streaks, apply_flattening=config["apply_flattening"]) + xr_dataset, config) + xr_dataset = xr_dataset.compute() #  add attributes firstMeasurementTime = None lastMeasurementTime = None @@ -1035,7 +1137,26 @@ def makeL2(filename, outdir, config_path, overwrite=False, generateCSV=True, add os.makedirs(os.path.dirname(out_file), exist_ok=True) + # Sauvegarde de xr_dataset dans le fichier de sortie final xr_dataset.to_netcdf(out_file, mode="w", encoding=encoding) + + # Vérifier si le dataset de streaks est présent + if xr_dataset_streaks is not None: + # Créer un fichier temporaire pour le dataset streaks + with tempfile.NamedTemporaryFile(delete=False, suffix=".nc") as tmp_file: + temp_out_file = tmp_file.name + + # Écrire xr_dataset_streaks dans le fichier temporaire + xr_dataset_streaks.to_netcdf( + temp_out_file, mode="w", group="owiWindStreaks") + + # Charger le fichier temporaire et l'ajouter au fichier final en tant que groupe + with xr.open_dataset(temp_out_file, group="owiWindStreaks") as ds_streaks: + ds_streaks.to_netcdf(out_file, mode="a", group="owiWindStreaks") + + # Supprimer le fichier temporaire après l'opération + os.remove(temp_out_file) + if generateCSV: df = xr_dataset.to_dataframe() df = df[df.owiMask == False] diff --git a/grdwindinversion/streaks.py b/grdwindinversion/streaks.py index 9588aa4..117b8de 100644 --- a/grdwindinversion/streaks.py +++ b/grdwindinversion/streaks.py @@ -1,79 +1,400 @@ -import xarray as xr import xsarsea.gradients +import cv2 +import xarray as xr import xarray as xr from scipy.ndimage import binary_dilation import numpy as np +import logging + +import logging +logger = logging.getLogger('grdwindinversion.streaks') +logger.addHandler(logging.NullHandler()) + + +class GradientFeatures: + def __init__(self, xr_dataset, xr_dataset_100, windows_sizes, downscales_factors, window_step=1): + """ + Initialize the GradientFeatures with the dataset and parameters. + + Parameters: + - xr_dataset: xarray.Dataset containing the SAR data. + - xr_dataset_100: xarray.Dataset containing the 100m resolution SAR data. + - windows_sizes: List of window sizes for gradient computation. + - downscales_factors: List of downscale factors for gradient computation. + - window_step: Step size for the window (default is 1). + """ + self.xr_dataset = xr_dataset + self.xr_dataset_100 = xr_dataset_100 + self.windows_sizes = windows_sizes + self.downscales_factors = downscales_factors + self.window_step = window_step + self.gradients = None + self.hist = None + self._compute_gradients() + + def _compute_gradients(self): + """ + Compute gradients and histograms. + + Parameters: + - xr_dataset: xarray.Dataset containing the SAR data. + - windows_sizes: List of window sizes for gradient computation. + - downscales_factors: List of downscale factors for gradient computation. + + """ + self.gradients = xsarsea.gradients.Gradients( + self.xr_dataset_100['sigma0_detrend'], + windows_sizes=self.windows_sizes, + downscales_factors=self.downscales_factors, + window_step=self.window_step + ) + self.hist = self.gradients.histogram + # Get orthogonal gradients + self.hist['angles'] = self.hist['angles'] + np.pi / 2 + + def get_heterogeneity_mask(self, config): + """ + Get the heterogeneity mask from the wind field. + """ + dual_pol = config["l2_params"]["dual_pol"] + + new_dataArrays = {} + + try: + + sigma0_400_co = [da.sigma0 for da in self.gradients.gradients_list if ( + da.sigma0["pol"] == config["l2_params"]["copol"] and da.sigma0.downscale_factor == 4)][0] + sigs = [sigma0_400_co] + + if dual_pol: + sigma0_800_cross = [da.sigma0 for da in self.gradients.gradients_list if ( + da.sigma0["pol"] == config["l2_params"]["crosspol"] and da.sigma0.downscale_factor == 8)][0] + sigs.append(sigma0_800_cross) + + filters = {} + for sig in sigs: + + pol = sig["pol"].values + res = 100 * sig.downscale_factor.values + + # delete useless coords : could be problematic to have it later + if 'downscale_factor' in sig.coords: + sig = sig.reset_coords("downscale_factor", drop=True) + + if 'window_size' in sig.coords: + sig = sig.reset_coords("window_size", drop=True) + # mask + sig = xr.where(sig <= 0, 1e-15, sig) + + # map incidence for detrend + incidence = xr.DataArray(data=cv2.resize( + self.xr_dataset_100.incidence.values, sig.shape[::-1], cv2.INTER_NEAREST), dims=sig.dims, coords=sig.coords) + + sigma0_detrend = xsarsea.sigma0_detrend(sig, incidence) + + filter_name = str(res)+"_"+str(pol) + I = sigma0_detrend + f1, f2, f3, f4, f = xsarsea.gradients.filtering_parameters(I) + filters[filter_name] = f + + thresholds = [0.78] # < is unusable + if dual_pol: + # Seuil pour crosspol si dual_pol est activé + thresholds.append(0.71) + + for idx_filter, filter in enumerate(filters): + # interp to user resolution and map on dataset grid + new_dataArrays[filter] = filters[filter].interp( + line=self.xr_dataset.line, sample=self.xr_dataset.sample, method="nearest") + new_dataArrays[filter+"_mask"] = xr.where( + new_dataArrays[filter] > thresholds[idx_filter], True, False) + + varname_400_copol_mask = f'400_{config["l2_params"]["copol"]}_mask' + varname_800_crosspol_mask = f'800_{config["l2_params"]["crosspol"]}_mask' + + # Cas 0 : no heterogeneity + new_dataArrays["heterogeneity_mask"] = xr.full_like( + new_dataArrays[varname_400_copol_mask], 0) + + if dual_pol: + # Cas 3 : Dual-polarization + new_dataArrays["heterogeneity_mask"] = xr.where( + new_dataArrays[varname_400_copol_mask] & new_dataArrays[varname_800_crosspol_mask], 3, new_dataArrays["heterogeneity_mask"]) + + # Cas 1 : Co-polarization only + new_dataArrays["heterogeneity_mask"] = xr.where( + new_dataArrays[varname_400_copol_mask] & ~new_dataArrays[varname_800_crosspol_mask], 1, new_dataArrays["heterogeneity_mask"]) + + # Cas 2 : Cross-polarization only + new_dataArrays["heterogeneity_mask"] = xr.where( + ~new_dataArrays[varname_400_copol_mask] & new_dataArrays[varname_800_crosspol_mask], 2, new_dataArrays["heterogeneity_mask"]) + + # Attributes + new_dataArrays["heterogeneity_mask"].attrs["valid_range"] = np.array([ + 0, 3]) + new_dataArrays["heterogeneity_mask"].attrs["flag_values"] = np.array([ + 0, 1, 2, 3]) + new_dataArrays["heterogeneity_mask"].attrs["flag_meanings"] = ( + "homogeneous_NRCS, heterogeneous_from_co-polarization_NRCS, " + "heterogeneous_from_cross-polarization_NRCS, heterogeneous_from_dual-polarization_NRCS" + ) + else: + # no crosspol + new_dataArrays["heterogeneity_mask"] = xr.where( + new_dataArrays[varname_400_copol_mask], 1, new_dataArrays["heterogeneity_mask"]) + + # Attributs pour le cas single-pol + new_dataArrays["heterogeneity_mask"].attrs["valid_range"] = np.array([ + 0, 1]) + new_dataArrays["heterogeneity_mask"].attrs["flag_values"] = np.array([ + 0, 1]) + new_dataArrays["heterogeneity_mask"].attrs["flag_meanings"] = ( + "homogeneous_NRCS, heterogeneous_from_co-polarization_NRCS" + ) + + # Attributs généraux + new_dataArrays["heterogeneity_mask"].attrs["long_name"] = "Quality flag taking into account the local heterogeneity" + return new_dataArrays + + except Exception as e: + logging.error("Error in get_heterogeneity_mask: %s", e) + + new_dataArrays["heterogeneity_mask"] = xr.DataArray(data=np.nan * np.ones([len(self.xr_dataset.coords[dim]) for dim in ['line', 'sample']]), + dims=[ + 'line', 'sample'], + coords=[self.xr_dataset.coords[dim] + for dim in ['line', 'sample']], + attrs={"comment": "no heterogeneity mask found"}) + + return new_dataArrays + + def _interpolate_and_mask(self, streaks, dataset): + """Interpolate streaks to match the xr_dataset grid and apply land mask.""" + # Interpolate to match the original dataset grid + streaks_dir_interp = streaks['angle'].interp( + line=dataset.line, + sample=dataset.sample + ) + # Apply land mask + streaks_dir_interp = xr.where( + dataset['land_mask'], np.nan, streaks_dir_interp + ) + return streaks_dir_interp + + def _remove_ambiguity(self, streaks, dataset_hr): + """ + Remove direction ambiguity using ancillary wind data. + + Parameters: + - streaks: xarray.Dataset containing the streaks. + - dataset_hr: xarray.Dataset containing the ancillary wind data. + """ + + # Select ancillary wind data matching the streaks coordinates + + try: + + ancillary_wind = dataset_hr['ancillary_wind'].sel( + line=streaks.line, + sample=streaks.sample, + method='nearest' + ).compute() + + except Exception as e: + logging.warn( + "Warning in _remove_ambiguity\nProbably because 1 is not in downscale factor and so lines are not aligned. interpolating ; error : %s", e) + + ancillary_wind = dataset_hr['ancillary_wind'].interp( + line=streaks.line, + sample=streaks.sample, + method='nearest' + ).compute() + + return streaks + + # Convert angles to complex numbers + streaks_c = streaks['weight'] * np.exp(1j * streaks['angle']) + # Calculate the difference in angle + diff_angle = xr.apply_ufunc(np.angle, ancillary_wind / streaks_c) + + # Remove ambiguity + streaks_c = xr.where(np.abs(diff_angle) > np.pi / + 2, -streaks_c, streaks_c) + + # Update streaks with corrected values + streaks['weight'] = np.abs(streaks_c) + streaks['angle'] = xr.apply_ufunc(np.angle, streaks_c) + return streaks + + def streaks_smooth_mean(self): + """Smooth the histograms individually, compute the mean""" + + try: + # Smooth the histograms individually + hist_smooth = self.hist.copy() + hist_smooth['weight'] = xsarsea.gradients.circ_smooth( + hist_smooth['weight']) + + # Compute the mean across 'downscale_factor', 'window_size', and 'pol' + hist_smooth_mean = hist_smooth.mean( + ['downscale_factor', 'window_size', 'pol']) + + # Select histogram peak + iangle_smooth_mean = hist_smooth_mean['weight'].fillna( + 0).argmax(dim='angles') + streaks_dir_smooth_mean = hist_smooth_mean['angles'].isel( + angles=iangle_smooth_mean) + streaks_weight_smooth_mean = hist_smooth_mean['weight'].isel( + angles=iangle_smooth_mean) + + # Combine angles and weights into a dataset + streaks_smooth_mean = xr.Dataset({ + 'angle': streaks_dir_smooth_mean, + 'weight': streaks_weight_smooth_mean + }) + + # Remove 'angles' coordinate + streaks_smooth_mean = streaks_smooth_mean.reset_coords( + 'angles', drop=True) + + # Remove ambiguity with ancillary wind + streaks_smooth_mean = self._remove_ambiguity( + streaks_smooth_mean, self.xr_dataset_100) + + # Interpolate to match the original dataset grid + streaks_dir_smooth_mean_interp = self._interpolate_and_mask( + streaks_smooth_mean, self.xr_dataset) + + # Set attributes + streaks_dir_smooth_mean_interp.attrs['comment'] = ( + 'Angle in radians, anticlockwise from 0=line; smoothed then mean solution' + ) + streaks_dir_smooth_mean_interp.attrs['description'] = ( + 'Wind direction estimated from local gradient; histograms smoothed first, then mean computed' + ) + + return streaks_dir_smooth_mean_interp + + except Exception as e: + logging.error("Error in streaks_smooth_mean: %s", e) + + streaks_dir_smooth_mean_interp = xr.DataArray(data=np.nan * np.ones([len(self.xr_dataset.coords[dim]) for dim in ['line', 'sample']]), + dims=[ + 'line', 'sample'], + coords=[self.xr_dataset.coords[dim] + for dim in ['line', 'sample']], + attrs={"comment": "no streaks_smooth_mean found"}) + + return streaks_dir_smooth_mean_interp + + def streaks_mean_smooth(self): + """Compute the mean of the histograms and then smooth the mean histogram.""" + try: + # Compute the mean of the histograms + hist_mean = self.hist.copy().mean( + ['downscale_factor', 'window_size', 'pol']) + + # Smooth the mean histogram + hist_mean_smooth = hist_mean.copy() + hist_mean_smooth['weight'] = xsarsea.gradients.circ_smooth( + hist_mean['weight']) + + # Select histogram peak + iangle_mean_smooth = hist_mean_smooth['weight'].fillna( + 0).argmax(dim='angles') + streaks_dir_mean_smooth = hist_mean_smooth['angles'].isel( + angles=iangle_mean_smooth) + streaks_weight_mean_smooth = hist_mean_smooth['weight'].isel( + angles=iangle_mean_smooth) + + # Combine angles and weights into a dataset + streaks_mean_smooth = xr.Dataset({ + 'angle': streaks_dir_mean_smooth, + 'weight': streaks_weight_mean_smooth + }) + + # Remove 'angles' coordinate + streaks_mean_smooth = streaks_mean_smooth.reset_coords( + 'angles', drop=True) + + # Remove ambiguity with ancillary wind + streaks_mean_smooth = self._remove_ambiguity( + streaks_mean_smooth, self.xr_dataset_100) + # Interpolate to match the original dataset grid + streaks_dir_mean_smooth_interp = self._interpolate_and_mask( + streaks_mean_smooth, self.xr_dataset) + + # Set attributes + streaks_dir_mean_smooth_interp.attrs['comment'] = ( + 'Angle in radians, anticlockwise from 0=line; mean then smoothed solution' + ) + streaks_dir_mean_smooth_interp.attrs['description'] = ( + 'Wind direction estimated from local gradient; mean computed first, then histogram smoothed' + ) + + return streaks_dir_mean_smooth_interp + except Exception as e: + logging.error("Error in streaks_mean_smooth: %s", e) + + streaks_dir_mean_smooth_interp = xr.DataArray(data=np.nan * np.ones([len(self.xr_dataset.coords[dim]) for dim in ['line', 'sample']]), + dims=[ + 'line', 'sample'], + coords=[self.xr_dataset.coords[dim] + for dim in ['line', 'sample']], + attrs={"comment": "no streaks_mean_smooth found"}) + + return streaks_dir_mean_smooth_interp + + def streaks_individual(self): + """Smooth the histograms individually. Only smooth the histograms, do not compute the mean.""" + + try: + # Compute the mean of the histograms + hist_smooth = self.hist.copy() + hist_smooth['weight'] = xsarsea.gradients.circ_smooth( + hist_smooth['weight']) + + # Select histogram peak for each individual solution + iangle_individual = hist_smooth['weight'].fillna( + 0).argmax(dim='angles') + streaks_dir_individual = hist_smooth['angles'].isel( + angles=iangle_individual) + streaks_weight_individual = hist_smooth['weight'].isel( + angles=iangle_individual) + # Combine angles and weights into a dataset + streaks_individual = xr.Dataset({ + 'angle': streaks_dir_individual, + 'weight': streaks_weight_individual + }) + # Remove 'angles' coordinate + streaks_individual = streaks_individual.reset_coords( + 'angles', drop=True) + + # Remove ambiguity with ancillary wind for each individual solution + streaks_individual = self._remove_ambiguity( + streaks_individual, self.xr_dataset_100) + + # Interpolate to match the original dataset grid + streaks_dir_individual_interp = self._interpolate_and_mask( + streaks_individual, self.xr_dataset) + + # Set attributes + streaks_dir_individual_interp.attrs['comment'] = ( + 'Angle in radians, anticlockwise from 0=line; individual solutions' + ) + streaks_dir_individual_interp.attrs['description'] = ( + 'Wind direction estimated from local gradient for each individual solution; histograms smoothed individually' + ) + return streaks_dir_individual_interp + + except Exception as e: + logging.error("Error in streaks_individual: %s", e) + streaks_dir_individual_interp = xr.DataArray(data=np.nan * np.ones([len(self.xr_dataset.coords[dim]) for dim in ['line', 'sample']]), + dims=[ + 'line', 'sample'], + coords=[self.xr_dataset.coords[dim] + for dim in ['line', 'sample']], + attrs={"comment": "no streaks_individual found"}) -def get_streaks(xr_dataset, xr_dataset_100): - """ - Get the streaks from the wind field. - - Parameters - ---------- - xr_dataset : xarray.Dataset - dataset at user resolution. - xr_dataset_100 : xarray.Dataset - dataset at 100m resolution. - - Returns - ------- - xarray.Dataset - Extract wind direction from Koch Method using xsarsea tools. - """ - - # return empy dataArray, waiting for solution - return xr.DataArray(data=np.nan * np.ones([len(xr_dataset.coords[dim]) for dim in ['line','sample']]), - dims=['line','sample'], - coords=[xr_dataset.coords[dim] for dim in ['line','sample']]) - # - - """ - gradients = xsarsea.gradients.Gradients(xr_dataset_100['sigma0_detrend'], windows_sizes=[ - 1600, 3200], downscales_factors=[1, 2], window_step=1) - - # get gradients histograms as an xarray dataset - hist = gradients.histogram - - # get orthogonals gradients - hist['angles'] = hist['angles'] + np.pi/2 - - # mean - hist_mean = hist.mean(['downscale_factor', 'window_size', 'pol']) - - # smooth - hist_mean_smooth = hist_mean.copy() - hist_mean_smooth['weight'] = xsarsea.gradients.circ_smooth( - hist_mean['weight']) - - # smooth only - # hist_smooth = hist.copy() - # hist_smooth['weight'] = xsarsea.gradients.circ_smooth(hist_smooth['weight']) - - # select histogram peak - iangle = hist_mean_smooth['weight'].fillna(0).argmax(dim='angles') - streaks_dir = hist_mean_smooth.angles.isel(angles=iangle) - streaks_weight = hist_mean_smooth['weight'].isel(angles=iangle) - streaks = xr.merge( - [dict(angle=streaks_dir, weight=streaks_weight)]).drop('angles') - - # streaks are [0, pi]. Remove ambiguity with anciallary wind - ancillary_wind = xr_dataset_100['ancillary_wind'].sel(line=streaks.line, - sample=streaks.sample, - method='nearest').compute() - streaks_c = streaks['weight'] * np.exp(1j * streaks['angle']) - diff_angle = xr.apply_ufunc(np.angle, ancillary_wind / streaks_c) - streaks_c = xr.where(np.abs(diff_angle) > np.pi/2, -streaks_c, streaks_c) - streaks['weight'] = np.abs(streaks_c) - streaks['angle'] = xr.apply_ufunc(np.angle, streaks_c) - - streaks_dir = xr.apply_ufunc( - np.angle, streaks_c.interp(line=xr_dataset.line, sample=xr_dataset.sample)) - streaks_dir = xr.where( - xr_dataset['land_mask'], np.nan, streaks_dir) - streaks_dir.attrs['comment'] = 'angle in radians, anticlockwise, 0=line' - streaks_dir.attrs['description'] = 'wind direction estimated from local gradient, and direction ambiguity removed with ancillary wind' - - return streaks_dir - - """ + return streaks_dir_individual_interp diff --git a/grdwindinversion/utils.py b/grdwindinversion/utils.py index 014ee0c..f986a69 100644 --- a/grdwindinversion/utils.py +++ b/grdwindinversion/utils.py @@ -1,7 +1,22 @@ +import os +import time import logging import xsarsea +logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(name)s - %(levelname)s - %(message)s') +logger = logging.getLogger('grdwindinversion') + + +mem_monitor = True +try: + from psutil import Process +except ImportError: + logger.warning("psutil module not found. Disabling memory monitor") + mem_monitor = False + + def check_incidence_range(incidence, models, **kwargs): """ Check if the incidence range of the dataset is within the range of the LUT of the model. @@ -81,3 +96,28 @@ def check_format(s): return "not_written_in_lut" else: return '/' + + +def timing(logger=logger.debug): + """provide a @timing decorator() for functions, that log time spent in it""" + + def decorator(f): + # @wraps(f) + def wrapper(*args, **kwargs): + mem_str = '' + process = None + if mem_monitor: + process = Process(os.getpid()) + startrss = process.memory_info().rss + starttime = time.time() + result = f(*args, **kwargs) + endtime = time.time() + if mem_monitor: + endrss = process.memory_info().rss + mem_str = 'mem: %+.1fMb' % ((endrss - startrss) / (1024 ** 2)) + logger( + 'timing %s : %.2fs. %s' % (f.__name__, endtime - starttime, mem_str)) + return result + wrapper.__doc__ = f.__doc__ + return wrapper + return decorator diff --git a/tests/config_test.yaml b/tests/config_test.yaml index 92337e8..e9fd0f5 100644 --- a/tests/config_test.yaml +++ b/tests/config_test.yaml @@ -1,5 +1,7 @@ no_subdir: True convention: "meteorological" +add_gradientsfeatures: True +add_nrcs_model: True S1A: GMF_VV_NAME: "gmf_cmod5n" GMF_VH_NAME: "gmf_s1_v2" diff --git a/tests/test_grdwindinversion_ci.py b/tests/test_grdwindinversion_ci.py index 4797a18..595de9f 100644 --- a/tests/test_grdwindinversion_ci.py +++ b/tests/test_grdwindinversion_ci.py @@ -40,7 +40,6 @@ def test_makeL2_generation(): config_path=config_path, overwrite=True, # Set to True to ensure a clean run generateCSV=False, # Disable CSV generation for now - add_streaks=False, resolution="1000m", )