TNO Intern

Commit 3a61db4f authored by Hen Brett's avatar Hen Brett 🐔
Browse files

adding a chance of success test with maps

parent c261621a
Loading
Loading
Loading
Loading
Loading
+49 −7
Original line number Diff line number Diff line
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"

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():
# sigmoid functions
temp_midpoint, temp_k = 80, 0.1
perm_midpoint, perm_k = 250, 0.01
thick_midpoint, thick_k = 50, 0.01
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)
@@ -115,3 +121,39 @@ def test_play_based_risk_single_location():
    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 = 10

    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)

    plot_grid(cos_mean_grid, add_netherlands_shapefile=True, show=True)
    plot_grid(cos_sd_grid, add_netherlands_shapefile=True, show=True)
 No newline at end of file