Skip to content

bug: UcnFile (load concetrations from .con file) does load in Flopy #2559

@CooperGarth

Description

@CooperGarth

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:

  1. Run transport model in Modflow 6 and save .con file for concentrations
  2. try Flopy.util.UcnFile to import concentration file (.con)
  3. 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")

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions