diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py index c050c586..744ab665 100644 --- a/nlmod/dims/layers.py +++ b/nlmod/dims/layers.py @@ -827,7 +827,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 6314c53c..097cd6ce 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -183,7 +183,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. @@ -201,6 +201,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 (kh/kv) used when there are no kv values in df. The + default is 1.0. Returns ------- @@ -209,8 +212,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 "kv" in gt: + logger.info( + f"Calculating kh of geotop by multiplying kv with an anisotropy of {anisotropy}" + ) + gt["kh"] = gt["kv"] * anisotropy + elif "kh" 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) @@ -232,9 +252,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)