Updating the example data and docs

```python

from pythermogis import calculate_doublet_performance

from pygridsio import read_grid, plot_grid # for reading and plotting grids

from matplotlib import pyplot as plt # for more customized plotting

import xarray as xr # to handle the in and output data

from pygridsio import read_grid, plot_grid, resample_xarray_grid, remove_padding_from_grid

from matplotlib import pyplot as plt

import numpy as np

from pathlib import Path # to handle pathing

import xarray as xr

from pathlib import Path

from os import path

# the location of the input data: the data can be found in the resources/example_data directory of the repo

input_data_path = Path("path/to/example_data")

input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data"

# create a directory to write the output files to

output_data_path = Path("path/to/output_data")

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

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

# if set to True then simulation is always run, otherwise pre-calculated results are read (if available)

run_simulation = True

# grids can be in .nc, .asc, .zmap or .tif format: see pygridsio package

reservoir_properties = read_grid(input_data_path / "SLDNA__thick.nc").to_dataset(name="thickness_mean") # thickness mean in [m] of the aquifer

reservoir_properties["thickness_sd"] = read_grid(input_data_path / "SLDNA__thick_sd.nc") # thickness standard deviation in [m] of the aquifer

reservoir_properties["ntg"] = read_grid(input_data_path / "SLDNA__ntg.nc") # the net-to-gross of the aquifer [0-1]

reservoir_properties["porosity"] = read_grid(input_data_path / "SLDNA__phi.nc") / 100 # porosity has to be [0-1]

reservoir_properties["depth"] = read_grid(input_data_path / "SLDNA__top.nc") # depth to the top of the aquifer [m]

reservoir_properties["ln_permeability_mean"] = np.log(read_grid(input_data_path / "SLDNA__k.nc")) # mean permeability has to be in logspace

reservoir_properties["ln_permeability_sd"] = read_grid(input_data_path / "SLDNA__k_ln_sd.nc")

new_cellsize = 5000  # in m, this sets the resolution of the model; to speed up calcualtions you can increase the cellsize

reservoir_properties = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")

reservoir_properties["thickness_sd"] = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)

reservoir_properties["ntg"] = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)

reservoir_properties["porosity"] = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"), new_cellsize=new_cellsize) / 100

reservoir_properties["depth"] = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize)

reservoir_properties["ln_permeability_mean"] = np.log(resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize))

reservoir_properties["ln_permeability_sd"] = resample_xarray_grid(read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"), new_cellsize=new_cellsize)

# simulate

# if doublet simulation has already been run then read in results, or run the simulation and write results out

if (output_data_path / "output_results.nc").exists and not run_simulation:

    results = xr.load_dataset(output_data_path / "output_results.nc")

else:

    results = calculate_doublet_performance(reservoir_properties, p_values=[10,20,30,40,50,60,70,80,90]) # This can take ~5 minutes for the sample dataset

    results = calculate_doublet_performance(reservoir_properties, p_values=[10,20,30,40,50,60,70,80,90])

results.to_netcdf(output_data_path / "output_results.nc") # write doublet simulation results to file

# plotting

# plot results as maps

variables_to_plot = ["power","capex", "npv"]

variables_to_plot = ["transmissivity","power","capex", "npv"]

p_values_to_plot = [10, 50, 90]

fig, axes = plt.subplots(nrows=len(p_values_to_plot), ncols=len(variables_to_plot), figsize=(5*len(variables_to_plot), 5*len(p_values_to_plot)), sharex=True, sharey=True)

fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=len(p_values_to_plot), figsize=(5*len(variables_to_plot), 5*len(p_values_to_plot)), sharex=True, sharey=True)

for j, p_value in enumerate(p_values_to_plot):

    results_p_value = results.sel(p_value=p_value)

    for i, variable in enumerate(variables_to_plot):

        plot_grid(results_p_value[variable], axes=axes[i,j], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile

plt.tight_layout()  # ensure there is enough spacing between subplots

plt.tight_layout()  # ensure there is enough spacing

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

# plot net-present value at a single location as a function of p-value and find the probability of success

x, y = 85e3, 440e3  # define location

x, y = 125e3, 525e3  # define location

results_loc = results.sel(x=x,y=y,method="nearest")

# probability of success is defined as the p-value where npv = 0.0, use interpolation to find pos:

# plot npv versus p-value and a map showing the location of interest

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

results_loc.npv.plot(y="p_value", ax=axes[0])

axes[0].set_title(f"Aquifer: SLDNA\nlocation: [{x:.0f}m, {y:.0f}m]")

axes[0].set_title(f"Aquifer: ROSL_ROSLU\nlocation: [{x:.0f}m, {y:.0f}m]")

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

axes[0].axvline(0.0,ls="--", c="tab:orange")

axes[0].legend()

# plot map

plot_grid(results.sel(p_value=50).npv,axes=axes[1], add_netherlands_shapefile=True)

axes[1].scatter(x,y,marker="x",color="tab:red")

plt.tight_layout()  # ensure there is enough spacing between subplots

plt.tight_layout()  # ensure there is enough spacing

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

```

def validate_input_data(input_data: xr.Dataset):

    """

    Ensure that the input_data Dataset contains the minimum required variables.

    Ensure that the input_data Dataset contains the minimum required variables. Otherwise raise an informative error.

    Parameters

    ----------

    Returns

    -------

    These parameters as a tuple:

    - "power"

    - "heat_pump_power"

    - "capex"

    - "opex"

    - "utc"

    - "npv"

    - "hprod"

    - "cop"

    - "cophp"

    - "pres"

    - "flow_rate"

    - "welld"

    power: float

    heat_pump_power: float

    capex: float

    opex: float

    utc: float

    npv: float

    hprod: float

    cop: float

    cophp: float

    pres: float

    flow_rate: float

    welld: float

    """

    if np.isnan(thickness):