Merge branch 'main' into 'update_doc'

Read the full documenation [here](https://pythermogis-15909e.ci.tno.nl/).

Or, you can deploy the documentation locally using:

If you wish to build the documentation locally, run:

```bash

pixi run mkdocs serve

The code will simulate a Geothermal doublet, utilizing ThermoGIS with DoubletCalc1D as the engine to produce values of:

- power [MW]

- heat pump power [MW]

- capex (Capital expenditure) [Million €/$]

- opex (Operational expenditure) in first year [Million €/$]

- utc (Unit Technical Cost [cts/kWH])

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

- hprod (Discounted Heat Produced) [MWh]

- cop (coefficient of performance of system) [-]

- cophp (coeffecient of performance of heat pump) [-]

- pressure for driving the thermal loop (wells+reservoir) [bar]

- power [Mega Watt Hour]

- heat pump power [Mega Watt Hour]

- capex (Capital expenditure) [Million €]

- opex (Operational expenditure) [€/kW]

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

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

- hprod (Discounted Heat Produced)

- cop 

- cophp

- pressure 

- flow rate [m³/hr]

- reservoir depth [m]

The [basic usage](basic_usage.md) page demonstrates how to run a simulation and retrieve results. This page provides an example of a more advanced workflow using **pyThermoGIS**, including some plotting examples to illustrate the outputs.

!!! info "Plotting"

    pyThermoGIS is designed to enable users to run geothermal doublet simulations. Customized plotting is intentionally limited, as users often have specific preferences for visual presentation. Rather than attempt to accommodate all 

    potential styles, our focus is on providing reliable simulation tools. For creating your own visualizations, we recommend using libraries such as [matplotlib](https://matplotlib.org/) and/or [xarray](https://docs.xarray.dev/en/stable/).

### This example demonstrates the following steps:

1. **Reading input grids** for:

     - Mean thickness  

     - Thickness standard deviation  

     - Net-to-gross ratio  

     - Porosity  

     - Mean permeability  

     - Log(permeability) standard deviation  

2. **Running doublet simulations** and calculating economic indicators (e.g., power production, CAPEX, NPV) across all non-NaN cells in the input grids for a range of p-values (10% to 90% in 10% increments).

3. **Exporting results** to file.

4. **Plotting result maps** for selected p-values (10%, 50%, and 90%) for:

     - Power output  

     - Capital expenditure (CAPEX)  

     - Net Present Value (NPV)

5. **Plotting an NPV curve** at a single specified location on the grid.

Example input data for this workflow is available in the `/resources/example_data` directory of the repository.

```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

import numpy as np

from pathlib import Path # to handle pathing

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

# create a directory to write the output files to

output_data_path = Path("path/to/output_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 between subplots

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 between subplots

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

```

And the output plots produced by the code are:

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

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