diff --git a/.gitignore b/.gitignore index e385dceba51eefdcdf61a82ce2f8a640be3f583c..c41961ec8226587f8a0546f5ba2b031cb8431544 100644 --- a/.gitignore +++ b/.gitignore @@ -10,5 +10,5 @@ dist tests/resources/test_output -/pythermogis/resources/JVM17 +src/pythermogis/jvm **/test_output/** \ No newline at end of file diff --git a/tests/test_play_based_risk.py b/tests/test_play_based_risk.py new file mode 100644 index 0000000000000000000000000000000000000000..cd90a4d59af1952dc2332021def48db2c2ef018c --- /dev/null +++ b/tests/test_play_based_risk.py @@ -0,0 +1,166 @@ +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