adding a more complex example to the docs

pixi run mkdocs serve

```

If you wish to build the documentation locally, run:

```bash

pixi run mkdocs serve

```

## What can this package do?

This package allows a user to simulate geothermal doublets providing the following parameters: 

- capex (Capital expenditure) [Million €]

- opex (Operational expenditure) [€/kW]

- utc (Unit Technical Cost [€cent/kWH])

- npv (Net-present-value)

- npv (Net-present-value) [Million €]

- hprod (Discounted Heat Produced)

- cop 

- cophp

The [basic usage](basic_usage.md) page shows how to run a simulation and get outputs, here we want to show how a script with some more advanced reading in and writing out of data would look like with some customized plotting.

If you wish to make your own plots we recommend the user to learn how to use [matplotlib](https://matplotlib.org/) and/or [xarrays](https://docs.xarray.dev/en/stable/).

**Note: pythermogis has been made to enable users to run doublet simulations; it is not our goal to make lots of customised plotting functions. 

The reason is that when it comes to plotting users have a specific idea of how they want their plots to look and it is hard to cater to everyone. It is better for us to give you the tools to run the simulations and let you design your plots 

yourselves.**

Here is an example which:

1. reads in grids of: mean thickness, thickness standard deviation, net-to-gross, porosity, mean permeability, ln(permeability) standard deviation

2. runs a doublet simulations and calculates economics across all the non-nan cells of the grids for p-values of 10%-90% (with increments of 10%)

3. writes the results out to file

4. plots maps of power, capex and npv for the pvalues: 10%, 50%, 90%

5. plots the npv curve for a single location on the grids

The example input data can be found in the `/resoureces/example_data` directory in the repository.

```python

from pythermogis import calculate_doublet_performance

from pygridsio import read_grid, plot_grid

from matplotlib import pyplot as plt

import xarray as xr

import numpy as np

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.dirname(__file__), "resources") / "test_input" / "example_data"

# create a directory to write the output files to

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")

# 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.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"]

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)

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

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

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:

pos = np.interp(0.0, results_loc.npv[::-1], results_loc.p_value[::-1]) # The order of the npv data has to be reversed as np.interp requires values to increase to function properly

# 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].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()

axes[0].set_xlabel("net-present-value [Million €]")

axes[0].set_ylabel("p-value [%]")

# 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

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

```

And the output plots are:

![maps.png](../images/maps.png)

![npv.png](../images/npv.png)

### 🧪 Basic Example

hello!

```python

from src.pythermogis import calculate_doublet_performance

import xarray as xr