From 8b3d6d2e6f3b3ad701d7df58c76f4a5dea41b56e Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 15 Nov 2024 16:28:27 +0100 Subject: [PATCH 01/19] Add from_imod5_data for meteo_mappings --- imod/msw/meteo_mapping.py | 61 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 0fa3c818c..776cdcea5 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,18 +1,49 @@ from copy import deepcopy +from pathlib import Path from typing import Any, Optional, TextIO import numpy as np import pandas as pd import xarray as xr +import imod from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.prepare import common -from imod.typing import GridDataArray, IntArray +from imod.typing import GridDataArray, GridDataDict, IntArray from imod.util.regrid_method_type import RegridMethodType +def find_in_file_list(filename: str, paths: list[str]) -> str: + for file in paths: + if filename == Path(file[0]).name.lower(): + return file[0] + raise ValueError(f"could not find {filename} in list of paths: {paths}") + + +def open_first_meteo_grid(mete_grid_path: str, column_nr: int) -> xr.DataArray: + """ + Find first meteo grid path in mete_grid.inp. Only read the first grid, so it + can be used to generate meteomappings. + """ + if column_nr not in [2, 3]: + raise ValueError("Column nr should be 2 or 3") + + f = open(mete_grid_path, "r") + lines = f.readlines() + meteo_filepath = Path(lines[0].split(",")[column_nr].replace('"', "")) + return imod.rasterio.open(meteo_filepath) + + +def open_first_meteo_grid_from_imod5_data( + imod5_data: dict[str, GridDataDict], column_nr: int +): + paths = imod5_data["extra"]["paths"] + metegrid_path = find_in_file_list("mete_grid.inp", paths) + return open_first_meteo_grid(metegrid_path, column_nr=column_nr) + + class MeteoMapping(MetaSwapPackage): """ This class provides common methods for creating mappings between @@ -125,6 +156,20 @@ def __init__( super().__init__() self.meteo = precipitation + @classmethod + def from_imod5_data( + cls, imod5_data: dict[str, GridDataDict] + ) -> "PrecipitationMapping": + """ + Construct precipitation mapping from imod5 data. Opens first ascii grid + in mete_grid.inp, which is used to construct mappings to svats. The + grids should not change in dimension over time. No checks are done + whether cells switch from inactive to active or vice versa. + """ + column_nr = 2 + meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) + return cls(meteo_grid) + class EvapotranspirationMapping(MeteoMapping): """ @@ -156,3 +201,17 @@ def __init__( ): super().__init__() self.meteo = evapotranspiration + + @classmethod + def from_imod5_data( + cls, imod5_data: dict[str, GridDataDict] + ) -> "EvapotranspirationMapping": + """ + Construct evapotranspiration mapping from imod5 data. Opens first ascii + grid in mete_grid.inp, which is used to construct mappings to svats. The + grids should not change in dimension over time. No checks are done + whether cells switch from inactive to active or vice versa. + """ + column_nr = 3 + meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) + return cls(meteo_grid) From 19d9cc3ea517139dbb0b7a05e78a9127eb005f54 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 15 Nov 2024 17:30:01 +0100 Subject: [PATCH 02/19] Reduce code duplication in tests --- .../test_msw/test_precipitation_mapping.py | 239 ++++++------------ 1 file changed, 79 insertions(+), 160 deletions(-) diff --git a/imod/tests/test_msw/test_precipitation_mapping.py b/imod/tests/test_msw/test_precipitation_mapping.py index 3d74ac306..ccb28eb39 100644 --- a/imod/tests/test_msw/test_precipitation_mapping.py +++ b/imod/tests/test_msw/test_precipitation_mapping.py @@ -2,6 +2,7 @@ from pathlib import Path import numpy as np +import pandas as pd import pytest import xarray as xr from numpy.testing import assert_equal @@ -9,30 +10,8 @@ from imod import msw -def test_precipitation_mapping_simple(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on - +@pytest.fixture(scope="function") +def svat_index() -> xr.DataArray: x_svat = [1.0, 2.0, 3.0] y_svat = [1.0, 2.0, 3.0] subunit_svat = [0, 1] @@ -57,31 +36,12 @@ def test_precipitation_mapping_simple(fixed_format_parser): ) # fmt: on index = (svat != 0).values.ravel() + return svat, index - precipitation_mapping = msw.PrecipitationMapping(precipitation) - - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - precipitation_mapping.write(output_dir, index, svat, None, None) - - results = fixed_format_parser( - output_dir / msw.PrecipitationMapping._file_name, - msw.PrecipitationMapping._metadata_dict, - ) - assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) - - -def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [4.5, 2.5, 0.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = -2.0 +def create_meteo_grid(x, y, subunit, dx, dy): # fmt: off - precipitation = xr.DataArray( + return xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -94,35 +54,21 @@ def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} ) # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], +def test_precipitation_mapping_simple(fixed_format_parser, svat_index): + svat, index = svat_index - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() + x = [-0.5, 1.5, 3.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) with tempfile.TemporaryDirectory() as output_dir: @@ -135,58 +81,22 @@ def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): ) assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([3, 2, 3, 2, 2, 2])) + assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [4.5, 2.5, 0.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = -2.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [3.0, 2.0, 1.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 +def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser, svat_index): + svat, index = svat_index - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) + with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) precipitation_mapping.write(output_dir, index, svat, None, None) @@ -197,60 +107,46 @@ def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser): ) assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([2, 3, 2, 2, 2, 2])) + assert_equal(results["row"], np.array([3, 2, 3, 2, 2, 2])) assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_precipitation_mapping_out_of_bound(): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [2.5, 4.5, 6.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 +def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser, svat_index): + svat, index = svat_index + svat = svat.assign_coords(y=[3.0, 2.0, 1.0], dy=-1.0) - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on + precipitation = create_meteo_grid(x, y, subunit, dx, dy) + precipitation_mapping = msw.PrecipitationMapping(precipitation) + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + precipitation_mapping.write(output_dir, index, svat, None, None) - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 + results = fixed_format_parser( + output_dir / msw.PrecipitationMapping._file_name, + msw.PrecipitationMapping._metadata_dict, + ) - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) + assert_equal(results["row"], np.array([2, 3, 2, 2, 2, 2])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() +def test_precipitation_mapping_out_of_bound(svat_index): + svat, index = svat_index + + x = [-0.5, 1.5, 3.5] + y = [2.5, 4.5, 6.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) with tempfile.TemporaryDirectory() as output_dir: @@ -258,3 +154,26 @@ def test_precipitation_mapping_out_of_bound(): # The grid is out of bounds, which is why we expect a ValueError to be raisen with pytest.raises(ValueError): precipitation_mapping.write(output_dir, index, svat, None, None) + + +def test_from_imod5(tmpdir_factory): + datadir = tmpdir_factory.mktemp("precipitation_mapping") + + x = [-0.5, 1.5, 3.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + + time = [np.datetime64(t) for t in["2001-01-01", "2001-01-02", "2001-01-03"]] + time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) + + precipitation = create_meteo_grid(x, y, subunit, dx, dy) + precipitation_times = time_da * precipitation + mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) + mete_grid.write(datadir) + + imod5_data = {"paths": [["foo"], [datadir / "mete_grid.inp"], ["bar"]]} + precipitation_mapping = msw.PrecipitationMapping.from_imod5_data(imod5_data) + + xr.testing.assert_equal(precipitation, precipitation_mapping.meteo) From 700b349a5222d261ac1bd214e7d75dc37d5ba682 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 15 Nov 2024 17:31:50 +0100 Subject: [PATCH 03/19] Add from_imod5_data test and search with proper meteo grid path --- imod/msw/meteo_mapping.py | 2 +- .../test_msw/test_precipitation_mapping.py | 22 ++++++++++++------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 776cdcea5..5bf5ca049 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -33,7 +33,7 @@ def open_first_meteo_grid(mete_grid_path: str, column_nr: int) -> xr.DataArray: f = open(mete_grid_path, "r") lines = f.readlines() meteo_filepath = Path(lines[0].split(",")[column_nr].replace('"', "")) - return imod.rasterio.open(meteo_filepath) + return imod.rasterio.open(mete_grid_path / ".." / meteo_filepath) def open_first_meteo_grid_from_imod5_data( diff --git a/imod/tests/test_msw/test_precipitation_mapping.py b/imod/tests/test_msw/test_precipitation_mapping.py index ccb28eb39..d5138aa5c 100644 --- a/imod/tests/test_msw/test_precipitation_mapping.py +++ b/imod/tests/test_msw/test_precipitation_mapping.py @@ -2,7 +2,6 @@ from pathlib import Path import numpy as np -import pandas as pd import pytest import xarray as xr from numpy.testing import assert_equal @@ -157,23 +156,30 @@ def test_precipitation_mapping_out_of_bound(svat_index): def test_from_imod5(tmpdir_factory): - datadir = tmpdir_factory.mktemp("precipitation_mapping") + datadir = Path(tmpdir_factory.mktemp("precipitation_mapping")) + # Arrange x = [-0.5, 1.5, 3.5] - y = [0.5, 2.5, 4.5] + y = [4.5, 2.5, 0.5] subunit = [0, 1] dx = 2.0 - dy = 2.0 + dy = -2.0 - time = [np.datetime64(t) for t in["2001-01-01", "2001-01-02", "2001-01-03"]] + time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) - precipitation = create_meteo_grid(x, y, subunit, dx, dy) + precipitation = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) precipitation_times = time_da * precipitation mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) mete_grid.write(datadir) - imod5_data = {"paths": [["foo"], [datadir / "mete_grid.inp"], ["bar"]]} + paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] + imod5_data = {"extra": {"paths": paths}} + + # Act precipitation_mapping = msw.PrecipitationMapping.from_imod5_data(imod5_data) + actual = precipitation_mapping.meteo - xr.testing.assert_equal(precipitation, precipitation_mapping.meteo) + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(precipitation, actual.isel(time=0, drop=True)) From 4fbe5c81ac9631e3a0077510b1c96e85dfc6afc5 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 18 Nov 2024 16:06:34 +0100 Subject: [PATCH 04/19] Better test name --- imod/tests/test_msw/test_precipitation_mapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/tests/test_msw/test_precipitation_mapping.py b/imod/tests/test_msw/test_precipitation_mapping.py index d5138aa5c..cc537541e 100644 --- a/imod/tests/test_msw/test_precipitation_mapping.py +++ b/imod/tests/test_msw/test_precipitation_mapping.py @@ -155,7 +155,7 @@ def test_precipitation_mapping_out_of_bound(svat_index): precipitation_mapping.write(output_dir, index, svat, None, None) -def test_from_imod5(tmpdir_factory): +def test_precipitation_from_imod5(tmpdir_factory): datadir = Path(tmpdir_factory.mktemp("precipitation_mapping")) # Arrange From 0ea21ed7c8d7c15a29f255e758574d690c25eb69 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 18 Nov 2024 16:07:04 +0100 Subject: [PATCH 05/19] Clean up tests and add from_imod5 test --- .../test_evapotranspiration_mapping.py | 191 +++++++----------- 1 file changed, 76 insertions(+), 115 deletions(-) diff --git a/imod/tests/test_msw/test_evapotranspiration_mapping.py b/imod/tests/test_msw/test_evapotranspiration_mapping.py index 0d5278f58..e969fe515 100644 --- a/imod/tests/test_msw/test_evapotranspiration_mapping.py +++ b/imod/tests/test_msw/test_evapotranspiration_mapping.py @@ -9,29 +9,8 @@ from imod import msw -def test_evapotranspiration_mapping_simple(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 - # fmt: off - evapotranspiration = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on +@pytest.fixture(scope="function") +def svat_index() -> xr.DataArray: x_svat = [1.0, 2.0, 3.0] y_svat = [1.0, 2.0, 3.0] subunit_svat = [0, 1] @@ -56,31 +35,11 @@ def test_evapotranspiration_mapping_simple(fixed_format_parser): ) # fmt: on index = (svat != 0).values.ravel() + return svat, index - evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) - - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - evapotranspiration_mapping.write(output_dir, index, svat, None, None) - - results = fixed_format_parser( - output_dir / msw.EvapotranspirationMapping._file_name, - msw.EvapotranspirationMapping._metadata_dict, - ) - - assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) - - -def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): - x_meteo = [3.5, 1.5, -0.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = -2.0 - dy_meteo = 2.0 +def create_meteo_grid(x, y, subunit, dx, dy): # fmt: off - evapotranspiration = xr.DataArray( + return xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -93,34 +52,19 @@ def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} ) # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() +def test_evapotranspiration_mapping_simple(fixed_format_parser, svat_index): + svat, index = svat_index + x = [-0.5, 1.5, 3.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) with tempfile.TemporaryDirectory() as output_dir: @@ -134,58 +78,45 @@ def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 1])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_evapotranspiration_mapping_out_of_bound(): - x_meteo = [-2.5, -0.5, 1.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 +def test_evapotranspiration_mapping_negative_dx(fixed_format_parser, svat_index): + svat, index = svat_index - # fmt: off - evapotranspiration = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], + x = [3.5, 1.5, -0.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = -2.0 + dy = 2.0 - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) + evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + evapotranspiration_mapping.write(output_dir, index, svat, None, None) + + results = fixed_format_parser( + output_dir / msw.EvapotranspirationMapping._file_name, + msw.EvapotranspirationMapping._metadata_dict, + ) + + assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) + assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 1])) - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() +def test_evapotranspiration_mapping_out_of_bound(svat_index): + svat, index = svat_index + + x = [-2.5, -0.5, 1.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) with tempfile.TemporaryDirectory() as output_dir: @@ -193,3 +124,33 @@ def test_evapotranspiration_mapping_out_of_bound(): # The grid is out of bounds, which is why we expect a ValueError to be raisen with pytest.raises(ValueError): evapotranspiration_mapping.write(output_dir, index, svat, None, None) + + +def test_evapotranspiration_from_imod5(tmpdir_factory): + datadir = Path(tmpdir_factory.mktemp("evapotranspiration_mapping")) + + # Arrange + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 + + time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] + time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) + + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) + evapotranspiration_times = time_da * evapotranspiration + mete_grid = msw.MeteoGrid(evapotranspiration_times, evapotranspiration_times) + mete_grid.write(datadir) + + paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] + imod5_data = {"extra": {"paths": paths}} + + # Act + evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data(imod5_data) + actual = evapotranspiration_mapping.meteo + + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(evapotranspiration, actual.isel(time=0, drop=True)) From f40adc414f7a4e0b35a3280c6c8324c3f7f31212 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 18 Nov 2024 16:08:05 +0100 Subject: [PATCH 06/19] Format --- imod/tests/test_msw/test_evapotranspiration_mapping.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_msw/test_evapotranspiration_mapping.py b/imod/tests/test_msw/test_evapotranspiration_mapping.py index e969fe515..798d8877b 100644 --- a/imod/tests/test_msw/test_evapotranspiration_mapping.py +++ b/imod/tests/test_msw/test_evapotranspiration_mapping.py @@ -37,6 +37,7 @@ def svat_index() -> xr.DataArray: index = (svat != 0).values.ravel() return svat, index + def create_meteo_grid(x, y, subunit, dx, dy): # fmt: off return xr.DataArray( @@ -56,6 +57,7 @@ def create_meteo_grid(x, y, subunit, dx, dy): ) # fmt: on + def test_evapotranspiration_mapping_simple(fixed_format_parser, svat_index): svat, index = svat_index @@ -139,7 +141,9 @@ def test_evapotranspiration_from_imod5(tmpdir_factory): time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) - evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy).isel( + subunit=0, drop=True + ) evapotranspiration_times = time_da * evapotranspiration mete_grid = msw.MeteoGrid(evapotranspiration_times, evapotranspiration_times) mete_grid.write(datadir) @@ -148,7 +152,9 @@ def test_evapotranspiration_from_imod5(tmpdir_factory): imod5_data = {"extra": {"paths": paths}} # Act - evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data(imod5_data) + evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data( + imod5_data + ) actual = evapotranspiration_mapping.meteo # Assert From 6c62638ea26e101a71a57a20db0f4c6ac139ea42 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 18 Nov 2024 17:01:18 +0100 Subject: [PATCH 07/19] move file_in_file_list to common utilities --- imod/msw/meteo_mapping.py | 8 +------- imod/msw/utilities/common.py | 9 +++++++++ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 5bf5ca049..6e7a376f6 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -10,18 +10,12 @@ from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage +from imod.msw.utilities.common import find_in_file_list from imod.prepare import common from imod.typing import GridDataArray, GridDataDict, IntArray from imod.util.regrid_method_type import RegridMethodType -def find_in_file_list(filename: str, paths: list[str]) -> str: - for file in paths: - if filename == Path(file[0]).name.lower(): - return file[0] - raise ValueError(f"could not find {filename} in list of paths: {paths}") - - def open_first_meteo_grid(mete_grid_path: str, column_nr: int) -> xr.DataArray: """ Find first meteo grid path in mete_grid.inp. Only read the first grid, so it diff --git a/imod/msw/utilities/common.py b/imod/msw/utilities/common.py index 836902662..db69aab33 100644 --- a/imod/msw/utilities/common.py +++ b/imod/msw/utilities/common.py @@ -1,6 +1,15 @@ +from pathlib import Path + from imod.typing import GridDataArray from imod.typing.grid import concat def concat_imod5(arg1: GridDataArray, arg2: GridDataArray) -> GridDataArray: return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1]) + + +def find_in_file_list(filename: str, paths: list[str]) -> str: + for file in paths: + if filename == Path(file[0]).name.lower(): + return file[0] + raise ValueError(f"could not find {filename} in list of paths: {paths}") From 53770b5eb6ca91ded8abbe4e299e2ec4e033a662 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Tue, 19 Nov 2024 09:37:56 +0100 Subject: [PATCH 08/19] Start working on MeteoGridCopy class --- imod/msw/__init__.py | 2 +- imod/msw/meteo_grid.py | 38 ++++++++++++++++++++++++++ imod/tests/test_msw/test_meteo_grid.py | 15 +++++++++- 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/imod/msw/__init__.py b/imod/msw/__init__.py index bcbe174aa..83f8df90b 100644 --- a/imod/msw/__init__.py +++ b/imod/msw/__init__.py @@ -9,7 +9,7 @@ InitialConditionsSavedState, ) from imod.msw.landuse import LanduseOptions -from imod.msw.meteo_grid import MeteoGrid +from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy from imod.msw.meteo_mapping import EvapotranspirationMapping, PrecipitationMapping from imod.msw.model import MetaSwapModel from imod.msw.output_control import TimeOutputControl, VariableOutputControl diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index d8c283252..7d166cb96 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,5 +1,6 @@ import csv from pathlib import Path +from shutil import copyfile from typing import Optional, Union import numpy as np @@ -11,6 +12,9 @@ from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod from imod.msw.timeutil import to_metaswap_timeformat +from imod.msw.utilities.common import find_in_file_list +from imod.typing import GridDataDict +from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType class MeteoGrid(MetaSwapPackage, IRegridPackage): @@ -188,3 +192,37 @@ def _pkgcheck(self): f"Received excess dims {excess_dims} in {self.__class__} for " f"{varname}, please provide data with {allowed_dims}" ) + + +class MeteoGridCopy(MetaSwapPackage, IRegridPackage): + """ + Class to copy existing ``mete_grid.inp``, which contains the meteorological + grid data. Next to a MeteoGridCopy instance, instances of + PrecipitationMapping and EvapotranspirationMapping are required as well to + specify meteorological information to MetaSWAP. + + Parameters + ---------- + path: Path to mete_grid.inp file + """ + + _file_name = "mete_grid.inp" + _meteo_dirname = "meteo_grids" + + _regrid_method: RegridMethodType = EmptyRegridMethod() + + def __init__(self, path: Path | str): + super().__init__() + self.dataset["path"] = path + + def write(self, directory: Path | str, *args): + path_metegrid = Path(self.dataset["path"]) + new_path = directory / self._file_name + copyfile(path_metegrid, new_path) + + @classmethod + def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "MeteoGridCopy": + paths = imod5_data["extra"]["paths"] + filepath = find_in_file_list(cls._file_name, paths) + + return cls(filepath) diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index 8cca2cc06..879afd0ad 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -11,7 +11,7 @@ from numpy.testing import assert_equal from imod.mf6.utilities.regrid import RegridderWeightsCache -from imod.msw import MeteoGrid +from imod.msw import MeteoGrid, MeteoGridCopy def setup_meteo_grid(): @@ -123,3 +123,16 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["x"].values == new_grid["x"].values) assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) + + +def test_meteocopy(): + # Arrange + meteo_grid = setup_meteo_grid() + + with tempfile.TemporaryDirectory() as output_dir: + grid_dir = Path(output_dir) / "grid" + meteo_grid.write(grid_dir) + + meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") + + From c0be8edae61fa955d3e31fe45b414ef57f34c19b Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 14:58:09 +0100 Subject: [PATCH 09/19] Add test for meteo grid copy --- imod/tests/test_msw/test_meteo_grid.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index 879afd0ad..00a392d01 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -1,4 +1,5 @@ import csv +import filecmp import os import tempfile from pathlib import Path @@ -125,7 +126,7 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) -def test_meteocopy(): +def test_meteo_grid_copy(): # Arrange meteo_grid = setup_meteo_grid() @@ -134,5 +135,7 @@ def test_meteocopy(): meteo_grid.write(grid_dir) meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") - + copy_dir = Path(output_dir) / "copied" + meteo_grid_copy.write(copy_dir) + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") From 7ad968f47270ed72dc4498d991c182864b1600f3 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 15:23:15 +0100 Subject: [PATCH 10/19] fix bug where path was not properly unpacked from dataset --- imod/msw/meteo_grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 7d166cb96..2a44972b2 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -216,7 +216,7 @@ def __init__(self, path: Path | str): self.dataset["path"] = path def write(self, directory: Path | str, *args): - path_metegrid = Path(self.dataset["path"]) + path_metegrid = Path(self.dataset["path"].values[()]) new_path = directory / self._file_name copyfile(path_metegrid, new_path) From 7acf0337e085249c0eb55ef37e3ea820747f108e Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 15:23:38 +0100 Subject: [PATCH 11/19] Add test for from_imod5_data --- imod/tests/test_msw/test_meteo_grid.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index 00a392d01..927a5097f 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -126,16 +126,40 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) -def test_meteo_grid_copy(): +def test_meteogridcopy_write(): # Arrange meteo_grid = setup_meteo_grid() with tempfile.TemporaryDirectory() as output_dir: grid_dir = Path(output_dir) / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) meteo_grid.write(grid_dir) meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") copy_dir = Path(output_dir) / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") + + +def test_meteogridcopy_from_imod5(): + meteo_grid = setup_meteo_grid() + with tempfile.TemporaryDirectory() as output_dir: + grid_dir = Path(output_dir) / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + + imod5_data = {} + imod5_data["extra"] = {} + imod5_data["extra"]["paths"] = [["foo"], [(grid_dir / "mete_grid.inp").resolve()], ["bar"]] + + meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) + copy_dir = Path(output_dir) / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") From cc2b71d6733173f2a4456fa0c93a328958170afb Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 15:24:20 +0100 Subject: [PATCH 12/19] format --- imod/tests/test_msw/test_meteo_grid.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index 927a5097f..b5108842e 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -154,7 +154,11 @@ def test_meteogridcopy_from_imod5(): imod5_data = {} imod5_data["extra"] = {} - imod5_data["extra"]["paths"] = [["foo"], [(grid_dir / "mete_grid.inp").resolve()], ["bar"]] + imod5_data["extra"]["paths"] = [ + ["foo"], + [(grid_dir / "mete_grid.inp").resolve()], + ["bar"], + ] meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) copy_dir = Path(output_dir) / "copied" From 316b4242c5258d5203c4f891a06d718ea51bf761 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 15:45:47 +0100 Subject: [PATCH 13/19] Update changelog --- docs/api/changelog.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 8c7aac2ed..d0ecceaef 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -9,6 +9,12 @@ The format is based on `Keep a Changelog`_, and this project adheres to [Unreleased] ------------ +Added +~~~~~ + +- :class:`imod.msw.MeteoGridCopy` to copy existing `mete_grid.inp` files, so + ASCII grids in large existing meteo databases do not have to be read. + Changed ~~~~~~~ From 9566f98579463f41057e281221fad6f33eadbdc3 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 22 Nov 2024 16:15:08 +0100 Subject: [PATCH 14/19] Fix mypy errors --- imod/msw/meteo_grid.py | 11 ++++++----- imod/msw/meteo_mapping.py | 22 +++++++++------------- imod/typing/__init__.py | 1 + 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 2a44972b2..70d3625a5 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,7 +1,7 @@ import csv from pathlib import Path from shutil import copyfile -from typing import Optional, Union +from typing import Optional, Union, cast import numpy as np import pandas as pd @@ -13,7 +13,7 @@ from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod from imod.msw.timeutil import to_metaswap_timeformat from imod.msw.utilities.common import find_in_file_list -from imod.typing import GridDataDict +from imod.typing import Imod5DataDict from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType @@ -216,13 +216,14 @@ def __init__(self, path: Path | str): self.dataset["path"] = path def write(self, directory: Path | str, *args): - path_metegrid = Path(self.dataset["path"].values[()]) + directory = Path(directory) + path_metegrid = Path(str(self.dataset["path"].values[()])) new_path = directory / self._file_name copyfile(path_metegrid, new_path) @classmethod - def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "MeteoGridCopy": - paths = imod5_data["extra"]["paths"] + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": + paths = cast(list[str], imod5_data["extra"]["paths"]) filepath = find_in_file_list(cls._file_name, paths) return cls(filepath) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 6e7a376f6..6b78f53b6 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,6 +1,6 @@ from copy import deepcopy from pathlib import Path -from typing import Any, Optional, TextIO +from typing import Any, Optional, TextIO, cast import numpy as np import pandas as pd @@ -12,11 +12,11 @@ from imod.msw.pkgbase import MetaSwapPackage from imod.msw.utilities.common import find_in_file_list from imod.prepare import common -from imod.typing import GridDataArray, GridDataDict, IntArray +from imod.typing import GridDataArray, Imod5DataDict, IntArray from imod.util.regrid_method_type import RegridMethodType -def open_first_meteo_grid(mete_grid_path: str, column_nr: int) -> xr.DataArray: +def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.DataArray: """ Find first meteo grid path in mete_grid.inp. Only read the first grid, so it can be used to generate meteomappings. @@ -24,16 +24,16 @@ def open_first_meteo_grid(mete_grid_path: str, column_nr: int) -> xr.DataArray: if column_nr not in [2, 3]: raise ValueError("Column nr should be 2 or 3") + mete_grid_path = Path(mete_grid_path) + f = open(mete_grid_path, "r") lines = f.readlines() meteo_filepath = Path(lines[0].split(",")[column_nr].replace('"', "")) return imod.rasterio.open(mete_grid_path / ".." / meteo_filepath) -def open_first_meteo_grid_from_imod5_data( - imod5_data: dict[str, GridDataDict], column_nr: int -): - paths = imod5_data["extra"]["paths"] +def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int): + paths = cast(list[str], imod5_data["extra"]["paths"]) metegrid_path = find_in_file_list("mete_grid.inp", paths) return open_first_meteo_grid(metegrid_path, column_nr=column_nr) @@ -151,9 +151,7 @@ def __init__( self.meteo = precipitation @classmethod - def from_imod5_data( - cls, imod5_data: dict[str, GridDataDict] - ) -> "PrecipitationMapping": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "PrecipitationMapping": """ Construct precipitation mapping from imod5 data. Opens first ascii grid in mete_grid.inp, which is used to construct mappings to svats. The @@ -197,9 +195,7 @@ def __init__( self.meteo = evapotranspiration @classmethod - def from_imod5_data( - cls, imod5_data: dict[str, GridDataDict] - ) -> "EvapotranspirationMapping": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMapping": """ Construct evapotranspiration mapping from imod5 data. Opens first ascii grid in mete_grid.inp, which is used to construct mappings to svats. The diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 760346894..8beb7fd13 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -12,6 +12,7 @@ GridDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] GridDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] GridDataDict: TypeAlias = dict[str, GridDataArray] +Imod5DataDict: TypeAlias = dict[str, GridDataDict | dict[str, list[str]]] ScalarAsDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] ScalarAsDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray] From 961b07a4a101b427dfa1c0255d5d9cb2a08d2acc Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 25 Nov 2024 14:04:03 +0100 Subject: [PATCH 15/19] Open mete_grid.inp in context and only read first line --- imod/msw/meteo_mapping.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 6b78f53b6..068214c63 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -18,17 +18,17 @@ def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.DataArray: """ - Find first meteo grid path in mete_grid.inp. Only read the first grid, so it - can be used to generate meteomappings. + Find and open first meteo grid path in mete_grid.inp. This grid is enough to + generate meteomappings. """ if column_nr not in [2, 3]: raise ValueError("Column nr should be 2 or 3") mete_grid_path = Path(mete_grid_path) - f = open(mete_grid_path, "r") - lines = f.readlines() - meteo_filepath = Path(lines[0].split(",")[column_nr].replace('"', "")) + with open(mete_grid_path, "r") as f: + first_line = f.readline() + meteo_filepath = Path(first_line.split(",")[column_nr].replace('"', "")) return imod.rasterio.open(mete_grid_path / ".." / meteo_filepath) From e823ffbc2bc6ba663b3518f6ef6a74548b8bcd4c Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 26 Nov 2024 16:33:16 +0100 Subject: [PATCH 16/19] Add test cases where mete_grid.inp is adapted, so that paths are replaced with floats. --- .../test_evapotranspiration_mapping.py | 63 ++++++++++++++----- .../test_msw/test_precipitation_mapping.py | 51 ++++++++++++--- 2 files changed, 91 insertions(+), 23 deletions(-) diff --git a/imod/tests/test_msw/test_evapotranspiration_mapping.py b/imod/tests/test_msw/test_evapotranspiration_mapping.py index 798d8877b..b2752a60f 100644 --- a/imod/tests/test_msw/test_evapotranspiration_mapping.py +++ b/imod/tests/test_msw/test_evapotranspiration_mapping.py @@ -1,3 +1,4 @@ +import re import tempfile from pathlib import Path @@ -5,6 +6,7 @@ import pytest import xarray as xr from numpy.testing import assert_equal +from pytest_cases import parametrize_with_cases from imod import msw @@ -128,9 +130,8 @@ def test_evapotranspiration_mapping_out_of_bound(svat_index): evapotranspiration_mapping.write(output_dir, index, svat, None, None) -def test_evapotranspiration_from_imod5(tmpdir_factory): - datadir = Path(tmpdir_factory.mktemp("evapotranspiration_mapping")) - +def setup_meteo_grid(datadir): + """Setup precipitation grid and write mete_grid.inp""" # Arrange x = [-0.5, 1.5, 3.5] y = [4.5, 2.5, 0.5] @@ -141,22 +142,54 @@ def test_evapotranspiration_from_imod5(tmpdir_factory): time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) - evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy).isel( - subunit=0, drop=True - ) - evapotranspiration_times = time_da * evapotranspiration - mete_grid = msw.MeteoGrid(evapotranspiration_times, evapotranspiration_times) + precipitation = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) + precipitation_times = time_da * precipitation + mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) mete_grid.write(datadir) + return precipitation + + +class MeteGridCases: + """ + Cases return strings to replace in mete_grid.inp, to set paths to floats. + """ + + def case_all_paths(self): + return r"nothing_to_replace" + + def case_some_paths(self): + return r'"meteo_grids\\(\w+)_20010101000000.asc"' + + def case_no_paths(self): + return r'"meteo_grids\\(\w+)_([0-9]+).asc"' + + +@parametrize_with_cases("replace_string", cases=MeteGridCases) +def test_evapotranspiration_from_imod5(tmpdir_factory, replace_string, request): + datadir = Path(tmpdir_factory.mktemp("evapotranspiration_mapping")) + evapotranspiration = setup_meteo_grid(datadir) paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] imod5_data = {"extra": {"paths": paths}} + # Replace text in existing mete_grid.inp + with open(datadir / "mete_grid.inp", "r") as f: + text = f.read() + text_replaced = re.sub(replace_string, '"0.0"', text) + with open(datadir / "mete_grid.inp", "w") as f: + f.write(text_replaced) + # Act - evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data( - imod5_data - ) - actual = evapotranspiration_mapping.meteo + if request.node.callspec.id == "no_paths": + with pytest.raises(ValueError): + msw.EvapotranspirationMapping.from_imod5_data(imod5_data) + + else: + evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data( + imod5_data + ) + actual = evapotranspiration_mapping.meteo - # Assert - assert len(actual.coords["time"]) == 1 - xr.testing.assert_equal(evapotranspiration, actual.isel(time=0, drop=True)) + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(evapotranspiration, actual.isel(time=0, drop=True)) diff --git a/imod/tests/test_msw/test_precipitation_mapping.py b/imod/tests/test_msw/test_precipitation_mapping.py index cc537541e..e9162fdf0 100644 --- a/imod/tests/test_msw/test_precipitation_mapping.py +++ b/imod/tests/test_msw/test_precipitation_mapping.py @@ -1,3 +1,4 @@ +import re import tempfile from pathlib import Path @@ -5,6 +6,7 @@ import pytest import xarray as xr from numpy.testing import assert_equal +from pytest_cases import parametrize_with_cases from imod import msw @@ -155,9 +157,8 @@ def test_precipitation_mapping_out_of_bound(svat_index): precipitation_mapping.write(output_dir, index, svat, None, None) -def test_precipitation_from_imod5(tmpdir_factory): - datadir = Path(tmpdir_factory.mktemp("precipitation_mapping")) - +def setup_meteo_grid(datadir): + """Setup precipitation grid and write mete_grid.inp""" # Arrange x = [-0.5, 1.5, 3.5] y = [4.5, 2.5, 0.5] @@ -172,14 +173,48 @@ def test_precipitation_from_imod5(tmpdir_factory): precipitation_times = time_da * precipitation mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) mete_grid.write(datadir) + return precipitation + + +class MeteGridCases: + """ + Cases return strings to replace in mete_grid.inp, to set paths to floats. + """ + + def case_all_paths(self): + return r"nothing_to_replace" + + def case_some_paths(self): + return r'"meteo_grids\\(\w+)_20010101000000.asc"' + + def case_no_paths(self): + return r'"meteo_grids\\(\w+)_([0-9]+).asc"' + +@parametrize_with_cases("replace_string", cases=MeteGridCases) +def test_precipitation_from_imod5(tmpdir_factory, replace_string, request): + # Arrange + datadir = Path(tmpdir_factory.mktemp("precipitation_mapping")) + precipitation = setup_meteo_grid(datadir) paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] imod5_data = {"extra": {"paths": paths}} + # Replace text in existing mete_grid.inp + with open(datadir / "mete_grid.inp", "r") as f: + text = f.read() + text_replaced = re.sub(replace_string, '"0.0"', text) + with open(datadir / "mete_grid.inp", "w") as f: + f.write(text_replaced) + # Act - precipitation_mapping = msw.PrecipitationMapping.from_imod5_data(imod5_data) - actual = precipitation_mapping.meteo + if request.node.callspec.id == "no_paths": + with pytest.raises(ValueError): + msw.PrecipitationMapping.from_imod5_data(imod5_data) + + else: + precipitation_mapping = msw.PrecipitationMapping.from_imod5_data(imod5_data) + actual = precipitation_mapping.meteo - # Assert - assert len(actual.coords["time"]) == 1 - xr.testing.assert_equal(precipitation, actual.isel(time=0, drop=True)) + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(precipitation, actual.isel(time=0, drop=True)) From a9d9b43bdef5b8510cfb436658a106327f76d697 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 26 Nov 2024 16:34:17 +0100 Subject: [PATCH 17/19] Look for first ascii file, skip floats --- imod/msw/meteo_mapping.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 068214c63..3743e832e 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,5 +1,6 @@ from copy import deepcopy from pathlib import Path +from textwrap import dedent from typing import Any, Optional, TextIO, cast import numpy as np @@ -16,6 +17,20 @@ from imod.util.regrid_method_type import RegridMethodType +def _is_parsable_and_existing_path(potential_path: str, mete_grid_path: Path): + """ + mete_grid.inp can contain values like "0.", which are converted to float by + MetaSWAP. String is converted to path and checked if existing path. + """ + try: + float(potential_path) + return False + except ValueError: + # Resolve paths relative to mete_grid.inp path. + path = mete_grid_path / ".." / Path(potential_path) + return path.is_file() + + def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.DataArray: """ Find and open first meteo grid path in mete_grid.inp. This grid is enough to @@ -25,11 +40,20 @@ def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.Data raise ValueError("Column nr should be 2 or 3") mete_grid_path = Path(mete_grid_path) - with open(mete_grid_path, "r") as f: - first_line = f.readline() - meteo_filepath = Path(first_line.split(",")[column_nr].replace('"', "")) - return imod.rasterio.open(mete_grid_path / ".." / meteo_filepath) + lines = f.readlines() + + potential_paths = [line.split(",")[column_nr].replace('"', "") for line in lines] + for potential_path in potential_paths: + if _is_parsable_and_existing_path(potential_path, mete_grid_path): + resolved_path = mete_grid_path / ".." / Path(potential_path) + return imod.rasterio.open(resolved_path) + + error_message = dedent(f""" + Did not find parsable path to existing .ASC file in column {column_nr}. Got + values (printing first 10): {potential_paths[:10]}.""") + + raise ValueError(error_message) def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int): From 29e584eba7daf2e989eb7fabfe65db2f7abb6f11 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 26 Nov 2024 16:35:46 +0100 Subject: [PATCH 18/19] type annotate --- imod/msw/meteo_mapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 3743e832e..6877ad4a5 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -17,7 +17,7 @@ from imod.util.regrid_method_type import RegridMethodType -def _is_parsable_and_existing_path(potential_path: str, mete_grid_path: Path): +def _is_parsable_and_existing_path(potential_path: str, mete_grid_path: Path) -> bool: """ mete_grid.inp can contain values like "0.", which are converted to float by MetaSWAP. String is converted to path and checked if existing path. From 47a1e7ba31fec141f34459a5dda2aa7e13decc07 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 26 Nov 2024 16:43:46 +0100 Subject: [PATCH 19/19] Extend docstring --- imod/msw/meteo_mapping.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 6877ad4a5..1de00fa69 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -34,7 +34,8 @@ def _is_parsable_and_existing_path(potential_path: str, mete_grid_path: Path) -> def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.DataArray: """ Find and open first meteo grid path in mete_grid.inp. This grid is enough to - generate meteomappings. + generate meteomappings. There can be floats before in the column which + should be skipped. """ if column_nr not in [2, 3]: raise ValueError("Column nr should be 2 or 3")