-
Notifications
You must be signed in to change notification settings - Fork 349
Description
Describe the bug
When I try to load in a concentration file from a transport simulation created in Modflow 6, I can't load it into Flopy using Flopy.util.UcnFile.
To Reproduce
Steps to reproduce the behavior:
- Run transport model in Modflow 6 and save .con file for concentrations
- try Flopy.util.UcnFile to import concentration file (.con)
- Works for heads file
Expected behavior
Python Script gets stuck on loading concentrations. Works ok for heads file
Additional context
Code here:
mport os
import pandas as pd
import numpy as np
from flopy.utils import HeadFile, UcnFile
from scipy.spatial import cKDTree
from tqdm import tqdm
=== USER INPUT ===
base_path = r"D:\Groundwater Models\580\PR_27_a\PR_27_a"
output_file_name = "PR_27_sal_tr.con" # or .hds
===================
os.chdir(base_path)
File paths
bore_file = os.path.join(base_path, "Bore_locs.csv")
node_file = os.path.join(base_path, "nodes_xyz.csv")
output_file = os.path.join(base_path, output_file_name)
base_name = os.path.splitext(output_file_name)[0]
output_csv = os.path.join(base_path, f"{base_name}_extracted.csv")
pivot_csv = os.path.join(base_path, f"{base_name}_pivot.csv")
Constants
total_nodes = 707520
layers = 8
nodes_per_layer = total_nodes // layers
Load data
bore_locs = pd.read_csv(bore_file)
nodes_xyz = pd.read_csv(node_file)
nodes_xyz["Node"] = nodes_xyz["Node"].astype(int)
bore_order = bore_locs["Name"].tolist()
Precompute node lookup
node_lookup = {}
for _, row in bore_locs.iterrows():
name, x, y, layer = row["Name"], row["X"], row["Y"], int(row["Layer"])
if not (1 <= layer <= layers):
continue
start = (layer - 1) * nodes_per_layer
end = layer * nodes_per_layer
layer_nodes = nodes_xyz.iloc[start:end].copy()
tree = cKDTree(layer_nodes[["X", "Y"]].values)
dist, idx = tree.query([x, y])
cell_node = int(layer_nodes.iloc[idx]["Node"])
node_index = cell_node - 1
node_lookup[name] = (layer, cell_node, node_index)
Detect file type
ext = output_file_name.lower()
if ext.endswith(".ucn") or ext.endswith(".con"):
datafile = UcnFile(output_file,precision='single')
elif ext.endswith(".hds"):
datafile = HeadFile(output_file)
else:
raise ValueError("Unsupported file type")