From 7d962fce76ed95b712967e11d6fb5c88c039ddc2 Mon Sep 17 00:00:00 2001
From: Joost Delsman <Joost.Delsman@deltares.nl>
Date: Tue, 15 Nov 2022 12:00:32 +0100
Subject: [PATCH 01/11] add ani hfb support for wq

---
 imod/wq/__init__.py |   2 +
 imod/wq/ani.py      | 145 ++++++++++++++++++++++++++++++++++++++++++++
 imod/wq/hfb.py      |  48 +++++++++++++++
 imod/wq/model.py    |   4 +-
 4 files changed, 197 insertions(+), 2 deletions(-)
 create mode 100644 imod/wq/ani.py
 create mode 100644 imod/wq/hfb.py

diff --git a/imod/wq/__init__.py b/imod/wq/__init__.py
index 6cbaf02ef..08e2d5b4f 100644
--- a/imod/wq/__init__.py
+++ b/imod/wq/__init__.py
@@ -14,6 +14,7 @@
     AdvectionModifiedMOC,
     AdvectionTVD,
 )
+from imod.wq.ani import HorizontalAnisotropy, HorizontalAnisotropyFile
 from imod.wq.bas import BasicFlow
 from imod.wq.btn import BasicTransport
 from imod.wq.chd import ConstantHead
@@ -26,6 +27,7 @@
     EvapotranspirationTopLayer,
 )
 from imod.wq.ghb import GeneralHeadBoundary
+from imod.wq.hfb import HorizontalFlowBarrier
 from imod.wq.lpf import LayerPropertyFlow
 from imod.wq.mal import MassLoading
 from imod.wq.model import SeawatModel
diff --git a/imod/wq/ani.py b/imod/wq/ani.py
new file mode 100644
index 000000000..9fc5a9ecd
--- /dev/null
+++ b/imod/wq/ani.py
@@ -0,0 +1,145 @@
+import numpy as np
+import pathlib
+import re
+import shutil
+
+import imod
+from imod.wq.pkgbase import Package
+
+class HorizontalAnisotropyFile(Package):
+    """
+    Horizontal Anisotropy package.
+
+    Parameters
+    ----------
+    anifile: str
+        is the file location of the imod-wq ani-file. This file contains the 
+        anisotropy factor and angle of each layer, either as a constant or a
+        reference to the file location of an '.arr' file. No checks are 
+        implemented for this file, user is responsible for consistency with 
+        model.
+    """
+
+    _pkg_id = "ani"
+
+    _template = (
+        "[ani]\n"
+        "    anifile = {anifile}\n"
+    )
+
+    def __init__(
+        self,
+        anifile,
+    ):
+        super().__init__()
+        self["anifile"] = anifile
+    
+    def _render(self, directory, *args, **kwargs):
+        path_ani = pathlib.Path(str(self['anifile'].values))
+        d = {"anifile": f"ani/{path_ani.name}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overload save function. 
+        Saves anifile to location, along with referenced .arr files"""
+        directory.mkdir(exist_ok=True)
+        
+        path_ani = pathlib.Path(str(self['anifile'].values))
+        
+        # regex adapted from stackoverflow: https://stackoverflow.com/questions/54990405/a-general-regex-to-extract-file-paths-not-urls
+        rgx = r'((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)'
+        with open(path_ani, "r") as f, open(directory / path_ani.name, "w") as f2:
+            for line in f:
+                p = re.search(rgx, line)
+                if p:
+                    # path to file detected, 
+                    # replace to relative and
+                    # copy to directory
+                    path = pathlib.Path(p[0])
+                    f2.write(line.replace(p[0], f"{directory.stem}/{path.name}"))
+                    if not path.is_absolute():
+                        path = path_ani.parent / path
+                    shutil.copyfile(path, directory / path.name)
+                else:
+                    f2.write(line)
+
+    def _pkgcheck(self, ibound=None):
+        pass
+
+
+class HorizontalAnisotropy(Package):
+    """
+    Horizontal Anisotropy package.
+    Anisotropy is a phenomenon for which the permeability k is not equal along the x- and y Cartesian axis.
+    
+    Parameters
+    ----------
+    factor : float or xr.DataArray of floats
+        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability). 
+        Factor between 0.0 (full anisotropic) and 1.0 (full isotropic).
+    angle : float or xr.DataArray of floats
+        The angle along the main principal axis (highest permeability) measured in degrees from north (0),
+        east (90), south (180) and west (270).
+    fn_ani : filename for created anifile, optional
+        Default: ani.ani
+    """
+
+    _pkg_id = "ani"
+
+    _template = (
+        "[ani]\n"
+        "    anifile = {anifile}\n"
+    )
+
+    def __init__(
+        self,
+        factor,
+        angle,
+        fn_ani="ani.ani"
+    ):
+        super().__init__()
+        self["factor"] = factor
+        self["factor"] = self["factor"].fillna(1.)
+        self["angle"] = angle
+        self["angle"] = self["angle"].fillna(0.)
+        self["fn_ani"] = fn_ani
+    
+    def _render(self, directory, *args, **kwargs):
+        d = {"anifile": f"ani/{str(self['fn_ani'].values)}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overload save function. 
+        Saves anifile to location, along with referenced .arr files"""
+        def _write(path, a, nodata=1.0e20, dtype=np.float32):
+            if not np.all(a == a[0][0]):
+                return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
+            else:
+                # write single value to ani file
+                pass
+
+        for name, da in self.dataset.data_vars.items():  # pylint: disable=no-member
+            if "y" in da.coords and "x" in da.coords:
+                path = pathlib.Path(directory).joinpath(f"{name}.arr")
+                if name == "factor":
+                    nodata = 1.
+                elif name == "angle":
+                    nodata = 0.
+                else:
+                    nodata = 1.0e20
+                imod.array_io.writing._save(path, da, nodata=nodata, pattern=None, dtype=np.float32, write=_write)
+
+        # write ani file
+        with open(directory / f"{str(self['fn_ani'].values)}", "w") as f:
+            for l in self.dataset.layer.values:
+                for prm in ["factor","angle"]:
+                    a = self.dataset[prm].sel(layer=l).values
+                    if not np.all(a == a[0][0]):
+                        f.write(f"open/close {directory.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n")
+                    else:
+                        f.write(f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n")
+
+    def _pkgcheck(self, ibound=None):
+        pass
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
new file mode 100644
index 000000000..712cb06b0
--- /dev/null
+++ b/imod/wq/hfb.py
@@ -0,0 +1,48 @@
+import pathlib
+import shutil
+from imod.wq.pkgbase import Package
+
+class HorizontalFlowBarrier(Package):
+    """
+    Horizontal Flow Barrier package.
+
+    Parameters
+    ----------
+    hfbfile: str
+        is the file location of the imod-wq hfb-file. This file contains cell- 
+        to-cell resistance values. The hfb file can be constructed from generate
+        files using imod-batch. No checks are implemented for this file, user is 
+        responsible for consistency with model.
+    """
+
+    _pkg_id = "hfb"
+
+    _template = (
+        "[hfb]\n"
+        "    hfbfile = {hfbfile}\n"
+    )
+
+    def __init__(
+        self,
+        hfbfile,
+    ):
+        super().__init__()
+        self["hfbfile"] = hfbfile
+    
+    def _render(self, directory, *args, **kwargs):
+        path_hfb = pathlib.Path(str(self["hfbfile"].values))
+        d = {"hfbfile": f"hfb/{path_hfb.name}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overloads save function. 
+        Saves hfbfile to directory"""
+        directory.mkdir(exist_ok=True)
+        
+        path_hfb = pathlib.Path(str(self["hfbfile"].values))
+        
+        shutil.copyfile(path_hfb, directory / path_hfb.name)
+
+    def _pkgcheck(self, ibound=None):
+        pass
diff --git a/imod/wq/model.py b/imod/wq/model.py
index 74a5ef9b7..d45cea69a 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -666,7 +666,7 @@ def render(self, directory, result_dir, writehelp):
                 directory=directory, globaltimes=globaltimes, nlayer=nlayer
             )
         )
-        for key in ("vdf", "adv", "dsp"):
+        for key in ("vdf", "adv", "dsp", "ani", "hfb"):
             content.append(
                 self._render_pkg(
                     key=key, directory=directory, globaltimes=globaltimes, nlayer=nlayer
@@ -793,7 +793,7 @@ def write(
             if (
                 "x" in pkg.dataset.coords
                 and "y" in pkg.dataset.coords
-                or pkg._pkg_id == "wel"
+                or pkg._pkg_id in ("wel", "ani", "hfb")
             ):
                 try:
                     pkg.save(directory=directory / pkgname)

From 4bdd4a1bd16e9c17e07f63a360777d0f1799642b Mon Sep 17 00:00:00 2001
From: Joost Delsman <Joost.Delsman@deltares.nl>
Date: Tue, 15 Nov 2022 14:12:53 +0100
Subject: [PATCH 02/11] store relative path in render to use in save

---
 imod/wq/ani.py   | 131 ++++++++++++++++++++++++++++-------------------
 imod/wq/hfb.py   |  33 ++++++------
 imod/wq/model.py |  19 ++++++-
 3 files changed, 113 insertions(+), 70 deletions(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 9fc5a9ecd..332b3369a 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -6,6 +6,7 @@
 import imod
 from imod.wq.pkgbase import Package
 
+
 class HorizontalAnisotropyFile(Package):
     """
     Horizontal Anisotropy package.
@@ -13,19 +14,16 @@ class HorizontalAnisotropyFile(Package):
     Parameters
     ----------
     anifile: str
-        is the file location of the imod-wq ani-file. This file contains the 
+        is the file location of the imod-wq ani-file. This file contains the
         anisotropy factor and angle of each layer, either as a constant or a
-        reference to the file location of an '.arr' file. No checks are 
-        implemented for this file, user is responsible for consistency with 
+        reference to the file location of an '.arr' file. No checks are
+        implemented for this file, user is responsible for consistency with
         model.
     """
 
     _pkg_id = "ani"
 
-    _template = (
-        "[ani]\n"
-        "    anifile = {anifile}\n"
-    )
+    _template = "[ani]\n" "    anifile = {anifile}\n\n"
 
     def __init__(
         self,
@@ -33,31 +31,40 @@ def __init__(
     ):
         super().__init__()
         self["anifile"] = anifile
-    
-    def _render(self, directory, *args, **kwargs):
-        path_ani = pathlib.Path(str(self['anifile'].values))
-        d = {"anifile": f"ani/{path_ani.name}"}
+
+    def _render(self, modelname, directory, nlayer):
+        # write ani file
+        # in render function, as here we know where the model is run from
+        # and how many layers: store this info for later use in save
+        self.anifile = f"{modelname}.ani"
+        self.rendir = directory
+        self.nlayer = nlayer
+
+        d = {"anifile": f"{directory.as_posix()}/{modelname}.ani"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overload save function. 
-        Saves anifile to location, along with referenced .arr files"""
-        directory.mkdir(exist_ok=True)
-        
-        path_ani = pathlib.Path(str(self['anifile'].values))
-        
+        """Overload save function.
+        Saves anifile to location, along with referenced .arr files
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
+        path_ani = pathlib.Path(str(self["anifile"].values))
+
         # regex adapted from stackoverflow: https://stackoverflow.com/questions/54990405/a-general-regex-to-extract-file-paths-not-urls
-        rgx = r'((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)'
-        with open(path_ani, "r") as f, open(directory / path_ani.name, "w") as f2:
+        rgx = r"((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)"
+        with open(path_ani, "r") as f, open(directory / self.anifile, "w") as f2:
             for line in f:
                 p = re.search(rgx, line)
                 if p:
-                    # path to file detected, 
+                    # path to file detected,
                     # replace to relative and
                     # copy to directory
                     path = pathlib.Path(p[0])
-                    f2.write(line.replace(p[0], f"{directory.stem}/{path.name}"))
+                    f2.write(
+                        line.replace(p[0], f"{self.rendir.as_posix()}/{path.name}")
+                    )
                     if not path.is_absolute():
                         path = path_ani.parent / path
                     shutil.copyfile(path, directory / path.name)
@@ -72,47 +79,50 @@ class HorizontalAnisotropy(Package):
     """
     Horizontal Anisotropy package.
     Anisotropy is a phenomenon for which the permeability k is not equal along the x- and y Cartesian axis.
-    
+
     Parameters
     ----------
     factor : float or xr.DataArray of floats
-        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability). 
+        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability).
         Factor between 0.0 (full anisotropic) and 1.0 (full isotropic).
     angle : float or xr.DataArray of floats
         The angle along the main principal axis (highest permeability) measured in degrees from north (0),
         east (90), south (180) and west (270).
-    fn_ani : filename for created anifile, optional
-        Default: ani.ani
     """
 
     _pkg_id = "ani"
 
-    _template = (
-        "[ani]\n"
-        "    anifile = {anifile}\n"
-    )
+    _template = "[ani]\n" "    anifile = {anifile}\n\n"
 
     def __init__(
         self,
         factor,
         angle,
-        fn_ani="ani.ani"
     ):
         super().__init__()
         self["factor"] = factor
-        self["factor"] = self["factor"].fillna(1.)
+        self["factor"] = self["factor"].fillna(1.0)
         self["angle"] = angle
-        self["angle"] = self["angle"].fillna(0.)
-        self["fn_ani"] = fn_ani
-    
-    def _render(self, directory, *args, **kwargs):
-        d = {"anifile": f"ani/{str(self['fn_ani'].values)}"}
+        self["angle"] = self["angle"].fillna(0.0)
+
+    def _render(self, modelname, directory, nlayer):
+        # write ani file
+        # in render function, as here we know where the model is run from
+        # and how many layers: store this info for later use in save
+        self.anifile = f"{modelname}.ani"
+        self.rendir = directory
+        self.nlayer = nlayer
+
+        d = {"anifile": f"{directory.as_posix()}/{modelname}.ani"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overload save function. 
-        Saves anifile to location, along with referenced .arr files"""
+        """Overload save function.
+        Saves anifile to location, along with created .arr files
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
         def _write(path, a, nodata=1.0e20, dtype=np.float32):
             if not np.all(a == a[0][0]):
                 return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
@@ -124,22 +134,39 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
             if "y" in da.coords and "x" in da.coords:
                 path = pathlib.Path(directory).joinpath(f"{name}.arr")
                 if name == "factor":
-                    nodata = 1.
+                    nodata = 1.0
                 elif name == "angle":
-                    nodata = 0.
+                    nodata = 0.0
                 else:
                     nodata = 1.0e20
-                imod.array_io.writing._save(path, da, nodata=nodata, pattern=None, dtype=np.float32, write=_write)
-
-        # write ani file
-        with open(directory / f"{str(self['fn_ani'].values)}", "w") as f:
-            for l in self.dataset.layer.values:
-                for prm in ["factor","angle"]:
-                    a = self.dataset[prm].sel(layer=l).values
-                    if not np.all(a == a[0][0]):
-                        f.write(f"open/close {directory.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n")
+                imod.array_io.writing._save(
+                    path,
+                    da,
+                    nodata=nodata,
+                    pattern=None,
+                    dtype=np.float32,
+                    write=_write,
+                )
+
+        # save anifile with data stored during _render
+        with open(directory / self.anifile, "w") as f:
+            for l in range(1, self.nlayer + 1):
+                for prm in ["factor", "angle"]:
+                    da = self.dataset[prm]
+                    if "layer" in da.coords and "y" in da.coords and "x" in da.coords:
+                        a = da.sel(layer=l).values
+                        if not np.all(a == a[0][0]):
+                            f.write(
+                                f"open/close {self.rendir.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n"
+                            )
+                        else:
+                            f.write(
+                                f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n"
+                            )
                     else:
-                        f.write(f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n")
+                        f.write(
+                            f"constant {float(self.dataset[prm].values):.5f} {prm}_l{float(l):.0f}\n"
+                        )
 
     def _pkgcheck(self, ibound=None):
         pass
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 712cb06b0..5dd069454 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -2,6 +2,7 @@
 import shutil
 from imod.wq.pkgbase import Package
 
+
 class HorizontalFlowBarrier(Package):
     """
     Horizontal Flow Barrier package.
@@ -9,18 +10,15 @@ class HorizontalFlowBarrier(Package):
     Parameters
     ----------
     hfbfile: str
-        is the file location of the imod-wq hfb-file. This file contains cell- 
+        is the file location of the imod-wq hfb-file. This file contains cell-
         to-cell resistance values. The hfb file can be constructed from generate
-        files using imod-batch. No checks are implemented for this file, user is 
+        files using imod-batch. No checks are implemented for this file, user is
         responsible for consistency with model.
     """
 
     _pkg_id = "hfb"
 
-    _template = (
-        "[hfb]\n"
-        "    hfbfile = {hfbfile}\n"
-    )
+    _template = "[hfb]\n" "    hfbfile = {hfbfile}\n\n"
 
     def __init__(
         self,
@@ -28,21 +26,22 @@ def __init__(
     ):
         super().__init__()
         self["hfbfile"] = hfbfile
-    
-    def _render(self, directory, *args, **kwargs):
-        path_hfb = pathlib.Path(str(self["hfbfile"].values))
-        d = {"hfbfile": f"hfb/{path_hfb.name}"}
+
+    def _render(self, modelname, directory, *args, **kwargs):
+        self.hfbfile = f"{modelname}.hfb"
+        d = {"hfbfile": f"{directory.as_posix()}/{modelname}.hfb"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overloads save function. 
-        Saves hfbfile to directory"""
-        directory.mkdir(exist_ok=True)
-        
+        """Overloads save function.
+        Saves hfbfile to directory
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
         path_hfb = pathlib.Path(str(self["hfbfile"].values))
-        
-        shutil.copyfile(path_hfb, directory / path_hfb.name)
+
+        shutil.copyfile(path_hfb, directory / self.hfbfile)
 
     def _pkgcheck(self, ibound=None):
         pass
diff --git a/imod/wq/model.py b/imod/wq/model.py
index d45cea69a..07bde9acb 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -605,6 +605,16 @@ def _bas_btn_rch_evt_mal_tvc_sinkssources(self):
 
         return n_extra
 
+    def _render_anihfb(self, modelname, directory, nlayer):
+        out = ""
+        for key in ("ani", "hfb"):
+            pkgkey = self._get_pkgkey(key)
+            if pkgkey is not None:
+                out += self[pkgkey]._render(
+                    modelname=modelname, directory=directory / pkgkey, nlayer=nlayer
+                )
+        return out
+
     def render(self, directory, result_dir, writehelp):
         """
         Render the runfile as a string, package by package.
@@ -636,6 +646,13 @@ def render(self, directory, result_dir, writehelp):
                 )
             )
 
+        # ani and hfb packages
+        content.append(
+            self._render_anihfb(
+                modelname=self.modelname, directory=directory, nlayer=nlayer
+            )
+        )
+
         # multi-system package group: chd, drn, ghb, riv, wel
         modflowcontent, ssm_content, n_sinkssources = self._render_groups(
             directory=directory, globaltimes=globaltimes
@@ -666,7 +683,7 @@ def render(self, directory, result_dir, writehelp):
                 directory=directory, globaltimes=globaltimes, nlayer=nlayer
             )
         )
-        for key in ("vdf", "adv", "dsp", "ani", "hfb"):
+        for key in ("vdf", "adv", "dsp"):
             content.append(
                 self._render_pkg(
                     key=key, directory=directory, globaltimes=globaltimes, nlayer=nlayer

From 6574710e60a95445e6604049a7f727fecdae088e Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Thu, 17 Nov 2022 12:15:14 +0100
Subject: [PATCH 03/11] fix constant or write arr

---
 imod/wq/ani.py | 39 ++++++++++++++++++++++++---------------
 1 file changed, 24 insertions(+), 15 deletions(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 332b3369a..965030860 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -101,9 +101,7 @@ def __init__(
     ):
         super().__init__()
         self["factor"] = factor
-        self["factor"] = self["factor"].fillna(1.0)
         self["angle"] = angle
-        self["angle"] = self["angle"].fillna(0.0)
 
     def _render(self, modelname, directory, nlayer):
         # write ani file
@@ -123,9 +121,22 @@ def save(self, directory):
         assumes _render() to have run previously"""
         directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
 
-        def _write(path, a, nodata=1.0e20, dtype=np.float32):
-            if not np.all(a == a[0][0]):
-                return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
+        nodata_val = {"factor": 1.0, "angle": 0.0}
+
+        def _check_all_equal(da):
+            return np.all(np.isnan(da)) or np.all(
+                da.values[~np.isnan(da)] == da.values[~np.isnan(da)][0]
+            )
+
+        def _write(path, da, nodata=1.0e20, dtype=np.float32):
+            if not _check_all_equal(da):
+                dx, xmin, xmax, dy, ymin, ymax = imod.util.spatial_reference(da)
+                ncol, nrow = da.shape
+                footer = f" DIMENSIONS\n{ncol}\n{nrow}\n{xmin}\n{ymin}\n{xmax}\n{ymax}\n{nodata}\n0\n{dx}\n{dx}"
+                a = np.nan_to_num(da.values, nan=nodata)
+                return np.savetxt(
+                    path, a, fmt="%.5f", delimiter=" ", footer=footer, comments=""
+                )
             else:
                 # write single value to ani file
                 pass
@@ -133,16 +144,10 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
         for name, da in self.dataset.data_vars.items():  # pylint: disable=no-member
             if "y" in da.coords and "x" in da.coords:
                 path = pathlib.Path(directory).joinpath(f"{name}.arr")
-                if name == "factor":
-                    nodata = 1.0
-                elif name == "angle":
-                    nodata = 0.0
-                else:
-                    nodata = 1.0e20
                 imod.array_io.writing._save(
                     path,
                     da,
-                    nodata=nodata,
+                    nodata=nodata_val[name],
                     pattern=None,
                     dtype=np.float32,
                     write=_write,
@@ -154,14 +159,18 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
                 for prm in ["factor", "angle"]:
                     da = self.dataset[prm]
                     if "layer" in da.coords and "y" in da.coords and "x" in da.coords:
-                        a = da.sel(layer=l).values
-                        if not np.all(a == a[0][0]):
+                        a = da.sel(layer=l)
+                        if not _check_all_equal(a):
                             f.write(
                                 f"open/close {self.rendir.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n"
                             )
                         else:
+                            if np.all(np.isnan(a)):
+                                val = nodata_val[prm]
+                            else:
+                                val = a[~np.isnan()][0]
                             f.write(
-                                f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n"
+                                f"constant {float(val):.5f} {prm}_l{float(l):.0f}\n"
                             )
                     else:
                         f.write(

From 6ecbf3986de5816fa9ecb839772047cf339d79d7 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Thu, 17 Nov 2022 13:40:12 +0100
Subject: [PATCH 04/11] hfb id changed to hfb6 as reqd by imodwq

---
 imod/wq/hfb.py   | 4 ++--
 imod/wq/model.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 5dd069454..927e835e7 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -16,9 +16,9 @@ class HorizontalFlowBarrier(Package):
         responsible for consistency with model.
     """
 
-    _pkg_id = "hfb"
+    _pkg_id = "hfb6"
 
-    _template = "[hfb]\n" "    hfbfile = {hfbfile}\n\n"
+    _template = "[hfb6]\n" "    hfbfile = {hfbfile}\n\n"
 
     def __init__(
         self,
diff --git a/imod/wq/model.py b/imod/wq/model.py
index 07bde9acb..efdd24e1e 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -607,7 +607,7 @@ def _bas_btn_rch_evt_mal_tvc_sinkssources(self):
 
     def _render_anihfb(self, modelname, directory, nlayer):
         out = ""
-        for key in ("ani", "hfb"):
+        for key in ("ani", "hfb6"):
             pkgkey = self._get_pkgkey(key)
             if pkgkey is not None:
                 out += self[pkgkey]._render(
@@ -810,7 +810,7 @@ def write(
             if (
                 "x" in pkg.dataset.coords
                 and "y" in pkg.dataset.coords
-                or pkg._pkg_id in ("wel", "ani", "hfb")
+                or pkg._pkg_id in ("wel", "ani", "hfb6")
             ):
                 try:
                     pkg.save(directory=directory / pkgname)

From 9842e272f2d46089e6de6e83aedf0b6bd4831a1c Mon Sep 17 00:00:00 2001
From: Joost Delsman <Joost.Delsman@deltares.nl>
Date: Tue, 15 Nov 2022 12:00:32 +0100
Subject: [PATCH 05/11] add ani hfb support for wq

---
 imod/wq/__init__.py |   2 +
 imod/wq/ani.py      | 145 ++++++++++++++++++++++++++++++++++++++++++++
 imod/wq/hfb.py      |  48 +++++++++++++++
 imod/wq/model.py    |   4 +-
 4 files changed, 197 insertions(+), 2 deletions(-)
 create mode 100644 imod/wq/ani.py
 create mode 100644 imod/wq/hfb.py

diff --git a/imod/wq/__init__.py b/imod/wq/__init__.py
index 6cbaf02ef..08e2d5b4f 100644
--- a/imod/wq/__init__.py
+++ b/imod/wq/__init__.py
@@ -14,6 +14,7 @@
     AdvectionModifiedMOC,
     AdvectionTVD,
 )
+from imod.wq.ani import HorizontalAnisotropy, HorizontalAnisotropyFile
 from imod.wq.bas import BasicFlow
 from imod.wq.btn import BasicTransport
 from imod.wq.chd import ConstantHead
@@ -26,6 +27,7 @@
     EvapotranspirationTopLayer,
 )
 from imod.wq.ghb import GeneralHeadBoundary
+from imod.wq.hfb import HorizontalFlowBarrier
 from imod.wq.lpf import LayerPropertyFlow
 from imod.wq.mal import MassLoading
 from imod.wq.model import SeawatModel
diff --git a/imod/wq/ani.py b/imod/wq/ani.py
new file mode 100644
index 000000000..9fc5a9ecd
--- /dev/null
+++ b/imod/wq/ani.py
@@ -0,0 +1,145 @@
+import numpy as np
+import pathlib
+import re
+import shutil
+
+import imod
+from imod.wq.pkgbase import Package
+
+class HorizontalAnisotropyFile(Package):
+    """
+    Horizontal Anisotropy package.
+
+    Parameters
+    ----------
+    anifile: str
+        is the file location of the imod-wq ani-file. This file contains the 
+        anisotropy factor and angle of each layer, either as a constant or a
+        reference to the file location of an '.arr' file. No checks are 
+        implemented for this file, user is responsible for consistency with 
+        model.
+    """
+
+    _pkg_id = "ani"
+
+    _template = (
+        "[ani]\n"
+        "    anifile = {anifile}\n"
+    )
+
+    def __init__(
+        self,
+        anifile,
+    ):
+        super().__init__()
+        self["anifile"] = anifile
+    
+    def _render(self, directory, *args, **kwargs):
+        path_ani = pathlib.Path(str(self['anifile'].values))
+        d = {"anifile": f"ani/{path_ani.name}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overload save function. 
+        Saves anifile to location, along with referenced .arr files"""
+        directory.mkdir(exist_ok=True)
+        
+        path_ani = pathlib.Path(str(self['anifile'].values))
+        
+        # regex adapted from stackoverflow: https://stackoverflow.com/questions/54990405/a-general-regex-to-extract-file-paths-not-urls
+        rgx = r'((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)'
+        with open(path_ani, "r") as f, open(directory / path_ani.name, "w") as f2:
+            for line in f:
+                p = re.search(rgx, line)
+                if p:
+                    # path to file detected, 
+                    # replace to relative and
+                    # copy to directory
+                    path = pathlib.Path(p[0])
+                    f2.write(line.replace(p[0], f"{directory.stem}/{path.name}"))
+                    if not path.is_absolute():
+                        path = path_ani.parent / path
+                    shutil.copyfile(path, directory / path.name)
+                else:
+                    f2.write(line)
+
+    def _pkgcheck(self, ibound=None):
+        pass
+
+
+class HorizontalAnisotropy(Package):
+    """
+    Horizontal Anisotropy package.
+    Anisotropy is a phenomenon for which the permeability k is not equal along the x- and y Cartesian axis.
+    
+    Parameters
+    ----------
+    factor : float or xr.DataArray of floats
+        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability). 
+        Factor between 0.0 (full anisotropic) and 1.0 (full isotropic).
+    angle : float or xr.DataArray of floats
+        The angle along the main principal axis (highest permeability) measured in degrees from north (0),
+        east (90), south (180) and west (270).
+    fn_ani : filename for created anifile, optional
+        Default: ani.ani
+    """
+
+    _pkg_id = "ani"
+
+    _template = (
+        "[ani]\n"
+        "    anifile = {anifile}\n"
+    )
+
+    def __init__(
+        self,
+        factor,
+        angle,
+        fn_ani="ani.ani"
+    ):
+        super().__init__()
+        self["factor"] = factor
+        self["factor"] = self["factor"].fillna(1.)
+        self["angle"] = angle
+        self["angle"] = self["angle"].fillna(0.)
+        self["fn_ani"] = fn_ani
+    
+    def _render(self, directory, *args, **kwargs):
+        d = {"anifile": f"ani/{str(self['fn_ani'].values)}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overload save function. 
+        Saves anifile to location, along with referenced .arr files"""
+        def _write(path, a, nodata=1.0e20, dtype=np.float32):
+            if not np.all(a == a[0][0]):
+                return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
+            else:
+                # write single value to ani file
+                pass
+
+        for name, da in self.dataset.data_vars.items():  # pylint: disable=no-member
+            if "y" in da.coords and "x" in da.coords:
+                path = pathlib.Path(directory).joinpath(f"{name}.arr")
+                if name == "factor":
+                    nodata = 1.
+                elif name == "angle":
+                    nodata = 0.
+                else:
+                    nodata = 1.0e20
+                imod.array_io.writing._save(path, da, nodata=nodata, pattern=None, dtype=np.float32, write=_write)
+
+        # write ani file
+        with open(directory / f"{str(self['fn_ani'].values)}", "w") as f:
+            for l in self.dataset.layer.values:
+                for prm in ["factor","angle"]:
+                    a = self.dataset[prm].sel(layer=l).values
+                    if not np.all(a == a[0][0]):
+                        f.write(f"open/close {directory.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n")
+                    else:
+                        f.write(f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n")
+
+    def _pkgcheck(self, ibound=None):
+        pass
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
new file mode 100644
index 000000000..712cb06b0
--- /dev/null
+++ b/imod/wq/hfb.py
@@ -0,0 +1,48 @@
+import pathlib
+import shutil
+from imod.wq.pkgbase import Package
+
+class HorizontalFlowBarrier(Package):
+    """
+    Horizontal Flow Barrier package.
+
+    Parameters
+    ----------
+    hfbfile: str
+        is the file location of the imod-wq hfb-file. This file contains cell- 
+        to-cell resistance values. The hfb file can be constructed from generate
+        files using imod-batch. No checks are implemented for this file, user is 
+        responsible for consistency with model.
+    """
+
+    _pkg_id = "hfb"
+
+    _template = (
+        "[hfb]\n"
+        "    hfbfile = {hfbfile}\n"
+    )
+
+    def __init__(
+        self,
+        hfbfile,
+    ):
+        super().__init__()
+        self["hfbfile"] = hfbfile
+    
+    def _render(self, directory, *args, **kwargs):
+        path_hfb = pathlib.Path(str(self["hfbfile"].values))
+        d = {"hfbfile": f"hfb/{path_hfb.name}"}
+
+        return self._template.format(**d)
+            
+    def save(self, directory):
+        """Overloads save function. 
+        Saves hfbfile to directory"""
+        directory.mkdir(exist_ok=True)
+        
+        path_hfb = pathlib.Path(str(self["hfbfile"].values))
+        
+        shutil.copyfile(path_hfb, directory / path_hfb.name)
+
+    def _pkgcheck(self, ibound=None):
+        pass
diff --git a/imod/wq/model.py b/imod/wq/model.py
index 74a5ef9b7..d45cea69a 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -666,7 +666,7 @@ def render(self, directory, result_dir, writehelp):
                 directory=directory, globaltimes=globaltimes, nlayer=nlayer
             )
         )
-        for key in ("vdf", "adv", "dsp"):
+        for key in ("vdf", "adv", "dsp", "ani", "hfb"):
             content.append(
                 self._render_pkg(
                     key=key, directory=directory, globaltimes=globaltimes, nlayer=nlayer
@@ -793,7 +793,7 @@ def write(
             if (
                 "x" in pkg.dataset.coords
                 and "y" in pkg.dataset.coords
-                or pkg._pkg_id == "wel"
+                or pkg._pkg_id in ("wel", "ani", "hfb")
             ):
                 try:
                     pkg.save(directory=directory / pkgname)

From ca3945d832b66e36b48a42c7f78406b693999178 Mon Sep 17 00:00:00 2001
From: Joost Delsman <Joost.Delsman@deltares.nl>
Date: Tue, 15 Nov 2022 14:12:53 +0100
Subject: [PATCH 06/11] store relative path in render to use in save

---
 imod/wq/ani.py   | 131 ++++++++++++++++++++++++++++-------------------
 imod/wq/hfb.py   |  33 ++++++------
 imod/wq/model.py |  19 ++++++-
 3 files changed, 113 insertions(+), 70 deletions(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 9fc5a9ecd..332b3369a 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -6,6 +6,7 @@
 import imod
 from imod.wq.pkgbase import Package
 
+
 class HorizontalAnisotropyFile(Package):
     """
     Horizontal Anisotropy package.
@@ -13,19 +14,16 @@ class HorizontalAnisotropyFile(Package):
     Parameters
     ----------
     anifile: str
-        is the file location of the imod-wq ani-file. This file contains the 
+        is the file location of the imod-wq ani-file. This file contains the
         anisotropy factor and angle of each layer, either as a constant or a
-        reference to the file location of an '.arr' file. No checks are 
-        implemented for this file, user is responsible for consistency with 
+        reference to the file location of an '.arr' file. No checks are
+        implemented for this file, user is responsible for consistency with
         model.
     """
 
     _pkg_id = "ani"
 
-    _template = (
-        "[ani]\n"
-        "    anifile = {anifile}\n"
-    )
+    _template = "[ani]\n" "    anifile = {anifile}\n\n"
 
     def __init__(
         self,
@@ -33,31 +31,40 @@ def __init__(
     ):
         super().__init__()
         self["anifile"] = anifile
-    
-    def _render(self, directory, *args, **kwargs):
-        path_ani = pathlib.Path(str(self['anifile'].values))
-        d = {"anifile": f"ani/{path_ani.name}"}
+
+    def _render(self, modelname, directory, nlayer):
+        # write ani file
+        # in render function, as here we know where the model is run from
+        # and how many layers: store this info for later use in save
+        self.anifile = f"{modelname}.ani"
+        self.rendir = directory
+        self.nlayer = nlayer
+
+        d = {"anifile": f"{directory.as_posix()}/{modelname}.ani"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overload save function. 
-        Saves anifile to location, along with referenced .arr files"""
-        directory.mkdir(exist_ok=True)
-        
-        path_ani = pathlib.Path(str(self['anifile'].values))
-        
+        """Overload save function.
+        Saves anifile to location, along with referenced .arr files
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
+        path_ani = pathlib.Path(str(self["anifile"].values))
+
         # regex adapted from stackoverflow: https://stackoverflow.com/questions/54990405/a-general-regex-to-extract-file-paths-not-urls
-        rgx = r'((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)'
-        with open(path_ani, "r") as f, open(directory / path_ani.name, "w") as f2:
+        rgx = r"((?:[a-zA-Z]:|(?<![:/\\])[\\\/](?!CLOSE)(?!close )|\~[\\\/]|(?:\.{1,2}[\\\/])+)[\w+\\\s_\-\(\)\/]*(?:\.\w+)*)"
+        with open(path_ani, "r") as f, open(directory / self.anifile, "w") as f2:
             for line in f:
                 p = re.search(rgx, line)
                 if p:
-                    # path to file detected, 
+                    # path to file detected,
                     # replace to relative and
                     # copy to directory
                     path = pathlib.Path(p[0])
-                    f2.write(line.replace(p[0], f"{directory.stem}/{path.name}"))
+                    f2.write(
+                        line.replace(p[0], f"{self.rendir.as_posix()}/{path.name}")
+                    )
                     if not path.is_absolute():
                         path = path_ani.parent / path
                     shutil.copyfile(path, directory / path.name)
@@ -72,47 +79,50 @@ class HorizontalAnisotropy(Package):
     """
     Horizontal Anisotropy package.
     Anisotropy is a phenomenon for which the permeability k is not equal along the x- and y Cartesian axis.
-    
+
     Parameters
     ----------
     factor : float or xr.DataArray of floats
-        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability). 
+        The anisotropic factor perpendicular to the main principal axis (axis of highest permeability).
         Factor between 0.0 (full anisotropic) and 1.0 (full isotropic).
     angle : float or xr.DataArray of floats
         The angle along the main principal axis (highest permeability) measured in degrees from north (0),
         east (90), south (180) and west (270).
-    fn_ani : filename for created anifile, optional
-        Default: ani.ani
     """
 
     _pkg_id = "ani"
 
-    _template = (
-        "[ani]\n"
-        "    anifile = {anifile}\n"
-    )
+    _template = "[ani]\n" "    anifile = {anifile}\n\n"
 
     def __init__(
         self,
         factor,
         angle,
-        fn_ani="ani.ani"
     ):
         super().__init__()
         self["factor"] = factor
-        self["factor"] = self["factor"].fillna(1.)
+        self["factor"] = self["factor"].fillna(1.0)
         self["angle"] = angle
-        self["angle"] = self["angle"].fillna(0.)
-        self["fn_ani"] = fn_ani
-    
-    def _render(self, directory, *args, **kwargs):
-        d = {"anifile": f"ani/{str(self['fn_ani'].values)}"}
+        self["angle"] = self["angle"].fillna(0.0)
+
+    def _render(self, modelname, directory, nlayer):
+        # write ani file
+        # in render function, as here we know where the model is run from
+        # and how many layers: store this info for later use in save
+        self.anifile = f"{modelname}.ani"
+        self.rendir = directory
+        self.nlayer = nlayer
+
+        d = {"anifile": f"{directory.as_posix()}/{modelname}.ani"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overload save function. 
-        Saves anifile to location, along with referenced .arr files"""
+        """Overload save function.
+        Saves anifile to location, along with created .arr files
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
         def _write(path, a, nodata=1.0e20, dtype=np.float32):
             if not np.all(a == a[0][0]):
                 return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
@@ -124,22 +134,39 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
             if "y" in da.coords and "x" in da.coords:
                 path = pathlib.Path(directory).joinpath(f"{name}.arr")
                 if name == "factor":
-                    nodata = 1.
+                    nodata = 1.0
                 elif name == "angle":
-                    nodata = 0.
+                    nodata = 0.0
                 else:
                     nodata = 1.0e20
-                imod.array_io.writing._save(path, da, nodata=nodata, pattern=None, dtype=np.float32, write=_write)
-
-        # write ani file
-        with open(directory / f"{str(self['fn_ani'].values)}", "w") as f:
-            for l in self.dataset.layer.values:
-                for prm in ["factor","angle"]:
-                    a = self.dataset[prm].sel(layer=l).values
-                    if not np.all(a == a[0][0]):
-                        f.write(f"open/close {directory.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n")
+                imod.array_io.writing._save(
+                    path,
+                    da,
+                    nodata=nodata,
+                    pattern=None,
+                    dtype=np.float32,
+                    write=_write,
+                )
+
+        # save anifile with data stored during _render
+        with open(directory / self.anifile, "w") as f:
+            for l in range(1, self.nlayer + 1):
+                for prm in ["factor", "angle"]:
+                    da = self.dataset[prm]
+                    if "layer" in da.coords and "y" in da.coords and "x" in da.coords:
+                        a = da.sel(layer=l).values
+                        if not np.all(a == a[0][0]):
+                            f.write(
+                                f"open/close {self.rendir.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n"
+                            )
+                        else:
+                            f.write(
+                                f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n"
+                            )
                     else:
-                        f.write(f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n")
+                        f.write(
+                            f"constant {float(self.dataset[prm].values):.5f} {prm}_l{float(l):.0f}\n"
+                        )
 
     def _pkgcheck(self, ibound=None):
         pass
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 712cb06b0..5dd069454 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -2,6 +2,7 @@
 import shutil
 from imod.wq.pkgbase import Package
 
+
 class HorizontalFlowBarrier(Package):
     """
     Horizontal Flow Barrier package.
@@ -9,18 +10,15 @@ class HorizontalFlowBarrier(Package):
     Parameters
     ----------
     hfbfile: str
-        is the file location of the imod-wq hfb-file. This file contains cell- 
+        is the file location of the imod-wq hfb-file. This file contains cell-
         to-cell resistance values. The hfb file can be constructed from generate
-        files using imod-batch. No checks are implemented for this file, user is 
+        files using imod-batch. No checks are implemented for this file, user is
         responsible for consistency with model.
     """
 
     _pkg_id = "hfb"
 
-    _template = (
-        "[hfb]\n"
-        "    hfbfile = {hfbfile}\n"
-    )
+    _template = "[hfb]\n" "    hfbfile = {hfbfile}\n\n"
 
     def __init__(
         self,
@@ -28,21 +26,22 @@ def __init__(
     ):
         super().__init__()
         self["hfbfile"] = hfbfile
-    
-    def _render(self, directory, *args, **kwargs):
-        path_hfb = pathlib.Path(str(self["hfbfile"].values))
-        d = {"hfbfile": f"hfb/{path_hfb.name}"}
+
+    def _render(self, modelname, directory, *args, **kwargs):
+        self.hfbfile = f"{modelname}.hfb"
+        d = {"hfbfile": f"{directory.as_posix()}/{modelname}.hfb"}
 
         return self._template.format(**d)
-            
+
     def save(self, directory):
-        """Overloads save function. 
-        Saves hfbfile to directory"""
-        directory.mkdir(exist_ok=True)
-        
+        """Overloads save function.
+        Saves hfbfile to directory
+        assumes _render() to have run previously"""
+        directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
+
         path_hfb = pathlib.Path(str(self["hfbfile"].values))
-        
-        shutil.copyfile(path_hfb, directory / path_hfb.name)
+
+        shutil.copyfile(path_hfb, directory / self.hfbfile)
 
     def _pkgcheck(self, ibound=None):
         pass
diff --git a/imod/wq/model.py b/imod/wq/model.py
index d45cea69a..07bde9acb 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -605,6 +605,16 @@ def _bas_btn_rch_evt_mal_tvc_sinkssources(self):
 
         return n_extra
 
+    def _render_anihfb(self, modelname, directory, nlayer):
+        out = ""
+        for key in ("ani", "hfb"):
+            pkgkey = self._get_pkgkey(key)
+            if pkgkey is not None:
+                out += self[pkgkey]._render(
+                    modelname=modelname, directory=directory / pkgkey, nlayer=nlayer
+                )
+        return out
+
     def render(self, directory, result_dir, writehelp):
         """
         Render the runfile as a string, package by package.
@@ -636,6 +646,13 @@ def render(self, directory, result_dir, writehelp):
                 )
             )
 
+        # ani and hfb packages
+        content.append(
+            self._render_anihfb(
+                modelname=self.modelname, directory=directory, nlayer=nlayer
+            )
+        )
+
         # multi-system package group: chd, drn, ghb, riv, wel
         modflowcontent, ssm_content, n_sinkssources = self._render_groups(
             directory=directory, globaltimes=globaltimes
@@ -666,7 +683,7 @@ def render(self, directory, result_dir, writehelp):
                 directory=directory, globaltimes=globaltimes, nlayer=nlayer
             )
         )
-        for key in ("vdf", "adv", "dsp", "ani", "hfb"):
+        for key in ("vdf", "adv", "dsp"):
             content.append(
                 self._render_pkg(
                     key=key, directory=directory, globaltimes=globaltimes, nlayer=nlayer

From d6c83c625d53fdc470c2625fdbb980e9f7698735 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Thu, 17 Nov 2022 12:15:14 +0100
Subject: [PATCH 07/11] fix constant or write arr

---
 imod/wq/ani.py | 39 ++++++++++++++++++++++++---------------
 1 file changed, 24 insertions(+), 15 deletions(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 332b3369a..965030860 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -101,9 +101,7 @@ def __init__(
     ):
         super().__init__()
         self["factor"] = factor
-        self["factor"] = self["factor"].fillna(1.0)
         self["angle"] = angle
-        self["angle"] = self["angle"].fillna(0.0)
 
     def _render(self, modelname, directory, nlayer):
         # write ani file
@@ -123,9 +121,22 @@ def save(self, directory):
         assumes _render() to have run previously"""
         directory.mkdir(exist_ok=True)  # otherwise handled by idf.save
 
-        def _write(path, a, nodata=1.0e20, dtype=np.float32):
-            if not np.all(a == a[0][0]):
-                return np.savetxt(path, a, fmt="%.5f", delimiter=" ")
+        nodata_val = {"factor": 1.0, "angle": 0.0}
+
+        def _check_all_equal(da):
+            return np.all(np.isnan(da)) or np.all(
+                da.values[~np.isnan(da)] == da.values[~np.isnan(da)][0]
+            )
+
+        def _write(path, da, nodata=1.0e20, dtype=np.float32):
+            if not _check_all_equal(da):
+                dx, xmin, xmax, dy, ymin, ymax = imod.util.spatial_reference(da)
+                ncol, nrow = da.shape
+                footer = f" DIMENSIONS\n{ncol}\n{nrow}\n{xmin}\n{ymin}\n{xmax}\n{ymax}\n{nodata}\n0\n{dx}\n{dx}"
+                a = np.nan_to_num(da.values, nan=nodata)
+                return np.savetxt(
+                    path, a, fmt="%.5f", delimiter=" ", footer=footer, comments=""
+                )
             else:
                 # write single value to ani file
                 pass
@@ -133,16 +144,10 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
         for name, da in self.dataset.data_vars.items():  # pylint: disable=no-member
             if "y" in da.coords and "x" in da.coords:
                 path = pathlib.Path(directory).joinpath(f"{name}.arr")
-                if name == "factor":
-                    nodata = 1.0
-                elif name == "angle":
-                    nodata = 0.0
-                else:
-                    nodata = 1.0e20
                 imod.array_io.writing._save(
                     path,
                     da,
-                    nodata=nodata,
+                    nodata=nodata_val[name],
                     pattern=None,
                     dtype=np.float32,
                     write=_write,
@@ -154,14 +159,18 @@ def _write(path, a, nodata=1.0e20, dtype=np.float32):
                 for prm in ["factor", "angle"]:
                     da = self.dataset[prm]
                     if "layer" in da.coords and "y" in da.coords and "x" in da.coords:
-                        a = da.sel(layer=l).values
-                        if not np.all(a == a[0][0]):
+                        a = da.sel(layer=l)
+                        if not _check_all_equal(a):
                             f.write(
                                 f"open/close {self.rendir.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n"
                             )
                         else:
+                            if np.all(np.isnan(a)):
+                                val = nodata_val[prm]
+                            else:
+                                val = a[~np.isnan()][0]
                             f.write(
-                                f"constant {float(a[0][0]):.5f} {prm}_l{float(l):.0f}\n"
+                                f"constant {float(val):.5f} {prm}_l{float(l):.0f}\n"
                             )
                     else:
                         f.write(

From f6bfd6dd09c57b706e5ab7c1971c313a93c3d480 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Thu, 17 Nov 2022 13:40:12 +0100
Subject: [PATCH 08/11] hfb id changed to hfb6 as reqd by imodwq

---
 imod/wq/hfb.py   | 4 ++--
 imod/wq/model.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 5dd069454..927e835e7 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -16,9 +16,9 @@ class HorizontalFlowBarrier(Package):
         responsible for consistency with model.
     """
 
-    _pkg_id = "hfb"
+    _pkg_id = "hfb6"
 
-    _template = "[hfb]\n" "    hfbfile = {hfbfile}\n\n"
+    _template = "[hfb6]\n" "    hfbfile = {hfbfile}\n\n"
 
     def __init__(
         self,
diff --git a/imod/wq/model.py b/imod/wq/model.py
index 07bde9acb..efdd24e1e 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -607,7 +607,7 @@ def _bas_btn_rch_evt_mal_tvc_sinkssources(self):
 
     def _render_anihfb(self, modelname, directory, nlayer):
         out = ""
-        for key in ("ani", "hfb"):
+        for key in ("ani", "hfb6"):
             pkgkey = self._get_pkgkey(key)
             if pkgkey is not None:
                 out += self[pkgkey]._render(
@@ -810,7 +810,7 @@ def write(
             if (
                 "x" in pkg.dataset.coords
                 and "y" in pkg.dataset.coords
-                or pkg._pkg_id in ("wel", "ani", "hfb")
+                or pkg._pkg_id in ("wel", "ani", "hfb6")
             ):
                 try:
                     pkg.save(directory=directory / pkgname)

From 3d32903c617c69bba31df060f250820fe214efc8 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Tue, 29 Nov 2022 16:16:03 +0100
Subject: [PATCH 09/11] isort

---
 imod/wq/ani.py | 3 ++-
 imod/wq/hfb.py | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 965030860..9d6609594 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -1,8 +1,9 @@
-import numpy as np
 import pathlib
 import re
 import shutil
 
+import numpy as np
+
 import imod
 from imod.wq.pkgbase import Package
 
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 927e835e7..53e007e2f 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -1,5 +1,6 @@
 import pathlib
 import shutil
+
 from imod.wq.pkgbase import Package
 
 

From da14c102af4a1fc03d46ec006bf39e3761ed2ef1 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Wed, 30 Nov 2022 09:48:25 +0100
Subject: [PATCH 10/11] lint cleanup

---
 imod/wq/ani.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 9d6609594..92a2cdaf5 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -156,14 +156,14 @@ def _write(path, da, nodata=1.0e20, dtype=np.float32):
 
         # save anifile with data stored during _render
         with open(directory / self.anifile, "w") as f:
-            for l in range(1, self.nlayer + 1):
+            for lay in range(1, self.nlayer + 1):
                 for prm in ["factor", "angle"]:
                     da = self.dataset[prm]
                     if "layer" in da.coords and "y" in da.coords and "x" in da.coords:
-                        a = da.sel(layer=l)
+                        a = da.sel(layer=lay)
                         if not _check_all_equal(a):
                             f.write(
-                                f"open/close {self.rendir.as_posix()}/{prm}_l{float(l):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(l):.0f}\n"
+                                f"open/close {self.rendir.as_posix()}/{prm}_l{float(lay):.0f}.arr 1.0D0 (FREE) -1 {prm}_l{float(lay):.0f}\n"
                             )
                         else:
                             if np.all(np.isnan(a)):
@@ -171,11 +171,11 @@ def _write(path, da, nodata=1.0e20, dtype=np.float32):
                             else:
                                 val = a[~np.isnan()][0]
                             f.write(
-                                f"constant {float(val):.5f} {prm}_l{float(l):.0f}\n"
+                                f"constant {float(val):.5f} {prm}_l{float(lay):.0f}\n"
                             )
                     else:
                         f.write(
-                            f"constant {float(self.dataset[prm].values):.5f} {prm}_l{float(l):.0f}\n"
+                            f"constant {float(self.dataset[prm].values):.5f} {prm}_l{float(lay):.0f}\n"
                         )
 
     def _pkgcheck(self, ibound=None):

From 21438f9f6a54cd73675926babd4bc4dec3444947 Mon Sep 17 00:00:00 2001
From: Joost Delsman <joost.delsman@deltares.nl>
Date: Fri, 2 Dec 2022 13:33:17 +0100
Subject: [PATCH 11/11] test render

---
 imod/tests/test_wq/test_wq_ani.py | 30 ++++++++++++++++++++++++++++++
 imod/tests/test_wq/test_wq_hfb.py | 17 +++++++++++++++++
 imod/wq/ani.py                    |  4 ++--
 imod/wq/hfb.py                    |  2 +-
 imod/wq/model.py                  | 29 ++++++++++++++++-------------
 5 files changed, 66 insertions(+), 16 deletions(-)
 create mode 100644 imod/tests/test_wq/test_wq_ani.py
 create mode 100644 imod/tests/test_wq/test_wq_hfb.py

diff --git a/imod/tests/test_wq/test_wq_ani.py b/imod/tests/test_wq/test_wq_ani.py
new file mode 100644
index 000000000..8e781efdc
--- /dev/null
+++ b/imod/tests/test_wq/test_wq_ani.py
@@ -0,0 +1,30 @@
+import textwrap
+from pathlib import Path
+
+from imod.wq import HorizontalAnisotropy, HorizontalAnisotropyFile
+
+
+def test_render_ani():
+    ani = HorizontalAnisotropy(factor=1.0, angle=0.0)
+
+    compare = textwrap.dedent(
+        """\
+        [ani]
+            anifile = ani/test.ani
+        """
+    )
+
+    assert ani._render(modelname="test", directory=Path("./ani"), nlayer=1) == compare
+
+
+def test_render_anifile():
+    ani = HorizontalAnisotropyFile(anifile="test.test")
+
+    compare = textwrap.dedent(
+        """\
+        [ani]
+            anifile = ani/test.ani
+        """
+    )
+
+    assert ani._render(modelname="test", directory=Path("./ani"), nlayer=1) == compare
diff --git a/imod/tests/test_wq/test_wq_hfb.py b/imod/tests/test_wq/test_wq_hfb.py
new file mode 100644
index 000000000..92dc6a4e5
--- /dev/null
+++ b/imod/tests/test_wq/test_wq_hfb.py
@@ -0,0 +1,17 @@
+import textwrap
+from pathlib import Path
+
+from imod.wq import HorizontalFlowBarrier
+
+
+def test_render():
+    hfb = HorizontalFlowBarrier(hfbfile=None)
+
+    compare = textwrap.dedent(
+        """\
+        [hfb6]
+            hfbfile = hfb/test.hfb
+        """
+    )
+
+    assert hfb._render(modelname="test", directory=Path("./hfb")) == compare
diff --git a/imod/wq/ani.py b/imod/wq/ani.py
index 92a2cdaf5..29381ee44 100644
--- a/imod/wq/ani.py
+++ b/imod/wq/ani.py
@@ -24,7 +24,7 @@ class HorizontalAnisotropyFile(Package):
 
     _pkg_id = "ani"
 
-    _template = "[ani]\n" "    anifile = {anifile}\n\n"
+    _template = "[ani]\n" "    anifile = {anifile}\n"
 
     def __init__(
         self,
@@ -93,7 +93,7 @@ class HorizontalAnisotropy(Package):
 
     _pkg_id = "ani"
 
-    _template = "[ani]\n" "    anifile = {anifile}\n\n"
+    _template = "[ani]\n" "    anifile = {anifile}\n"
 
     def __init__(
         self,
diff --git a/imod/wq/hfb.py b/imod/wq/hfb.py
index 53e007e2f..f856970f3 100644
--- a/imod/wq/hfb.py
+++ b/imod/wq/hfb.py
@@ -19,7 +19,7 @@ class HorizontalFlowBarrier(Package):
 
     _pkg_id = "hfb6"
 
-    _template = "[hfb6]\n" "    hfbfile = {hfbfile}\n\n"
+    _template = "[hfb6]\n" "    hfbfile = {hfbfile}\n"
 
     def __init__(
         self,
diff --git a/imod/wq/model.py b/imod/wq/model.py
index efdd24e1e..aeff8a94e 100644
--- a/imod/wq/model.py
+++ b/imod/wq/model.py
@@ -605,15 +605,14 @@ def _bas_btn_rch_evt_mal_tvc_sinkssources(self):
 
         return n_extra
 
-    def _render_anihfb(self, modelname, directory, nlayer):
-        out = ""
-        for key in ("ani", "hfb6"):
-            pkgkey = self._get_pkgkey(key)
-            if pkgkey is not None:
-                out += self[pkgkey]._render(
-                    modelname=modelname, directory=directory / pkgkey, nlayer=nlayer
-                )
-        return out
+    def _render_anihfb(self, key, modelname, directory, nlayer):
+        pkgkey = self._get_pkgkey(key)
+        if pkgkey is not None:
+            return self[pkgkey]._render(
+                modelname=modelname, directory=directory / pkgkey, nlayer=nlayer
+            )
+        else:
+            return ""
 
     def render(self, directory, result_dir, writehelp):
         """
@@ -647,11 +646,15 @@ def render(self, directory, result_dir, writehelp):
             )
 
         # ani and hfb packages
-        content.append(
-            self._render_anihfb(
-                modelname=self.modelname, directory=directory, nlayer=nlayer
+        for key in ("ani", "hfb6"):
+            content.append(
+                self._render_anihfb(
+                    key=key,
+                    modelname=self.modelname,
+                    directory=directory,
+                    nlayer=nlayer,
+                )
             )
-        )
 
         # multi-system package group: chd, drn, ghb, riv, wel
         modflowcontent, ssm_content, n_sinkssources = self._render_groups(