From d11b76d1d3ed5240f151c859c11b903f14ee09ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Thu, 24 Aug 2023 14:22:47 +0200 Subject: [PATCH 1/3] Stop overwriting kh and gv of geotop in add_geotop_to_regis_layers --- nlmod/dims/layers.py | 2 +- nlmod/read/geotop.py | 2 +- nlmod/read/regis.py | 29 +++++++++++++++++++++++------ 3 files changed, 25 insertions(+), 8 deletions(-) diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py index c400abff..62f2a1bd 100644 --- a/nlmod/dims/layers.py +++ b/nlmod/dims/layers.py @@ -769,7 +769,7 @@ def get_kh_kv(kh, kv, anisotropy, fill_value_kh=1.0, fill_value_kv=0.1, idomain= logger.info(f"kv and kh both undefined in layer {layer}") # fill kh by kv * anisotropy - msg_suffix = f" of kh by multipying kv by an anisotropy of {anisotropy}" + msg_suffix = f" of kh by multipying kv with an anisotropy of {anisotropy}" kh = _fill_var(kh, kv * anisotropy, idomain, msg_suffix) # fill kv by kh / anisotropy diff --git a/nlmod/read/geotop.py b/nlmod/read/geotop.py index 3b9f183b..568cb2a9 100644 --- a/nlmod/read/geotop.py +++ b/nlmod/read/geotop.py @@ -360,7 +360,7 @@ def add_kh_and_kv( kh : str, optional THe name of the new variable with kh values in gt. The default is "kh". kv : str, optional - THe name of the new variable with kv values in gt. The default is "kv". + The name of the new variable with kv values in gt. The default is "kv". kh_df : str, optional The name of the column with kh values in df. The default is "kh". kv_df : str, optional diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py index 80b88b75..74b51a3f 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -185,7 +185,7 @@ def get_regis( def add_geotop_to_regis_layers( - rg, gt, layers="HLc", geotop_k=None, remove_nan_layers=True + rg, gt, layers="HLc", geotop_k=None, remove_nan_layers=True, anisotropy=1.0 ): """Combine geotop and regis in such a way that the one or more layers in Regis are replaced by the geo_eenheden of geotop. @@ -203,6 +203,9 @@ def add_geotop_to_regis_layers( DataFrame must at least contain columns 'lithok' and 'kh'. remove_nan_layers : bool, optional When True, layers with only 0 or NaN thickness are removed. The default is True. + anisotropy : float, optional + The anisotropy value used when there are no kv values in df. The default is + 1.0. Returns ------- @@ -211,8 +214,25 @@ def add_geotop_to_regis_layers( """ if isinstance(layers, str): layers = [layers] - if geotop_k is None: - geotop_k = geotop.get_lithok_props() + + # make sure geotop dataset contains kh and kv + if "kh" not in gt or "kv" not in gt: + if "kh" not in gt: + logger.info( + f"Calculating kh of geotop by multipying kv with an anisotropy of {anisotropy}" + ) + gt["kh"] = gt["kv"] * anisotropy + elif "kv" not in gt: + logger.info( + f"Calculating kv of geotop by dividing kh by an anisotropy of {anisotropy}" + ) + gt["kv"] = gt["kh"] / anisotropy + else: + # add kh and kv to gt + if geotop_k is None: + geotop_k = geotop.get_lithok_props() + gt = geotop.add_kh_and_kv(gt, geotop_k, anisotropy=anisotropy) + for layer in layers: # transform geotop data into layers gtl = geotop.to_model_layers(gt) @@ -234,9 +254,6 @@ def add_geotop_to_regis_layers( th = calculate_thickness(gtl) gtl = gtl.sel(layer=(th > 0).any(th.dims[1:])) - # add kh and kv to gt - gt = geotop.add_kh_and_kv(gt, geotop_k) - # add kh and kv from gt to gtl gtl = geotop.aggregate_to_ds(gt, gtl) From 0d9806ef52ccffc9459400649976853ea71176d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Fri, 25 Aug 2023 09:51:30 +0200 Subject: [PATCH 2/3] fix typos --- nlmod/read/regis.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py index 74b51a3f..6e836faf 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -204,8 +204,9 @@ def add_geotop_to_regis_layers( remove_nan_layers : bool, optional When True, layers with only 0 or NaN thickness are removed. The default is True. anisotropy : float, optional - The anisotropy value used when there are no kv values in df. The default is - 1.0. + + The anisotropy value (kh/kv) used when there are no kv values in df. The + default is 1.0. Returns ------- @@ -219,7 +220,7 @@ def add_geotop_to_regis_layers( if "kh" not in gt or "kv" not in gt: if "kh" not in gt: logger.info( - f"Calculating kh of geotop by multipying kv with an anisotropy of {anisotropy}" + f"Calculating kh of geotop by multiplying kv with an anisotropy of {anisotropy}" ) gt["kh"] = gt["kv"] * anisotropy elif "kv" not in gt: From 959d436c4ebe2bd4ab09942f1756482ed2b6a214 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Fri, 25 Aug 2023 10:06:45 +0200 Subject: [PATCH 3/3] Fix bug found by tests --- nlmod/read/regis.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py index 6e836faf..41036354 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -204,8 +204,7 @@ def add_geotop_to_regis_layers( remove_nan_layers : bool, optional When True, layers with only 0 or NaN thickness are removed. The default is True. anisotropy : float, optional - - The anisotropy value (kh/kv) used when there are no kv values in df. The + The anisotropy value (kh/kv) used when there are no kv values in df. The default is 1.0. Returns @@ -218,12 +217,12 @@ def add_geotop_to_regis_layers( # make sure geotop dataset contains kh and kv if "kh" not in gt or "kv" not in gt: - if "kh" not in gt: + if "kv" in gt: logger.info( f"Calculating kh of geotop by multiplying kv with an anisotropy of {anisotropy}" ) gt["kh"] = gt["kv"] * anisotropy - elif "kv" not in gt: + elif "kh" in gt: logger.info( f"Calculating kv of geotop by dividing kh by an anisotropy of {anisotropy}" )