pilot point parametrization unsupported for disv grid with pstFrom?

[sweco-seviwo]

I have been trying to get pilot points set up with pstFrom to no avail. It seems like pyemus write_list_tpl in the PstFrom class isnt currently supporting disv-style grids and requires cols and rows basically. Is it possible to implement disv support?

import rasterio
import numpy as np
import importlib
import helpers
importlib.reload(helpers)
from helpers import add_rch_pars

# Physical bounds on recharge field: 10–300 mm/yr, converted to m/s
seconds_per_year = 365.25 * 24 * 3600.0
rch_rate_lb = 10.0 / 1000.0 / seconds_per_year
rch_rate_ub = 300.0 / 1000.0 / seconds_per_year

rch_factor_lb = 0.3
rch_factor_ub = 1.5

# Sample recharge zone raster at DISV cell centroids
raster_path = r"E:\OLP\Loddby_Klinga\VOR\GIS\Recharge_Classes_99_1630_DS_260217.asc"
xc = gwf.disv.cell2d.array["xc"]
yc = gwf.disv.cell2d.array["yc"]

with rasterio.open(raster_path) as src:
    nodata = src.nodata
    zone_vals = np.array([v[0] for v in src.sample(zip(xc, yc))], dtype=float)

if nodata is not None:
    zone_vals[zone_vals == nodata] = 0
zone_vals = zone_vals.astype(int)

# Mask: active (1) and pass-through (-1) cells get raster class, inactive (0) get 0
rch_zone_arr = np.where(gwf.disv.idomain.array[0, :] != 0, zone_vals, 0)
print(f"Recharge zones: {np.unique(rch_zone_arr)}")
print(f"  0 = inactive, others = recharge zone classes from raster")
print(f"  Active recharge cells: {(rch_zone_arr > 0).sum()} / {rch_zone_arr.size}")

# Build zone array shaped (nlay, ncpl) — recharge in layers 1 and 4
nlay = gwf.disv.idomain.array.shape[0]
ncpl = gwf.disv.idomain.array.shape[1]
rch_zone_arr_2d = np.zeros((nlay, ncpl), dtype=int)
rch_zone_arr_2d[0, :] = rch_zone_arr  # layer 1
rch_zone_arr_2d[3, :] = rch_zone_arr  # layer 4

# Collect and sort SP files
rch_files = sorted([
    f for f in os.listdir(template_ws)
    if f.lower().endswith(".txt")
    and ("rch_stress_period_data" in f.lower()
         or ("rch" in f.lower() and "stress_period" in f.lower()))
])

if not rch_files:
    print("No recharge stress period data files found")
else:
    print(f"\nFound {len(rch_files)} recharge file(s) — parameterizing first SP only")
    for f in rch_files:
        print(f"  - {f}")
    first_file = rch_files[0]

    helpers.add_rch_pars(
        pf=pf,
        files=first_file,
        zone_array=rch_zone_arr_2d,
        pp_gs=pp_gs,
        lb=rch_factor_lb,
        ub=rch_factor_ub,
        ulb=rch_rate_lb,
        uub=rch_rate_ub,
        pp_space=5,
    )

gives me

Recharge zones: [ 0  2  3  4  8 15]
  0 = inactive, others = recharge zone classes from raster
  Active recharge cells: 28051 / 28700

Found 3 recharge file(s) — parameterizing first SP only
  - vor00_260216.rch_stress_period_data_1.txt
  - vor00_260216.rch_stress_period_data_2.txt
  - vor00_260216.rch_stress_period_data_3.txt
  Recharge (pilot points + zones): vor00_260216.rch_stress_period_data_1.txt
   Error adding recharge parameters: get_tpl_or_ins_df() error: unrecognized 'typ', if not 'obs', should be 'constant','zone', or 'grid', not 'pilotpoints'

helper function:

def add_rch_pars(
    pf,
    files,
    zone_array,
    pp_gs,
    lb=0.2, ub=5.0,
    ulb=0.0, uub=1e-3,
    pp_space=5,
):
    base = "rchrch"
    print(f"  Recharge (pilot points + zones): {files}")

    try:
        # --- Stage 1: Zone constants ---
        df_zn = pf.add_parameters(
            files,
            zone_array=zone_array,
            par_type="zone",
            par_name_base=base + "_zn",
            pargp=base + "_zn",
            index_cols=[0, 1],
            use_cols=[2],
            lower_bound=lb, upper_bound=ub,
            ult_lbound=ulb, ult_ubound=uub,
        )

        # --- Stage 2: Pilot points ---
        df_pp = pf.add_parameters(
            files,
            zone_array=zone_array,
            par_type="pilotpoints",
            geostruct=pp_gs,
            par_name_base=base + "_pp",
            pargp=base + "_pp",
            index_cols=[0, 1],
            use_cols=[2],
            lower_bound=0.5, upper_bound=2.0,
            ult_lbound=ulb, ult_ubound=uub,
            pp_space=pp_space,
        )

        n_zones = len(np.unique(zone_array[zone_array > 0]))
        print(f"    ✓ Stage 1: {n_zones} zone constants  ({lb}–{ub})")
        print(f"    ✓ Stage 2: {len(df_pp)} pilot points (0.5–2.0× multiplier)")
        print(f"    → Final RCH = zone_constant × pp_multiplier")

        return df_zn, df_pp

    except Exception as e:
        print(f"Error adding recharge parameters: {e}")
        return None

[p-ortega]

EDIT: I see that your issue is with recharge, dunno about that tbh .. keeping the rest if it helps!

@briochh can chime in if I'm wrong but I think the current recipe is as follows:

        df_pp = pf.add_parameters(
            files,
            zone_array=zone_array,
            par_type="pilotpoints",
            geostruct=pp_gs,
            par_name_base=base + "_pp",
            pargp=base + "_pp",
            index_cols=[0, 1],
            use_cols=[2],
            lower_bound=0.5, upper_bound=2.0,
            ult_lbound=ulb, ult_ubound=uub,
            pp_options={"prep_hyperpars":False,
                                "pp_space":25} #in this case is meters
        )

pp_space for disv would be the distance between ppoints .. It works for me

[sweco-seviwo]

well, the problem doesn't seem to be the spacing. but that its a disv-grid.

Error adding recharge parameters: get_tpl_or_ins_df() error: unrecognized 'typ', if not 'obs', should be 'constant','zone', or 'grid', not 'pilotpoints'.

it doesnt like index_cols and use_cols in the recepie for disv-grids.

[jtwhite79]

I think if you switch recharge by arrays (rcha) you can use pilotpoints. Pstfrom doesn't yet support pilotpoints for list style parameterization

[d-strom]

Ahha! I wasnt aware that you could! Is that available through flopy?

[briochh]

@d-strom, yep flopy supports rcha (and evta) just fine.