TNO Intern

Commit 3db530ec authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Merge branch '45-change-the-example-data-to-be-the-rotliegend-fix-docstrings' into 'main' 

Resolve "Change the example data to be the rotliegend, fix docstrings"

Closes #45

See merge request AGS/pythermogis!54
parents 5ef5914d 2cf5ea88

docs/images/maps.png

+70.9 KiB (526 KiB)
Loading image diff...

docs/images/npv.png

−962 B (91.2 KiB)
Loading image diff...

docs/usage/advanced_examples.md

+25 −23
Original line number Diff line number Diff line
@@ -30,53 +30,55 @@ Example input data for this workflow is available in the `/resources/example_dat 


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



# create a directory to write the output files to

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

output_data_path = Path("path/to/example_data/output")

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

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

# 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:
@@ -85,7 +87,7 @@ pos = np.interp(0.0, results_loc.npv[::-1], results_loc.p_value[::-1]) # The ord 
# 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()
@@ -95,7 +97,7 @@ 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.tight_layout()  # ensure there is enough spacing

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

```

resources/example_data/example_data.zip

+728 KiB (877 KiB)

File changed.

No diff preview for this file type.

View original file
View changed file

src/pythermogis/thermogis_classes/doublet.py

+137 −66
Original line number Diff line number Diff line
@@ -15,22 +15,56 @@ def calculate_doublet_performance(input_data: xr.Dataset, 
                                  p_values: List[float] = [50.0],

                                  ) -> xr.Dataset:

    """

    Perform a ThermoGIS Doublet performance simulation. This will occur across all dimensions of the input_data (ie. input data can have a single value for each required variable, or it can be 1Dimensional or a 2Dimensional grid)



    :param input_data:

        A xr.Dataset (input_data) with the required variables: "thickness_mean", "thickness_sd", "porosity", "ntg", "depth", "ln_permeability_mean", "ln_permeability_sd", Performance will be calculated across all dimensions.

        Optional Extra parameters: "temperature", "mask".

        If no temperature values are provided, temperature will be calculated using a gradient specified by the input_params dictionary and the depth variable.

        If mask values are provided, then any non-nan values in the mask variable will be set to zero across all variables in the returned output_data object.

    :param utc_properties:

    :param rng_seed:

    :param p_values:

        A list of p_values for the doublet calculation to perform over; if no p_values are provided then the default value of P50 is used.

    :return output_data:

        A xr.Dataset (input_data) with the resulting variables from the doublet calculation:

        "power", "heat_pump_power", "capex", "opex", "utc", "npv", "hprod", "cop", "cophp", "pres", "flow_rate", "welld"



        The output will have the same dimensions as the input_data class, with the additional p_value dimension

    Perform a ThermoGIS Doublet performance simulation.



    This function computes doublet performance metrics across all dimensions of the input dataset.

    The input can be scalar, 1D, or 2D gridded data. If no temperature values are provided, they

    are estimated from a gradient defined in `input_params`. If a mask is provided, any non-NaN

    values will result in zeroing the output values at those locations.



    Parameters

    ----------

    input_data : xr.Dataset

        An xarray Dataset containing the required input variables:

        - "thickness_mean"

        - "thickness_sd"

        - "porosity"

        - "ntg"

        - "depth"

        - "ln_permeability_mean"

        - "ln_permeability_sd"



        Optional variables:

        - "temperature" : If not provided, temperature is estimated using the depth and a temperature gradient from `input_params`.

        - "mask" : If provided, all non-NaN values will result in setting corresponding output values to zero.



    utc_properties : dict

        A dictionary of UTC (Underground Thermal Capacity) properties used for simulation.



    rng_seed : int

        Random seed used for stochastic components of the simulation.



    p_values : list of float, optional

        List of probability values (e.g., [0.1, 0.5, 0.9]) for the performance evaluation.

        If not provided, the default value of P50 (0.5) is used.



    Returns

    -------

    output_data : xr.Dataset

        An xarray Dataset with the same spatial dimensions as `input_data`, plus an added `p_value` dimension.

        Contains the following output variables:

        - "power"

        - "heat_pump_power"

        - "capex"

        - "opex"

        - "utc"

        - "npv"

        - "hprod"

        - "cop"

        - "cophp"

        - "pres"

        - "flow_rate"

        - "welld"

    """



    # Check that all essential variables are provided
@@ -112,9 +146,16 @@ def calculate_doublet_performance(input_data: xr.Dataset, 


def validate_input_data(input_data: xr.Dataset):

    """

    Ensure that the input_data Dataset contains the minimum required variables

    :param input_data:

    :return:

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



    Parameters

    ----------

    input_data : xr.Dataset

        Input dataset that must contain the required ThermoGIS variables.



    Returns

    -------

    None

    """

    missing_variables = []

    for variable in ["thickness_mean", "thickness_sd", "porosity", "ntg", "depth", "ln_permeability_mean", "ln_permeability_sd"]:
@@ -135,29 +176,45 @@ def validate_input_data(input_data: xr.Dataset): 
def calculate_performance_of_single_location(mask: float, depth: float, thickness: float, porosity: float, ntg: float, temperature: float, transmissivity: float, transmissivity_with_ntg: float, doublet: JClass = None ,

                                             utc_properties: JClass = None) -> float:

    """

        Calculate the performance of a doublet at a single location

    :param utc_properties:

    :param mask:

        mask value, if not nan set all output to 0.0

    :param depth:

        depth of the aquifer

    :param thickness:

        thickness of the aquifer

    :param porosity:

        porosity of the aquifer

    :param ntg:

        net-to-gross of the aquifer

    :param temperature:

        temperature of the aquifer

    :param transmissivity:

        transmissivity of the aquifer

    :param transmissivity_with_ntg:

        transmissivity * ntg of the aquifer

    :param doublet:

        a ThermoGIS doublet class

    :return:

        the values "power", "heat_pump_power", "capex", "opex", "utc", "npv", "hprod", "cop", "cophp", "pres", "flow_rate", "welld" from the doublet calculation



    Calculate the performance of a doublet at a single location.



    Parameters

    ----------

    utc_properties : dict

        Dictionary containing UTC (Underground Thermal Capacity) properties.

    mask : float

        Mask value; if not NaN, all output values will be set to 0.0.

    depth : float

        Depth of the aquifer in meters.

    thickness : float

        Thickness of the aquifer in meters.

    porosity : float

        Porosity of the aquifer (fraction).

    ntg : float

        Net-to-gross ratio of the aquifer (fraction).

    temperature : float

        Temperature of the aquifer in degrees Celsius.

    transmissivity : float

        Transmissivity of the aquifer.

    transmissivity_with_ntg : float

        Product of transmissivity and net-to-gross.

    doublet : object

        An instance of the ThermoGIS doublet class.



    Returns

    -------

    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):
@@ -206,12 +263,21 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 


def instantiate_thermogis_doublet(utc_properties, rng_seed: int = None) -> JClass:

    """

    Instantiate a ThermoGIS Doublet class, with a set random seed if provided

    :param utc_properties: A UTCProperties instance, with different input parameters

    :param rng_seed: An integer to set the random number generator

    :return:

        doublet, a JPype instantiated object of the ThermoGISDoublet class

    Instantiate a ThermoGIS Doublet class, optionally using a specified random seed.



    Parameters

    ----------

    utc_properties : UTCProperties

        Instance containing the input parameters required for the doublet simulation.

    rng_seed : int

        Optional random seed for reproducibility.



    Returns

    -------

    doublet : object

        A JPype-instantiated object of the ThermoGISDoublet class.

    """



    # Instantiate doublet class

    Logger = JClass("logging.Logger")

    Mockito = JClass("org.mockito.Mockito")
@@ -235,16 +301,30 @@ def instantiate_thermogis_doublet(utc_properties, rng_seed: int = None) -> JClas 


def set_doublet_parameters(doublet, transmissivity_with_ntg: float, depth: float, porosity: float, ntg: float, temperature: float, utc_properties: JClass):

    """

    For a single location sets the necessary data on the doublet class, to then run a doublet simulation

    :param utc_properties:

    :param doublet:

    :param transmissivity_with_ntg:

    :param depth:

    :param porosity:

    :param ntg:

    :param temperature:

    :return:

    Set necessary data on the doublet class for a single location prior to simulation.



    Parameters

    ----------

    utc_properties : dict

        Dictionary containing UTC properties.

    doublet : object

        An instance of the ThermoGIS doublet class.

    transmissivity_with_ntg : float

        Product of transmissivity and net-to-gross.

    depth : float

        Depth of the aquifer in meters.

    porosity : float

        Porosity of the aquifer (fraction).

    ntg : float

        Net-to-gross ratio of the aquifer (fraction).

    temperature : float

        Temperature of the aquifer in degrees Celsius.



    Returns

    -------

    None

    """



    if not utc_properties.useStimulation() or transmissivity_with_ntg > utc_properties.stimKhMax():

        doublet.setNoStimulation()
@@ -260,12 +340,3 @@ def set_doublet_parameters(doublet, transmissivity_with_ntg: float, depth: float 
            doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - utc_properties.maxCoolingTempRange(), utc_properties.hpMinimumInjectionTemperature()]))

    else:

        doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - utc_properties.maxCoolingTempRange(), utc_properties.dhReturnTemp()]))
 
 No newline at end of file
 


    # # old:

    # if not utc_properties.useHeatPump():

    #     doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - utc_properties.maxCoolingTempRange(), utc_properties.dhReturnTemp()]))

    # elif utc_properties.useHeatPump() and utc_properties.calculateCop() and not utc_properties.hpApplicationMode():

    #     doublet.doubletCalc1DData.setInjectionTemp(doublet.calculateInjectionTempWithHeatPump(temperature, utc_properties.hpDirectHeatInputTemp()))

    # else:

    #     doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - utc_properties.maxCoolingTempRange(), utc_properties.hpMinimumInjectionTemperature()]))