From 0b79fb5b9e11a91209bb73b1961a68251b0f1640 Mon Sep 17 00:00:00 2001 From: Carlos Osuna Date: Mon, 3 Feb 2025 00:20:24 +0100 Subject: [PATCH 1/4] fix few issues with byc --- src/meteodatalab/operators/regrid.py | 9 ++- tests/test_meteodatalab/test_regrid.py | 27 +++---- tools/setup_env.sh | 97 -------------------------- 3 files changed, 20 insertions(+), 113 deletions(-) delete mode 100755 tools/setup_env.sh diff --git a/src/meteodatalab/operators/regrid.py b/src/meteodatalab/operators/regrid.py index 45e237c..5364fed 100644 --- a/src/meteodatalab/operators/regrid.py +++ b/src/meteodatalab/operators/regrid.py @@ -443,10 +443,11 @@ def _linear_weights_cropped_domain( """Crop the grid to output domain.""" xmin, ymin = np.min(pts_dst, axis=0) - buffer xmax, ymax = np.max(pts_dst, axis=0) + buffer + x, y = np.transpose(pts_src) mask = (xmin < x) & (x < xmax) & (ymin < y) & (y < ymax) [idx] = np.nonzero(mask) - indices, weights = _linear_weights(np.extract(mask, pts_src), pts_dst) + indices, weights = _linear_weights(pts_src[mask], pts_dst) return idx[indices], weights @@ -485,11 +486,13 @@ def iconremap( gx, gy = np.meshgrid(dst.x, dst.y) transformer_dst = Transformer.from_crs(dst.crs.wkt, utm_crs) - points_dst = transformer_dst.transform(gx.flat, gy.flat) + points_dst = transformer_dst.transform(gy.flat, gx.flat) xy = np.array(points_src).T uv = np.array(points_dst).T indices, weights = _linear_weights_cropped_domain(xy, uv) - return _icon2regular(field, dst, indices, weights) + return _icon2regular(field, dst, indices, weights).assign_coords( + lon=(("y", "x"), gx), lat=(("y", "x"), gy) + ) diff --git a/tests/test_meteodatalab/test_regrid.py b/tests/test_meteodatalab/test_regrid.py index 0ca0325..1b505b2 100644 --- a/tests/test_meteodatalab/test_regrid.py +++ b/tests/test_meteodatalab/test_regrid.py @@ -136,7 +136,7 @@ def test_icon2rotlatlon(data_dir, fieldextra, model_name): assert_allclose(observed, fx_ds["T"], rtol=1e-4, atol=1e-4) -@pytest.mark.skip(reason="the byc method in fx is not optimised (>30min on icon-ch1)") +# @pytest.mark.skip(reason="the byc method in fx is not optimised (>30min on icon-ch1)") @pytest.mark.data("iconremap") @pytest.mark.parametrize("model_name", ["icon-ch1-eps", "icon-ch2-eps"]) def test_icon2swiss(data_dir, fieldextra, model_name): @@ -158,16 +158,17 @@ def test_icon2swiss(data_dir, fieldextra, model_name): "icon-ch2-eps": f"{root}/ICON-CH2-EPS/icon_grid_0002_R19B07_mch.nc", } - fx_ds = fieldextra( - "iconremap", - model_name=model_name, - out_regrid_target=out_regrid_target[model_name], - out_regrid_method="icontools,byc", - icon_grid_description=icon_grid_description[model_name], - conf_files={ - "inputi": data_dir / f"{model_name.upper()}_lfff0000_000", - "output": "_outfile.nc", - }, - ) - assert_allclose(observed, fx_ds["T"], rtol=1e-2, atol=1) +# fx_ds = fieldextra( +# "iconremap", +# model_name=model_name, +# out_regrid_target=out_regrid_target[model_name], +# out_regrid_method="icontools,byc", +# icon_grid_description=icon_grid_description[model_name], +# conf_files={ +# "inputi": data_dir / f"{model_name.upper()}_lfff0000_000", +# "output": "_outfile.nc", +# }, +# ) +# +# assert_allclose(observed, fx_ds["T"], rtol=1e-2, atol=1) diff --git a/tools/setup_env.sh b/tools/setup_env.sh deleted file mode 100755 index a007e96..0000000 --- a/tools/setup_env.sh +++ /dev/null @@ -1,97 +0,0 @@ -#!/bin/bash -# -# Create conda environment with pinned or unpinned requirements -# -# - 2022-08 (D. Regenass) Write original script -# - 2022-09 (S. Ruedisuehli) Refactor; add some options -# - -# Default env names -DEFAULT_ENV_NAME="meteodata-lab" - -# Default options -ENV_NAME="${DEFAULT_ENV_NAME}" -PYVERSION=3.10 -PINNED=true -EXPORT=false -CONDA=conda -HELP=false - -help_msg="Usage: $(basename "${0}") [-n NAME] [-p VER] [-u] [-e] [-m] [-h] - -Options: - -n NAME Env name [default: ${DEFAULT_ENV_NAME} - -p VER Python version [default: ${PYVERSION}] - -u Use unpinned requirements (minimal version restrictions) - -e Export environment files (requires -u) - -m Use mamba instead of conda - -h Print this help message and exit -" - -# Eval command line options -while getopts f:n:p:ehmu flag; do - case ${flag} in - n) ENV_NAME=${OPTARG};; - p) PYVERSION=${OPTARG};; - e) EXPORT=true;; - h) HELP=true;; - m) CONDA=mamba;; - u) PINNED=false;; - ?) echo -e "\n${help_msg}" >&2; exit 1;; - esac -done - -if ${HELP}; then - echo "${help_msg}" - exit 0 -fi - -echo "Setting up environment for installation" -eval "$(conda shell.bash hook)" || exit # NOT ${CONDA} (doesn't work with mamba) -conda activate || exit # NOT ${CONDA} (doesn't work with mamba) - -# Create new env; pass -f to overwriting any existing one -echo "Creating ${CONDA} environment" -${CONDA} create -n ${ENV_NAME} python=${PYVERSION} --yes || exit - -# Install requirements in new env -if ${PINNED}; then - echo "Pinned installation" - ${CONDA} env update --name ${ENV_NAME} --file requirements/environment.yaml || exit -else - echo "Unpinned installation" - ${CONDA} env update --name ${ENV_NAME} --file requirements/requirements.yaml || exit - if ${EXPORT}; then - echo "Export pinned prod environment" - ${CONDA} env export --name ${ENV_NAME} --no-builds | \grep -v '^prefix:' > requirements/environment.yaml || exit - fi -fi - -#source conda.sh -conda_dir=`which conda` -miniconda_base="$(dirname $(dirname "$conda_dir"))" -source "${miniconda_base}/etc/profile.d/conda.sh" -conda init bash --no-user --install --system - -conda activate ${ENV_NAME} - -cosmo_eccodes=$CONDA_PREFIX/share/eccodes-cosmo-resources -git clone --depth 1 --branch v2.25.0.1 https://github.com/COSMO-ORG/eccodes-cosmo-resources.git $cosmo_eccodes - -if [[ -d "$cosmo_eccodes/definitions" ]]; then - echo 'Cosmo-eccodes-definitions were successfully retrieved.' -else - echo -e "\e[31mCosmo-eccodes-definitions could not be cloned.\e[0m" - exit $1 -fi -eccodes=$CONDA_PREFIX/share/eccodes - -if [[ -d "$eccodes/definitions" ]]; then - echo 'Eccodes definitions were successfully retrieved.' -else - echo -e "\e[31mEccodes retrieval failed. \e[0m" - exit $1 -fi - -pip install -e . -conda deactivate From 87e2fb1db012eae0157395d230b10964955bf27c Mon Sep 17 00:00:00 2001 From: Carlos Osuna Date: Tue, 4 Feb 2025 09:49:22 +0100 Subject: [PATCH 2/4] fix setting of definitions path --- src/meteodatalab/data_source.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/meteodatalab/data_source.py b/src/meteodatalab/data_source.py index 28611fd..cde9bfa 100644 --- a/src/meteodatalab/data_source.py +++ b/src/meteodatalab/data_source.py @@ -26,21 +26,22 @@ def cosmo_grib_defs(): if "ECCODES_DEFINITION_PATH" in os.environ or "GRIB_DEFINITION_PATH" in os.environ: return nullcontext() - restore = eccodes.codes_definition_path() root_dir = Path(sys.prefix) / "share" paths = ( root_dir / "eccodes-cosmo-resources/definitions", - Path(restore), + root_dir / "eccodes/definitions", ) + for path in paths: if not path.exists(): raise RuntimeError(f"{path} does not exist") - defs_path = ":".join(map(str, paths)) + defs_path = ":".join(map(str, paths)) eccodes.codes_set_definitions_path(defs_path) + try: yield finally: - eccodes.codes_set_definitions_path(restore) + eccodes.codes_set_definitions_path(str(root_dir / "eccodes/definitions")) def grib_def_ctx(grib_def: str): From 5f276a4ad3fcf1958656096678e0c7d4590e409e Mon Sep 17 00:00:00 2001 From: Carlos Osuna Date: Tue, 4 Feb 2025 17:41:31 +0100 Subject: [PATCH 3/4] hack so that the hgrid is passed from the load function --- src/meteodatalab/grib_decoder.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/meteodatalab/grib_decoder.py b/src/meteodatalab/grib_decoder.py index 401cd2b..aed6c95 100644 --- a/src/meteodatalab/grib_decoder.py +++ b/src/meteodatalab/grib_decoder.py @@ -83,6 +83,7 @@ def _is_ensemble(field) -> bool: def _get_hcoords(field): if field.metadata("gridType") == "unstructured_grid": + return {} grid_uuid = field.metadata("uuidOfHGrid") return icon_grid.get_icon_grid(grid_uuid) return { @@ -124,7 +125,7 @@ def __post_init__(self, is_ensemble: bool): if not is_ensemble: self.dims = self.dims[1:] - def load(self, field: GribField) -> None: + def load(self, field: GribField, hcoords=None) -> None: key = _get_key(field, self.dims) name = field.metadata(NAME_KEY) logger.debug("Received field for param: %s, key: %s", name, key) @@ -142,7 +143,7 @@ def load(self, field: GribField) -> None: } if not self.hcoords: - self.hcoords = _get_hcoords(field) + self.hcoords = hcoords if hcoords else _get_hcoords(field) def _gather_coords(self): coord_values = zip(*self.values) @@ -170,8 +171,7 @@ def to_xarray(self) -> xr.DataArray: ref_time = xr.DataArray(coords["ref_time"], dims="ref_time") lead_time = xr.DataArray(coords["lead_time"], dims="lead_time") tcoords = {"valid_time": ref_time + lead_time} - hdims = self.hcoords["lon"].dims - + hdims = self.hcoords["lon"].dims if "lon" in self.hcoords else ("cell",) array = xr.DataArray( data=np.array( [self.values.pop(key) for key in sorted(self.values)] @@ -188,8 +188,7 @@ def to_xarray(self) -> xr.DataArray: def _load_buffer_map( - source: data_source.DataSource, - request: Request, + source: data_source.DataSource, request: Request, hcoords=None ) -> dict[str, _FieldBuffer]: logger.info("Retrieving request: %s", request) fs = source.retrieve(request) @@ -202,7 +201,7 @@ def _load_buffer_map( buffer = buffer_map[name] else: buffer = buffer_map[name] = _FieldBuffer(_is_ensemble(field)) - buffer.load(field) + buffer.load(field, hcoords) return buffer_map @@ -246,8 +245,7 @@ def load_single_param( def load( - source: data_source.DataSource, - request: Request, + source: data_source.DataSource, request: Request, hcoords=None ) -> dict[str, xr.DataArray]: """Request data from a data source. @@ -269,7 +267,7 @@ def load( A mapping of shortName to data arrays of the requested fields. """ - buffer_map = _load_buffer_map(source, request) + buffer_map = _load_buffer_map(source, request, hcoords) result = {} for name, buffer in buffer_map.items(): try: From c7d4a58e39b409745da8b7a1ca3389408d0befeb Mon Sep 17 00:00:00 2001 From: Carlos Osuna Date: Wed, 5 Feb 2025 19:05:34 +0100 Subject: [PATCH 4/4] remove protection --- src/meteodatalab/operators/wind.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/meteodatalab/operators/wind.py b/src/meteodatalab/operators/wind.py index b84b1f9..4218d84 100644 --- a/src/meteodatalab/operators/wind.py +++ b/src/meteodatalab/operators/wind.py @@ -35,9 +35,6 @@ def speed(u: xr.DataArray, v: xr.DataArray) -> xr.DataArray: the horizontal wind speed [m/s]. """ - if u.origin_x != 0.0 or v.origin_y != 0.0: - raise ValueError("The wind components should not be staggered.") - name = {"U": "SP", "U_10M": "SP_10M"}[u.parameter["shortName"]] return xr.DataArray( np.sqrt(u**2 + v**2),