TNO Intern

Commit 6a954780 authored by Florian Knappers's avatar Florian Knappers
Browse files

Resolve "fix examples after removing utc_cutoff due to SDE classes update"

parent 43904527
Loading
Loading
Loading
Loading
+4 −1
Original line number Diff line number Diff line
# Change log

## v3.0.2 (06-7-2026)
## v3.0.3 (07-7-2026)
- removed utc_cutoff from the UTCSettings class, as it is no longer used in the ThermoGIS Java code.

## v3.0.2 (07-7-2026)
- Fix installation issues for macOS users with a silicon chip.
- Add pydantic to project dependencies

+0 −1
Original line number Diff line number Diff line
%% 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,
    utc_cutoff=10,
)
```

%% 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)
+0 −1
Original line number Diff line number Diff line
%% 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,
    utc_cutoff=6.0,
)
```

%% 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")
```
+0 −1
Original line number Diff line number Diff line
%% Cell type:markdown id:e30949364f0b2cad tags:

This example 11 builds on example 6
it deploys an alternative energy conversion scenario, a chilling model (for data centres)
The power conversion is marked by an conversion efficiency of 0.4,  and a return temperature from the network of 40C.

%% Cell type:code id:e7680770ddf31b41 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") / "example11"
output_data_path.mkdir(parents=True, exist_ok=True)
```

%% Cell type:code id:23668b960e88cbc3 tags:

``` 
settings = UTCSettings(
    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,
    dh_return_temp=60,
    load_hours=8760,
    use_kestin_viscosity=True,
    utc_cutoff=10,
)
```

%% Cell type:code id:3857ce2abeb1110 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:6a91eb33cb456d18 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:dac50bb4b203aee8 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:c5640da280d2e009 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:bb9c9ec01a427cf0 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")
```
+0 −1
Original line number Diff line number Diff line
@@ -60,7 +60,6 @@
    "    dh_return_temp=60,\n",
    "    load_hours=8760,\n",
    "    use_kestin_viscosity=True,\n",
    "    utc_cutoff=25,\n",
    ")"
   ],
   "outputs": [],
Loading