Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Wq add anihfb huite #691

Open
wants to merge 3 commits into
base: wq_add_anihfb
Choose a base branch
from
Open

Wq add anihfb huite #691

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
11199ad
Simplify and refactor WQ ANI
Huite Feb 15, 2023
1a30c41
Split WQ MODFLOW render packages into required and optional
Huite Feb 15, 2023
2f11b18
Use xugrid in hfb, format
Huite Feb 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
188 changes: 53 additions & 135 deletions imod/wq/ani.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,79 +1,21 @@
import pathlib
import re
import shutil

import numpy as np

import imod
from imod.wq.pkgbase import Package

FLOAT_FORMAT = "%.18G"

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.
def _write_arr(path, da, nodata=1.0e20, *_):
"""

_pkg_id = "ani"

_template = "[ani]\n" " anifile = {anifile}\n"

def __init__(
self,
anifile,
):
super().__init__()
self["anifile"] = anifile

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
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 / self.anifile, "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"{self.rendir.as_posix()}/{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
Write an ".ARR" file: a plaintext array file which can be understood by MODFLOW.
"""
dx, xmin, xmax, _, 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)
np.savetxt(path, a, fmt=FLOAT_FORMAT, delimiter=" ", footer=footer, comments="")
return


class HorizontalAnisotropy(Package):
Expand All @@ -84,11 +26,11 @@ class HorizontalAnisotropy(Package):
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).
The anisotropic factor perpendicular to the main principal axis which
is the axis of highest permeability. Factor between 0.0 and 1.0 (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).
The angle along the main principal axis (highest permeability) measured
in degrees from north (0), east (90), south (180) and west (270).
"""

_pkg_id = "ani"
Expand All @@ -104,79 +46,55 @@ def __init__(
self["factor"] = factor
self["angle"] = angle

def _render(self, modelname, directory, nlayer):
def _render(self, directory, *args, **kwargs):
# 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"}

d = {"anifile": (directory / "horizontal_anistropy.ani").as_posix()}
return self._template.format(**d)

def save(self, directory):
"""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

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=""
)
def _render_anifile(self, directory):
"""
Unfortunately, the iMOD-WQ anisotropy works through an .ANI file, which
then refers to .ARR files rather than a single indirection.

So the runfile section points to the .ANI file, which in turn points to
the .ARR files, or contains constants.
"""
content = []
for varname, nodata in zip(("factor", "angle"), (1.0, 0.0)):
variable = self.dataset[varname]
if "x" in da.coords and "y" in da.coords:
for layer, da in variable.groupby("layer"):
path = (directory / f"{varname}_l{layer}.arr").as_posix()
content.append(
f"open/close {path} 1.0D0 (FREE) -1 {variable}_l{layer}"
)
else:
# write single value to ani file
pass
for layer, da in variable.groupby("layer"):
value = da.item()
if np.isnan(value):
value = nodata
content.append(f"constant {value:FLOAT_FORMAT} {varname}_l{layer}")
return "\n".join(content)

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")
def save(self, directory):
ani_content = self._render_anifile(directory)
path = (directory / "horizontal_anistropy.ani").as_posix()
with open(path, "w") as f:
f.write(ani_content)

# Write .ARR files
for varname, nodata in zip(("factor", "angle"), (1.0, 0.0)):
da = self.dataset[varname]
if "x" in da.coords and "y" in da.coords:
path = directory / f"{varname}.arr"
imod.array_io.writing._save(
path,
da,
nodata=nodata_val[name],
nodata=nodata,
pattern=None,
dtype=np.float32,
write=_write,
write=_write_arr,
)

# save anifile with data stored during _render
with open(directory / self.anifile, "w") as f:
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=lay)
if not _check_all_equal(a):
f.write(
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)):
val = nodata_val[prm]
else:
val = a[~np.isnan()][0]
f.write(
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(lay):.0f}\n"
)

def _pkgcheck(self, ibound=None):
pass
return
98 changes: 78 additions & 20 deletions imod/wq/hfb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,48 +1,106 @@
import pathlib
import shutil
import numpy as np
import xugrid as xu

from imod.wq.pkgbase import Package


def create_hfb_array(notnull, resistance, layer, ibound, grid):
nlayer = ibound["ibound_layer"].size
idomain2d = ibound.values.reshape((nlayer, -1))
no_layer_dim = notnull.ndim == 1
edge_faces = grid.edge_face_connectivity

# Fill in the indices
# For every edge, find the connected faces.
if no_layer_dim:
edge = np.argwhere(notnull).transpose()[0]
layer = layer - 1
cell2d = edge_faces[edge]
valid = ((cell2d != grid.fill_value) & (idomain2d[layer, cell2d] > 0)).all(
axis=1
)
else:
layer, edge = np.argwhere(notnull).transpose()
layer2d = np.repeat(layer, 2).reshape((-1, 2))
cell2d = edge_faces[edge]
valid = ((cell2d != grid.fill_value) & (idomain2d[layer2d, cell2d] > 0)).all(
axis=1
)
layer = layer[valid]

# Skip the exterior edges (marked with a fill value).
cell2d = cell2d[valid]
c = resistance[notnull][valid]

# Define the numpy structured array dtype
field_spec = [
("layer_1", np.int32),
("row_1", np.int32),
("column_1", np.int32),
("layer_2", np.int32),
("row_2", np.int32),
("column_2", np.int32),
("resistance", np.float64),
]
dtype = np.dtype(field_spec)
shape = (ibound["y"].size, ibound["x"].size)
row_1, column_1 = np.unravel_index(cell2d[:, 0], shape)
row_2, column_2 = np.unravel_index(cell2d[:, 1], shape)
# Set the indices
recarr = np.empty(len(cell2d), dtype=dtype)
recarr["layer_1"] = layer + 1
recarr["row_1"] = row_1 + 1
recarr["column_1"] = column_1 + 1
recarr["row_2"] = row_2 + 1
recarr["column_2"] = column_2 + 1
recarr["resistance"] = c
return recarr


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.
resistance: xu.UgridDataArray
ibound: xr.DataArray
"""

_pkg_id = "hfb6"

_template = "[hfb6]\n" " hfbfile = {hfbfile}\n"
_template = "[hfb6]\n hfbfile = {hfbfile}\n"

def __init__(
self,
hfbfile,
resistance,
ibound,
):
super().__init__()
self["hfbfile"] = hfbfile

def _render(self, modelname, directory, *args, **kwargs):
self.hfbfile = f"{modelname}.hfb"
d = {"hfbfile": f"{directory.as_posix()}/{modelname}.hfb"}
self.dataset = xu.UgridDataset()
self["resistance"] = resistance
self["ibound"] = ibound.rename({"layer": "ibound_layer"})

def _render(self, directory, *args, **kwargs):
d = {"hfbfile": (directory / "horizontal_flow_barrier.hfb").as_posix()}
return self._template.format(**d)

def save(self, directory):
"""Overloads save function.
Saves hfbfile to directory
assumes _render() to have run previously"""
directory.mkdir(exist_ok=True) # otherwise handled by idf.save
path = (directory / "horizontal_flow_barrier.hfb").as_posix()

path_hfb = pathlib.Path(str(self["hfbfile"].values))
hfb_array = create_hfb_array(
notnull=self["resistance"].notnull().values,
resistance=self["resistance"].values,
layer=self["resistance"].coords["layer"].values,
ibound=self["ibound"],
grid=self.dataset.ugrid.grid,
)
nhfbnp = len(hfb_array)
header = f"0 0 {nhfbnp} 0"
fmt = ("%i",) * 5 + ("%.18G",)

shutil.copyfile(path_hfb, directory / self.hfbfile)
with open(path, "w") as f:
np.savetxt(fname=f, X=hfb_array, fmt=fmt, header=header)

def _pkgcheck(self, ibound=None):
pass
Loading