Draft: Resolve "Produce a probabilistic play based risk analysis approach poc"

from pathlib import Path

import numpy as np

from matplotlib import pyplot as plt

from pygridsio import read_grid, plot_grid

from tqdm import tqdm

import copy

test_files_path = Path(__file__).parent / "resources"

test_input = test_files_path / "test_input" / "example_data"

# sigmoid functions

temp_midpoint, temp_k = 80, 0.1

perm_midpoint, perm_k = 250, 0.01

thick_midpoint, thick_k = 50, 0.1

ntg_midpoint, ntg_k = 0.5, 10

def smooth_logit(x, midpoint=7.5, k=1.5):

    return 1 / (1 + np.exp(-k * (x - midpoint)))

def test_play_based_risk_single_location():

    temp_range = (0, 200)

    perm_range = (0, 500)

    thick_range = (0, 500)

    ntg_range = (0, 1)

    n_samples = 1000

    temperature = np.random.normal(100, 5, n_samples)

    permeability = np.random.normal(300, 50, n_samples)

    thickness = np.random.normal(200, 20, n_samples)

    ntg = np.random.normal(0.8, 0.1, n_samples)

    # chance of success factors

    temperature_cos = smooth_logit(temperature, midpoint=temp_midpoint, k=temp_k)

    permeability_cos = smooth_logit(permeability, midpoint=perm_midpoint, k=perm_k)

    thickness_cos = smooth_logit(thickness, midpoint=thick_midpoint, k=thick_k)

    ntg_cos = smooth_logit(ntg, midpoint=ntg_midpoint, k=ntg_k)

    # combined chance of success

    combined_cos = temperature_cos * permeability_cos * thickness_cos * ntg_cos

    # plot Alles

    fig, axd = plt.subplot_mosaic(

        [

            ["temp", "temp_sig", "temp_cos", "combined"],

            ["perm", "perm_sig", "perm_cos", "combined"],

            ["thick", "thick_sig", "thick_cos", "combined"],

            ["ntg", "ntg_sig", "ntg_cos", "combined"],

        ],

        figsize=(30, 10),

        constrained_layout=True,

    )

    # temperature

    axd["temp"].set_title("Temperature samples [C]")

    axd["temp"].hist(temperature, bins=20)

    axd["temp_sig"].set_title("Temperature Sigmoid")

    axd["temp_sig"].plot(

        np.linspace(temp_range[0], temp_range[1]),

        smooth_logit(

            np.linspace(temp_range[0], temp_range[1]), midpoint=temp_midpoint, k=temp_k

        ),

    )

    axd["temp_cos"].set_title("Temperature COS")

    axd["temp_cos"].hist(temperature_cos, bins=20)

    # thickness

    axd["thick"].set_title("Thickness samples [m]")

    axd["thick"].hist(thickness, bins=20)

    axd["thick_sig"].set_title("Thickness Sigmoid")

    axd["thick_sig"].plot(

        np.linspace(thick_range[0], thick_range[1]),

        smooth_logit(

            np.linspace(thick_range[0], thick_range[1]),

            midpoint=thick_midpoint,

            k=thick_k,

        ),

    )

    axd["thick_cos"].set_title("Thickness COS")

    axd["thick_cos"].hist(thickness_cos, bins=20)

    # permeability

    axd["perm"].set_title("Permeability samples [mD]")

    axd["perm"].hist(permeability, bins=20)

    axd["perm_sig"].set_title("Permeability Sigmoid")

    axd["perm_sig"].plot(

        np.linspace(perm_range[0], perm_range[1]),

        smooth_logit(

            np.linspace(perm_range[0], perm_range[1]), midpoint=perm_midpoint, k=perm_k

        ),

    )

    axd["perm_cos"].set_title("Permeability COS")

    axd["perm_cos"].hist(permeability_cos, bins=20)

    # ntg

    axd["ntg"].set_title("NTG samples [0-1]")

    axd["ntg"].hist(ntg, bins=20)

    axd["ntg_sig"].set_title("NTG Sigmoid")

    axd["ntg_sig"].plot(

        np.linspace(ntg_range[0], ntg_range[1]),

        smooth_logit(

            np.linspace(ntg_range[0], ntg_range[1]), midpoint=ntg_midpoint, k=ntg_k

        ),

    )

    axd["ntg_cos"].set_title("NTG COS")

    axd["ntg_cos"].hist(ntg_cos, bins=20)

    # combined cos

    axd["combined"].set_title("Combined COS")

    axd["combined"].hist(combined_cos, bins=20)

    plt.show()

def test_cos_grid():

    aquifer = "ROSL_ROSLU"

    ntg = read_grid(test_input / f"{aquifer}__ntg.zmap")

    ln_perm = np.log(read_grid(test_input / f"{aquifer}__perm.zmap"))

    ln_perm_sd = read_grid(test_input / f"{aquifer}__ln_perm_sd.zmap")

    top = read_grid(test_input / f"{aquifer}__top.zmap")

    thick = read_grid(test_input / f"{aquifer}__thick.zmap")

    thick_sd = read_grid(test_input / f"{aquifer}__thick_sd.zmap")

    n_samples = 1000

    cos_grids = []

    cos_mean_grid = copy.deepcopy(ntg)

    cos_sd_grid = copy.deepcopy(ntg)

    for _ in tqdm(range(n_samples)):

        ntg_sample = np.random.normal(ntg.data, 0.1)

        perm_sample = np.exp(np.random.normal(ln_perm.data, ln_perm_sd.data))

        top_sample = np.random.normal(top.data, thick_sd.data)

        temp_sample = 8 + 0.032 * top_sample

        thick_sample = np.random.normal(thick.data, thick_sd.data)

        temp_sig =  smooth_logit(temp_sample, midpoint=temp_midpoint, k=temp_k)

        thick_sig = smooth_logit(thick_sample, midpoint=thick_midpoint, k=thick_k)

        perm_sig = smooth_logit(perm_sample, midpoint=perm_midpoint, k=perm_k)

        ntg_sig = smooth_logit(ntg_sample, midpoint=ntg_midpoint, k=ntg_k)

        cos_sample = temp_sig * thick_sig * perm_sig * ntg_sig

        cos_grids.append(copy.deepcopy(cos_sample))

    # plot the mean and SD COS grid

    cos_mean_grid.data = np.mean(cos_grids, axis=0)

    cos_sd_grid.data = np.std(cos_grids, axis=0)

    fig, ax = plt.subplots(nrows=1, ncols = 2, figsize=(20, 10))

    plot_grid(cos_mean_grid, axes=ax[0], add_netherlands_shapefile=True)

    plot_grid(cos_sd_grid, axes=ax[1], add_netherlands_shapefile=True)

    ax[0].set_title("Mean COS")

    ax[1].set_title("SD COS")

    plt.show()

 No newline at end of file

tests/resources/test_output

/pythermogis/resources/JVM17

src/pythermogis/jvm

**/test_output/**

 No newline at end of file