Creating a new stochastic example

## Other statistical Implementations

The ThermoGIS methodology `calculate_doublet_performance_stochastic` assumes that we know the net-to-gross, porosity and depth of the aquifer perfectly and simulates a range of values only over permeability, thickness (and thus; transmissivity).

This is because transmissivity is usually the most important aquifer property when determining performance and for the ThermoGIS project we make country-wide maps with simulations of 1km x 1km for the Netherlands; this means choices have to be made when selecting what properties to sample over.

However, if conducting a local study you may well want to also incorporate statistical sampling across other reservoir properties, it is relatively easy (and insightful) to use pythermogis to generate samples and make your own stochastic framework.

```python

from pythermogis import calculate_doublet_performance, calculate_pos

import xarray as xr

from matplotlib import pyplot as plt

from pathlib import Path

import numpy as np

from os import path

# output directory to write output to

output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "exceedance_example"

output_data_path.mkdir(parents=True, exist_ok=True)

# generate simulation samples across desired reservoir properties

Nsamples = 1000

thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)

porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples)

ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples)

depth_samples = np.random.normal(loc=2000, scale=100, size=Nsamples)

permeability_samples = np.random.normal(loc=400, scale=100, size=Nsamples)

reservoir_properties = xr.Dataset(

    {

        "thickness": (["sample"], thickness_samples),

        "porosity": (["sample"], porosity_samples),

        "ntg": (["sample"], ntg_samples),

        "depth": (["sample"], depth_samples),

        "permeability": (["sample"], permeability_samples),

    },

    coords={"sample": np.arange(Nsamples)}

)

# For every sample, run a doublet simulation store the output values

simulations = calculate_doublet_performance(reservoir_properties, print_execution_duration=True)

# post process the samples to calculate the percentiles of each variable

percentiles = np.arange(1,99)

results = simulations.quantile(q=percentiles/100, dim="sample")

results = results.assign_coords(quantile=("quantile", percentiles)) # Overwrite the 'quantile' coordinate with integer percentiles

results = results.rename({"quantile": "p_value"}) # Rename dimension to 'p_value' to match nomenclature in the rest of pythermogis

results["pos"] = calculate_pos(results)

# plotting

fig, axes = plt.subplots(3,2, figsize=(10,15))

plt.suptitle("Reservoir Property Distributions")

reservoir_properties.thickness.plot.hist(ax=axes[0,0])

reservoir_properties.porosity.plot.hist(ax=axes[0,1])

reservoir_properties.ntg.plot.hist(ax=axes[1,0])

reservoir_properties.depth.plot.hist(ax=axes[1,1])

reservoir_properties.permeability.plot.hist(ax=axes[2,0])

axes[2,1].set_visible(False)

plt.savefig(output_data_path / "reservoir_distributions.png")

fig, axes = plt.subplots(2,2, figsize=(10,10))

plt.suptitle(f"Simulation results\nNo. of Samples: {Nsamples}")

simulations.plot.scatter(x="permeability", y="transmissivity", ax=axes[0,0])

simulations.plot.scatter(x="permeability", y="utc", ax=axes[0,1])

simulations.plot.scatter(x="permeability", y="power", ax=axes[1,0])

simulations.plot.scatter(x="permeability", y="npv", ax=axes[1,1])

plt.savefig(output_data_path / "simulation_results.png")

fig, axes = plt.subplots(2,2, figsize=(10,10))

plt.suptitle(f"Exceedance probability curves\nNo. of Samples: {Nsamples}")

plot_exceedance(results["permeability"], ax=axes[0,0])

plot_exceedance(results["utc"], ax=axes[0,1])

plot_exceedance(results["power"], ax=axes[1,0])

plot_exceedance(results["npv"], ax=axes[1,1])

axes[1, 1].axhline(results["pos"], ls="--", c="tab:orange", label=f"POS: {results["pos"]:.1f}%")   # add the probability of success

axes[1, 1].axvline(0, ls="--", c="tab:orange")  # add the probability of success

axes[1, 1].legend()

plt.savefig(output_data_path / "exceedance_probabilities.png")

```

 No newline at end of file

      - customized properties: usage/customized_props.md

      - map run and analysis: usage/maprun_analysis.md

      - potrfolio run and analysis: usage/portfoliorun_analysis.md

      - general stochastic methodology: usage/general_stochastic_methodology.md

      #- out maps and locations: usage/pvalues_map_locations.md

  - API Reference:

from pythermogis.postprocessing.pos import *

from doublet_simulation.deterministic_doublet import *

from doublet_simulation.stochastic_doublet import *

from pythermogis.plotting.plot_exceedance import plot_exceedance

 No newline at end of file

from pathlib import Path

import xarray as xr

import matplotlib.pyplot as plt

def plot_exceedance(variable: xr.DataArray, ax = None, outfile: str | Path = None, show = False):

    """

    Provided a variable which has dimensions p_value then plot the exceedance probability of that variable

    """

    if "p_value" not in variable.dims: raise ValueError("Variable does not have dimension 'p_value'")

    if len(variable.dims) > 1: raise ValueError("Variable contains more than one dimension")

    values = variable.values

    p_values = variable.p_value

    survival_probs = 100 - p_values

    if ax is None:

        fig, ax = plt.subplots(1,1,figsize=(8,8))

    ax.plot(values, survival_probs)

    ax.set_xlabel(f"{variable.name}")

    ax.set_ylabel("Exceedance probability (P[X > x])")

    ax.grid(True)

    if outfile is not None:

        plt.savefig(outfile)

    if show:

        plt.show()

 No newline at end of file

    if npv.ndim != 1 or p_values.ndim != 1:

        raise ValueError("Both npv and p_values must be 1D arrays.")

    # Reverse the order of npv and p_values for interpolation

    pos = np.interp(0.0, npv[::-1], p_values[::-1])

    return pos

    if np.all(np.diff(npv) > 0): # if values are strictly increasing

        return 100 - np.interp(0.0, npv, p_values)

    # Reverse the order of npv and p_values for interpolation

    return np.interp(0.0, npv[::-1], p_values[::-1])

def calculate_pos(results:xr.Dataset) -> xr.Dataset:

    """