Merge branch...

%% Cell type:markdown id:afb3177a698c883b tags:

This example shows how to set up the UTC properties and use them in the doublet performance calculation.

The UTC properties are used to calculate the performance of a doublet system with a heat pump.

The UTC properties modified here activate the heatpump model, the goal and return temperature of the heat network

 and the use of the Kestin viscosity model.

It also modifies the cost parameters of the heat pump, which are used to calculate the UTC cost of the system.

%% Cell type:code id:initial_id tags:

``` python

from pythermogis.aquifer import Aquifer

from pythermogis.doublet import ThermoGISDoublet

from pythermogis.settings import UTCSettings

```

%% Cell type:code id:d78e3e6f8909074d tags:

``` python

aquifer = Aquifer(

    thickness=300,

    ntg=0.5,

    porosity=0.2,

    depth=2000,

    permeability=300,

)

```

%% Cell type:code id:59bb4943298aecec tags:

``` python

settings = UTCSettings(

    use_heatpump=True,

    goal_temperature=100,

    dh_return_temp=60,

    use_kestin_viscosity=True,

    heat_pump_capex=600,

    heat_pump_opex=60,

    heat_pump_alternative_heating_price=2.8,

    capex_constant=0.5,

    capex_variable=1100,

    heat_exchanger_efficiency=0.4,

    heat_exchanger_parasitic=0.1,

    min_production_temp=70,

    max_cooling_temperature_range=40,

    load_hours=8760,

)

```

%% Cell type:code id:8a75abccf9d7e13e tags:

``` python

doublet = ThermoGISDoublet(aquifer=aquifer, settings=settings)

results = doublet.simulate().to_dataset()

```

%% Cell type:code id:a01bab0bd12b4a7b tags:

``` python

print(results.data_vars.keys())

```

%% Output

    KeysView(Data variables:

        power                    float64 8B 7.527

        heat_pump_power          float64 8B 11.26

        capex                    float64 8B 21.23

        opex                     float64 8B -6.208

        utc                      float64 8B 14.31

        npv                      float64 8B -16.19

        hprod                    float64 8B 1.391e+06

        cop                      float64 8B 1.797

        cophp                    float64 8B 3.254

        pres                     float64 8B 60.0

        flow_rate                float64 8B 439.2

        welld                    float64 8B 1.25e+03

        inj_temp                 float64 8B 36.65

        prd_temp                 float64 8B 75.92

        transmissivity_with_ntg  float64 8B 45.0

        thickness                int64 8B 300

        transmissivity           int64 8B 90000

        permeability             int64 8B 300

        temperature              float64 8B 76.65)

%% Cell type:markdown id:6878f1eee3b9dc14 tags:

this example builds on example 6

it deploys an alternative energy conversion scenario, a heat pump model

marked by a potential increase of temperature (Goaltemperaure =90C) supplied the heat network, and a return temperature from the network of 50C.

%% Cell type:code id:77da872b2557c249 tags:

``` 

from pathlib import Path

import numpy as np

import xarray as xr

from matplotlib import pyplot as plt

from pygridsio import plot_grid, read_grid, resample_grid

from pythermogis.aquifer import StochasticAquifer

from pythermogis.doublet import ThermoGISDoublet

from pythermogis.postprocess import calculate_pos

from pythermogis.settings import UTCSettings

# the location of the input data: the data can be found in the

# resources/example_data directory of the repo

input_data_path = Path("test_input") / "example_data"

# create a directory to write the output files to

output_data_path = Path("test_output") / "example6"

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

```

%% Cell type:code id:d6fb885005380965 tags:

``` 

settings = UTCSettings(

    use_heatpump=True,

    max_cooling_temperature_range=40,

    minimum_injection_temperature=15,

    heat_pump_capex=600,

    heat_pump_opex=60,

    heat_pump_alternative_heating_price=2.8,

    goal_temperature=90,

    dh_return_temp=50,

    use_kestin_viscosity=True,

)

```

%% Cell type:code id:e4df70f914b6bf1b tags:

``` 

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

new_cellsize = 5000  # in m, this sets the resolution of the model;

# to speed up calcualtions you can increase the cellsize

thickness_mean = resample_grid(

    read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize

)

thickness_sd = resample_grid(

    read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize

)

ntg = resample_grid(

    read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize

)

porosity = (

    resample_grid(

        read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"), new_cellsize=new_cellsize

    )

    / 100

)

depth = resample_grid(

    read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize

)

ln_permeability_mean = np.log(

    resample_grid(

        read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize

    )

)

ln_permeability_sd = resample_grid(

    read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"),

    new_cellsize=new_cellsize,

)

aquifer = StochasticAquifer(

    thickness_mean=thickness_mean,

    thickness_sd=thickness_sd,

    ntg=ntg,

    porosity=porosity,

    depth=depth,

    ln_permeability_mean=ln_permeability_mean,

    ln_permeability_sd=ln_permeability_sd,

)

```

%% Cell type:code id:35aff3b6bdc55248 tags:

``` 

# output inputs

variables_to_plot = ["depth", "thickness_mean", "thickness_sd"]

fig, axes = plt.subplots(

    nrows=len(variables_to_plot),

    ncols=1,

    figsize=(7, 5 * len(variables_to_plot)),

    sharex=True,

    sharey=True,

)

for i, variable in enumerate(variables_to_plot):

    plot_grid(

        getattr(aquifer, variable), axes=axes[i], 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 / "input_maps_1.png")

variables_to_plot = ["porosity", "ln_permeability_mean", "ln_permeability_sd"]

fig, axes = plt.subplots(

    nrows=len(variables_to_plot),

    ncols=1,

    figsize=(7, 5 * len(variables_to_plot)),

    sharex=True,

    sharey=True,

)

for i, variable in enumerate(variables_to_plot):

    plot_grid(

        getattr(aquifer, variable), axes=axes[i], 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 / "input_maps_2.png")

```

%% Cell type:code id:a183a183d1e7d5aa tags:

``` 

# if set to True then simulation is always run,

# otherwise pre-calculated results are read (if available)

run_simulation = True

# 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:

    doublet = ThermoGISDoublet(aquifer=aquifer, settings=settings)

    results = doublet.simulate(

        p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90],

        chunk_size=100,

    ).to_dataset()

    results.to_netcdf(

        output_data_path / "output_results.nc"

    )  # write doublet simulation results to file

print("ready with simulation")

```

%% Cell type:code id:132e54899b47cae7 tags:

``` 

# plot results as maps

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

p_values_to_plot = [10, 50, 90]

fig, axes = plt.subplots(

    nrows=len(variables_to_plot),

    ncols=len(p_values_to_plot),

    figsize=(3 * 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):

        if variable == "utc":

            results_p_value[variable] = results_p_value[variable].clip(0, 25)

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

```

%% Cell type:code id:8c189af10ff9fb42 tags:

``` 

# plot net-present value at a single location as a function of

# p-value and find the probability of success

results["pos"] = calculate_pos(

    results

)  # calculate probability of success across whole aquifer

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

results_loc = results.sel(

    x=x, y=y, method="nearest"

)  # select only the location of interest

pos = (

    results_loc.pos.data

)  # get probability of success at location of interest as a single value

# 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: 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()

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

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

plot_grid(results.pos, 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 / "pos.png")

```