TNO Intern

Commit 2791e941 authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Merge branch '55-fix-example-6' into 'main' 

Resolve "fix example 6"

Closes #55

See merge request AGS/pythermogis!70
parents 0b3d3bf5 ea9238a7

docs/images/npv.png

−800 B (213 KiB)
Loading image diff...

docs/index.md

+2 −0
Original line number Diff line number Diff line
@@ -4,9 +4,11 @@ 


<div style="display: flex; gap: 10px;">

  <img src="../images/pyThermoGIS_transparent.png" style="width: 15%; max-width: 400px; height: auto;">

  <img src="../images/GeologischeDienstlogo.png" alt="Image 1" style="width: 50%; max-width: 400px; height: auto;">

  <img src="../images/TNOLogo.jpg" alt="Image 1" style="width: 25%; max-width: 400px; height: auto;">

</div>



The intended usage of this pacakge is to import into your own projects so you can run your own doublet simulations using your own input parameters or statistical frameworks.



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

docs/usage/contents.md

+1 −1
Original line number Diff line number Diff line 
#  Pythermogis Usage - Examples



Through a number examples the core functionality of **pyThermoGIS** is explained.  Many aspetcs 

Through a number examples the core functionality of **pyThermoGIS** is explained.  Many aspects 

of the these examples can also be found in the  tests.



The underlying theory of the geothermal doublet simulation is explained in the [theory section](../theory/introduction.md).

docs/usage/portfoliorun_analysis.md

+71 −63
Original line number Diff line number Diff line
@@ -34,17 +34,18 @@ This example corresponds to test case `test_example6` in `test_doc_examples.py` 
directory of the repository.



```python

from pythermogis import calculate_doublet_performance, calculate_pos_pvalues_singlelocation

from pythermogis import calculate_doublet_performance, calculate_pos

from pygridsio import read_grid

from matplotlib import pyplot as plt

from numpy import np

import xarray as xr

from pathlib import Path

from os import path



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)



# read in reservoir property maps

reservoir_properties = read_grid(input_data_path / "ROSL_ROSLU__thick.zmap").to_dataset(name="thickness_mean")

reservoir_properties["thickness_sd"] = read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap")

reservoir_properties["ntg"] = read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap")
@@ -53,70 +54,77 @@ reservoir_properties["depth"] = read_grid(input_data_path / "ROSL_ROSLU__top.zma 
reservoir_properties["ln_permeability_mean"] = np.log( read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"))

reservoir_properties["ln_permeability_sd"] =  read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap")



portfolioloc = [(125e3, 525e3), (100e3, 525e3), (110e3, 525e3), (125e3, 515e3),

              (125e3, 520e3)]  # define location

# Select locations of interest from the reservoir properties and simulate doublet performance

portfolio_coords = [(125e3, 525e3), (100e3, 525e3), (110e3, 525e3), (125e3, 515e3), (125e3, 520e3)]



# Collect selected locations, create a new Dataset with the coordinate 'location'

locations = []

for i, (x, y) in enumerate(portfolio_coords):

  point = reservoir_properties.sel(x=x, y=y, method="nearest")

  point = point.expand_dims(location=[i])  # Add new dim for stacking

  locations.append(point)

portfolio_reservoir_properties = xr.concat(locations, dim="location")



fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))

colors = plt.cm.tab10.colors

ic = 0

for x, y in portfolioloc:

   reservoir_properties_loc = reservoir_properties.sel(x=x, y=y, method="nearest")

results_portfolio = calculate_doublet_performance(portfolio_reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90])



   results_loc = calculate_doublet_performance(reservoir_properties_loc,

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

# post process: clip the npv and calculate the probability of success

AEC = -3.5

   results_loc['npv'] = results_loc.npv.clip(min=AEC)

results_portfolio['npv'] = results_portfolio.npv.clip(min=AEC)

results_portfolio["pos"] = calculate_pos(results_portfolio)



# plot results

fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))

colors = plt.cm.tab10.colors

for i_loc, loc in enumerate(portfolio_coords):



   # probability of success is defined as the p-value where npv = 0.0, use interpolation to find pos:

   pos = calculate_pos_pvalues_singlelocation(results_loc.npv,

                                                 results_loc.p_value)  # The order of the npv data has to be reversed as np.interp requires values to increase to function properly

  # select single portfolio location to plot

  results_loc = results_portfolio.isel(location=i_loc)

  pos = results_loc.pos.data

  x, y = loc



  # plot npv versus p-value and a map showing the location of interest

   axe = axes[0, 0]

   results_loc.npv.plot(y="p_value", ax=axe, color=colors[ic])

   axe.set_title(f"Aquifer: ROSL_ROSLU\n  POS ")

   axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

   axe.axvline(0.0, ls="--", c=colors[ic])

   axe.legend()

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

   axe.set_ylabel("p-value [%]")

   

   axe = axes[0, 1]

   kh = results_loc.transmissivity_with_ntg * 1.0

   kh.plot(y="p_value", ax=axe, color=colors[ic])

   axe.set_title(f"Aquifer: ROSL_ROSLU\n kH")

   temp = reservoir_properties_loc.temperature

   inj_temp = results_loc.inj_temp.sel(p_value=50, method="nearest")

   prd_temp = results_loc.prd_temp.sel(p_value=50, method="nearest")

   axe.axhline(50.0, label=f"POS: {pos:.1f}%  TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ",

              ls="--", c=colors[ic])

   kh50 = kh.sel(p_value=50, method="nearest")

   axe.axvline(kh50, ls="--", c=colors[ic])

   axe.legend()

   axe.set_xlabel("kH [Dm]")

   axe.set_ylabel("p-value [%]")

   

   axe = axes[1, 0]

   results_loc.power.plot(y="p_value", ax=axe, color=colors[ic])

   axe.set_title(f"Aquifer: ROSL_ROSLU\n Power")

   axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

   valp50 = results_loc.power.sel(p_value=50, method="nearest")

   axe.axvline(valp50, ls="--", c=colors[ic])

   axe.legend()

   axe.set_xlabel("power [MW]")

   axe.set_ylabel("p-value [%]")

   

   axe = axes[1, 1]

   results_loc.cop.plot(y="p_value", ax=axe, color=colors[ic])

   axe.set_title(f"Aquifer: ROSL_ROSLU\n COP")

   axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

   valp50 = results_loc.cop.sel(p_value=50, method="nearest")

   axe.axvline(valp50, ls="--", c=colors[ic])

   axe.legend()

   axe.set_xlabel("COP [-]")

   axe.set_ylabel("p-value [%]")

   ic += 1

  ax = axes[0, 0]

  results_loc.npv.plot(y="p_value", ax=ax, color=colors[i_loc])

  ax.set_title(f"Aquifer: ROSL_ROSLU\n  POS ")

  ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

  ax.axvline(0.0, ls="--", c=colors[i_loc])

  ax.legend()

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

  ax.set_ylabel("p-value [%]")



  ax = axes[0, 1]

  kh = results_loc.transmissivity_with_ntg

  kh.plot(y="p_value", ax=ax, color=colors[i_loc])

  ax.set_title(f"Aquifer: ROSL_ROSLU\n kH")

  temp = results_loc.temperature.sel(p_value=50)

  inj_temp = results_loc.inj_temp.sel(p_value=50)

  prd_temp = results_loc.prd_temp.sel(p_value=50)

  ax.axhline(50.0, label=f"POS: {pos:.1f}%  TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ",

              ls="--", c=colors[i_loc])

  ax.axvline(kh.sel(p_value=50), ls="--", c=colors[i_loc])

  ax.legend()

  ax.set_xlabel("kH [Dm]")

  ax.set_ylabel("p-value [%]")



  ax = axes[1, 0]

  results_loc.power.plot(y="p_value", ax=ax, color=colors[i_loc])

  ax.set_title(f"Aquifer: ROSL_ROSLU\n Power")

  ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

  valp50 = results_loc.power.sel(p_value=50)

  ax.axvline(valp50, ls="--", c=colors[i_loc])

  ax.legend()

  ax.set_xlabel("power [MW]")

  ax.set_ylabel("p-value [%]")



  ax = axes[1, 1]

  results_loc.cop.plot(y="p_value", ax=ax, color=colors[i_loc])

  ax.set_title(f"Aquifer: ROSL_ROSLU\n COP")

  ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

  valp50 = results_loc.cop.sel(p_value=50)

  ax.axvline(valp50, ls="--", c=colors[i_loc])

  ax.legend()

  ax.set_xlabel("COP [-]")

  ax.set_ylabel("p-value [%]")



plt.tight_layout()  # ensure there is enough spacing

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

tests/test_doc_examples.py

+67 −60
Original line number Diff line number Diff line
@@ -161,18 +161,18 @@ def test_example_5(): 




def test_example_6():

    from pythermogis import calculate_doublet_performance, calculate_pos_pvalues_singlelocation

    from pythermogis import calculate_doublet_performance, calculate_pos

    from pygridsio import read_grid

    from matplotlib import pyplot as plt

    import xarray as xr

    from pathlib import Path

    from os import path

    from numpy import np



    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)



    # read in reservoir property maps

    reservoir_properties = read_grid(input_data_path / "ROSL_ROSLU__thick.zmap").to_dataset(name="thickness_mean")

    reservoir_properties["thickness_sd"] = read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap")

    reservoir_properties["ntg"] = read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap")
@@ -181,70 +181,77 @@ def test_example_6(): 
    reservoir_properties["ln_permeability_mean"] = np.log( read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"))

    reservoir_properties["ln_permeability_sd"] =  read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap")



    portfolioloc = [(125e3, 525e3), (100e3, 525e3), (110e3, 525e3), (125e3, 515e3),

                    (125e3, 520e3)]  # define location

    # Select locations of interest from the reservoir properties and simulate doublet performance

    portfolio_coords = [(125e3, 525e3), (100e3, 525e3), (110e3, 525e3), (125e3, 515e3), (125e3, 520e3)]



    # Collect selected locations, create a new Dataset with the coordinate 'location'

    locations = []

    for i, (x, y) in enumerate(portfolio_coords):

        point = reservoir_properties.sel(x=x, y=y, method="nearest")

        point = point.expand_dims(location=[i])  # Add new dim for stacking

        locations.append(point)

    portfolio_reservoir_properties = xr.concat(locations, dim="location")



    fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))

    colors = plt.cm.tab10.colors

    ic = 0

    for x, y in portfolioloc:

        reservoir_properties_loc = reservoir_properties.sel(x=x, y=y, method="nearest")

    results_portfolio = calculate_doublet_performance(portfolio_reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90])



        results_loc = calculate_doublet_performance(reservoir_properties_loc,

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

    # post process: clip the npv and calculate the probability of success

    AEC = -3.5

        results_loc['npv'] = results_loc.npv.clip(min=AEC)



        # probability of success is defined as the p-value where npv = 0.0, use interpolation to find pos:

        pos = calculate_pos_pvalues_singlelocation(results_loc.npv,

                                                       results_loc.p_value)  # The order of the npv data has to be reversed as np.interp requires values to increase to function properly

    results_portfolio['npv'] = results_portfolio.npv.clip(min=AEC)

    results_portfolio["pos"] = calculate_pos(results_portfolio)



        # plot npv versus p-value and a map showing the location of interest

        axe = axes[0, 0]

        results_loc.npv.plot(y="p_value", ax=axe, color=colors[ic])

        axe.set_title(f"Aquifer: ROSL_ROSLU\n  POS ")

        axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

        axe.axvline(0.0, ls="--", c=colors[ic])

        axe.legend()

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

        axe.set_ylabel("p-value [%]")



        axe = axes[0, 1]

        kh = results_loc.transmissivity_with_ntg * 1.0

        kh.plot(y="p_value", ax=axe, color=colors[ic])

        axe.set_title(f"Aquifer: ROSL_ROSLU\n kH")

        temp = reservoir_properties_loc.temperature

        inj_temp = results_loc.inj_temp.sel(p_value=50, method="nearest")

        prd_temp = results_loc.prd_temp.sel(p_value=50, method="nearest")

        axe.axhline(50.0, label=f"POS: {pos:.1f}%  TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ",

                    ls="--", c=colors[ic])

        kh50 = kh.sel(p_value=50, method="nearest")

        axe.axvline(kh50, ls="--", c=colors[ic])

        axe.legend()

        axe.set_xlabel("kH [Dm]")

        axe.set_ylabel("p-value [%]")

    # plot results

    fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))

    colors = plt.cm.tab10.colors

    for i_loc, loc in enumerate(portfolio_coords):



        axe = axes[1, 0]

        results_loc.power.plot(y="p_value", ax=axe, color=colors[ic])

        axe.set_title(f"Aquifer: ROSL_ROSLU\n Power")

        axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

        valp50 = results_loc.power.sel(p_value=50, method="nearest")

        axe.axvline(valp50, ls="--", c=colors[ic])

        axe.legend()

        axe.set_xlabel("power [MW]")

        axe.set_ylabel("p-value [%]")

        # select single portfolio location to plot

        results_loc = results_portfolio.isel(location=i_loc)

        pos = results_loc.pos.data

        x, y = loc



        axe = axes[1, 1]

        results_loc.cop.plot(y="p_value", ax=axe, color=colors[ic])

        axe.set_title(f"Aquifer: ROSL_ROSLU\n COP")

        axe.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[ic])

        valp50 = results_loc.cop.sel(p_value=50, method="nearest")

        axe.axvline(valp50, ls="--", c=colors[ic])

        axe.legend()

        axe.set_xlabel("COP [-]")

        axe.set_ylabel("p-value [%]")

        ic += 1

        # plot npv versus p-value and a map showing the location of interest

        ax = axes[0, 0]

        results_loc.npv.plot(y="p_value", ax=ax, color=colors[i_loc])

        ax.set_title(f"Aquifer: ROSL_ROSLU\n  POS ")

        ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

        ax.axvline(0.0, ls="--", c=colors[i_loc])

        ax.legend()

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

        ax.set_ylabel("p-value [%]")



        ax = axes[0, 1]

        kh = results_loc.transmissivity_with_ntg

        kh.plot(y="p_value", ax=ax, color=colors[i_loc])

        ax.set_title(f"Aquifer: ROSL_ROSLU\n kH")

        temp = results_loc.temperature.sel(p_value=50)

        inj_temp = results_loc.inj_temp.sel(p_value=50)

        prd_temp = results_loc.prd_temp.sel(p_value=50)

        ax.axhline(50.0, label=f"POS: {pos:.1f}%  TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ",

                    ls="--", c=colors[i_loc])

        ax.axvline(kh.sel(p_value=50), ls="--", c=colors[i_loc])

        ax.legend()

        ax.set_xlabel("kH [Dm]")

        ax.set_ylabel("p-value [%]")



        ax = axes[1, 0]

        results_loc.power.plot(y="p_value", ax=ax, color=colors[i_loc])

        ax.set_title(f"Aquifer: ROSL_ROSLU\n Power")

        ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

        valp50 = results_loc.power.sel(p_value=50)

        ax.axvline(valp50, ls="--", c=colors[i_loc])

        ax.legend()

        ax.set_xlabel("power [MW]")

        ax.set_ylabel("p-value [%]")



        ax = axes[1, 1]

        results_loc.cop.plot(y="p_value", ax=ax, color=colors[i_loc])

        ax.set_title(f"Aquifer: ROSL_ROSLU\n COP")

        ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc])

        valp50 = results_loc.cop.sel(p_value=50)

        ax.axvline(valp50, ls="--", c=colors[i_loc])

        ax.legend()

        ax.set_xlabel("COP [-]")

        ax.set_ylabel("p-value [%]")



    plt.tight_layout()  # ensure there is enough spacing

    plt.savefig(output_data_path / "npv.png")
 
 No newline at end of file