getting there

    p = np.asarray(p_values)

    q = 1.0 - p / 100.0

    # P50 analytic shortcut

    if p.size == 1 and p[0] == 50:

        t = thickness_mean

        k = np.exp(ln_perm_mean)

        return np.array([t]), np.array([k]), np.array([t * k])

    # Thickness

    if thickness_sd == 0:

        thickness_p = np.full(p.size, thickness_mean)

    base = Path(config.input_data_dir)

    temperature_model = xr.open_dataset(base / "NetherlandsTemperatureModel_extended.nc")

    temperature_model = temperature_model.transpose("y", "x", "z")

    if coarsen_dataset_by is not None:

        temperature_model = temperature_model.coarsen(x=coarsen_dataset_by, y=coarsen_dataset_by, boundary="trim").mean()

    for aquifer in config.aquifers:

        aquifer_ds = xr.Dataset()

    ]

    def cell_calculation(

            x,

            y,

            unknown_input_value,

            thickness,

            transmissivity,

        if any(np.isnan(v) for v in inputs):

            return tuple(np.nan for _ in output_names)

        if transmissivity_with_ntg < config.kh_cutoff:

            return tuple(0 for _ in output_names)

        doublet_input = DoubletInput(

            unknown_input_value=float(unknown_input_value),

            thickness=float(thickness),

            transmissivity_with_ntg=float(transmissivity_with_ntg),

            ntg=float(ntg),

            depth=float(depth),

            porosity=float(porosity),

            porosity=float(porosity / 100),

            temperature=float(temperature),

        )

        try:

            result = calculate_doublet_performance(config, doublet_input)

        except ValueError:

            print(

                f"Failure at x={x:.1f}, y={y:.1f} | "

                f"thick={thickness:.3f}, temp={temperature:.2f}, "

                f"trans={transmissivity:.3f}, trans_ntg={transmissivity_with_ntg:.3e}"

                f"ntg={ntg}, depth={depth:.3f}, porosity={porosity:.3f}"

            )

            return tuple(np.nan for _ in output_names)

        if result.utc > 4000:

            print(

                f"Hmmmmm at x={x:.1f}, y={y:.1f} | "

                f"thick={thickness:.3f}, temp={temperature:.2f}, "

                f"trans={transmissivity:.3f}, trans_ntg={transmissivity_with_ntg:.3e}, "

                f"ntg={ntg}, depth={depth:.3f}, porosity={porosity:.3f}, "

                f"utc:{result.utc}"

            )

        return tuple(result) if result is not None else tuple(np.nan for _ in output_names)

    vectorized_calc = np.vectorize(

    # calc transmissivity(_with_ntg)

    p_values = xr.DataArray([50], dims="p")

    hydro = compute_hydro_properties_mc(aquifer_ds, p_values)

    transmissivity = hydro["transmissivity"]

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"]

    transmissivity = hydro["transmissivity"].squeeze("p", drop=True)

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"].squeeze("p", drop=True)

    aquifer_ds["transmissivity"] = transmissivity

    aquifer_ds["transmissivity_with_ntg"] = transmissivity_with_ntg

        method="linear"

    )

    aquifer_ds["temperature"] = temperature

    y_coord, x_coord = xr.broadcast(aquifer_ds["y"], aquifer_ds["x"])

    # use hc_accum mask

    hc_mask = aquifer_ds["hc_accum"] != 0

    vars_to_mask = [

        "thickness",

        "ntg",

        "porosity",

        "transmissivity",

        "transmissivity_with_ntg",

        "temperature",

    ]

    for v in vars_to_mask:

        aquifer_ds[v] = aquifer_ds[v].where(hc_mask)

    aquifer_ds["temperature"] = temperature

    results = xr.apply_ufunc(

        vectorized_calc,

        -9999,

        x_coord,

        y_coord,

        1.0E30,

        aquifer_ds["thickness"],

        transmissivity,

        transmissivity_with_ntg,

        aquifer_ds["transmissivity"],

        aquifer_ds["transmissivity_with_ntg"],

        aquifer_ds["ntg"],

        aquifer_ds["depth"],

        aquifer_ds["porosity"],

        temperature,

        aquifer_ds["temperature"],

        output_core_dims=[[] for _ in output_names],

        vectorize=False,

        dask="parallelized",

        }

    )

if __name__ == "__main__":

    config = UTCConfiguration(

        input_data_dir="C:/Users/knappersfy/work/thermogis/pythermogis/workdir/InputData",

        results_dir=".",

        aquifers=["ROSL_ROSLU"],

        segment_length=1,

    )

    run_utc_workflow(config, 2)

        f"{n_sims} simulations took: {time_elapsed:.1f} seconds\n"

        f"{n_sims/time_elapsed:.1f} simulations per second"

    )

def test_why_does_this_fail():

    props = UTCConfiguration(

        segment_length=1,

        dh_return_temp=40

    )

    # Create a minimal valid DoubletInput

    input_data = DoubletInput(

        unknown_input_value=-9999,

        thickness=5.259,

        transmissivity=1239.250,

        transmissivity_with_ntg=1.032,

        ntg=0.8331,

        depth=852.969,

        porosity=0.22168,

        temperature=39.30,

    )

    # Act

    results = calculate_doublet_performance(props, input_data)

    print(results)

import matplotlib.pyplot as plt

from pythermogis.workflow.utc.configuration import UTCConfiguration

from pythermogis.workflow.utc.utc import run_utc_workflow

    )

    result = run_utc_workflow(config)

    result.ROSL_ROSLU.utc.plot()

    plt.show()

    print(result)

 No newline at end of file