TNO Intern

Commit a7e46c4f authored by Florian Knappers's avatar Florian Knappers
Browse files

ptit ruffje

parent af9f5fa5

src/pythermogis/workflow/utc/configuration.py

+1 −3
Original line number Diff line number Diff line
@@ -88,9 +88,7 @@ class UTCConfiguration: 
    exclude_hc_accum: bool = True

    use_bounding_shape: bool = False

    grid_ext: str = ".nc"

    p_values: list[float] = field(

        default_factory=lambda: [10.0, 30.0, 50.0, 90.0]

    )

    p_values: list[float] = field(default_factory=lambda: [10.0, 30.0, 50.0, 90.0])

    # temp_voxet: 'Voxet' = None

    surface_temperature: float = 10.0

    temp_gradient: float = 31.0

src/pythermogis/workflow/utc/stats.py

+4 −2
Original line number Diff line number Diff line 
from scipy.stats import norm

import numpy as np

from scipy.stats import norm





def lognormal_quantile(mean, sd, p):

    z = norm.ppf(1 - p / 100)

    return mean * np.exp(sd * z)





def normal_quantile(mean, sd, p):

    z = norm.ppf(1 - p / 100)

    return mean + sd * z

src/pythermogis/workflow/utc/utc.py

+41 −37
Original line number Diff line number Diff line
@@ -4,16 +4,18 @@ from pathlib import Path 
import numpy as np

import xarray as xr



from pythermogis.transmissivity.calculate_thick_perm_trans import _tpkt_kernel

from pythermogis.workflow.utc.configuration import UTCConfiguration

from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance

from pythermogis.transmissivity.calculate_thick_perm_trans import _tpkt_kernel





def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree:

    root = xr.DataTree()

    base = Path(config.input_data_dir)



    temperature_model = xr.open_dataset(base / "NetherlandsTemperatureModel_extended.nc").transpose("y", "x", "z")



    temperature_model = xr.open_dataset(

        base / "NetherlandsTemperatureModel_extended.nc"

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



    for aquifer in config.aquifers:

        aquifer_ds = xr.Dataset()
@@ -43,14 +45,15 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 


        root[aquifer] = xr.DataTree(dataset=aquifer_ds, name=aquifer)





    for aquifer in root.children:

        print(f"\nProcessing aquifer: {aquifer}")



        start = time.perf_counter()



        aquifer_ds = root[aquifer].to_dataset()

        updated_ds, n_cells = compute_results_for_aquifer(aquifer_ds, aquifer, config, temperature_model)

        updated_ds, n_cells = compute_results_for_aquifer(

            aquifer_ds, aquifer, config, temperature_model

        )

        root[aquifer] = xr.DataTree(dataset=updated_ds, name=aquifer)



        end = time.perf_counter()
@@ -58,8 +61,10 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 


        cells_per_sec = n_cells / elapsed





        print(f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds. {cells_per_sec} cells per second")

        print(

            f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds. "

            f"{cells_per_sec} cells per second"

        )



    if config.results_dir != "":

        zarr_path = f"{config.results_dir}/ROSL_ROSLU_test.zarr"
@@ -69,7 +74,10 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 




def compute_results_for_aquifer(

        aquifer_ds: xr.Dataset, aquifer_name: str, config: UTCConfiguration, temperature_model: xr.Dataset

    aquifer_ds: xr.Dataset,

    aquifer_name: str,

    config: UTCConfiguration,

    temperature_model: xr.Dataset,

):

    output_names = [

        "power",
@@ -137,11 +145,12 @@ def compute_results_for_aquifer( 
        if result.utc > 1000:

            result = result._replace(utc=1000)





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

        return (

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

        )



    # TODO: loop over pvalues and save grids with pvalue noted

    # TODO: Maybe save the results as datasets iso dataarrays, and have p value be a dim?

    # TODO: Maybe save the results as datasets iso dataarrays & have p value be a dim?



    # calc transmissivity(_with_ntg)

    hydro = compute_hydro_properties_mc(aquifer_ds, 90)
@@ -153,14 +162,9 @@ def compute_results_for_aquifer( 
    # calc temperature grid

    mid_depth = -(aquifer_ds["depth"] + 0.5 * aquifer_ds["thickness"])

    temp_xy_aligned = temperature_model["temperature"].interp(

        x=aquifer_ds.x,

        y=aquifer_ds.y,

        method="linear"

    )

    temperature = temp_xy_aligned.interp(

        z=mid_depth,

        method="linear"

        x=aquifer_ds.x, y=aquifer_ds.y, method="linear"

    )

    temperature = temp_xy_aligned.interp(z=mid_depth, method="linear")



    aquifer_ds["temperature"] = temperature
@@ -223,7 +227,6 @@ def compute_hydro_properties_mc( 
    p_value: float,

    n_samples: int = 10_000,

) -> xr.Dataset:



    thickness_mean = aquifer_ds["thickness"]

    thickness_sd = aquifer_ds["thickness_sd"]

    ln_perm_mean = np.log(aquifer_ds["permeability"])
@@ -245,7 +248,9 @@ def compute_hydro_properties_mc( 
    )



    transmissivity_with_ntg = transmissivity_p * aquifer_ds["ntg"] / 1000.0

    transmissivity_with_ntg = transmissivity_with_ntg.where(transmissivity_with_ntg >= 0)

    transmissivity_with_ntg = transmissivity_with_ntg.where(

        transmissivity_with_ntg >= 0

    )



    return xr.Dataset(

        {
@@ -255,4 +260,3 @@ def compute_hydro_properties_mc( 
            "transmissivity_with_ntg": transmissivity_with_ntg,

        }

    )