Merge pull request #243 from ArtesiaWater/stop_overwriting_kh_and_kv_…

…geotop Stop overwriting kh and kv of geotop in add_geotop_to_regis_layers
gwmod · Aug 25, 2023 · 3e61b49 · 3e61b49
2 parents c73f524 + 959d436
commit 3e61b49
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 8 deletions.
diff --git 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
@@ -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
@@ -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)