Updating documentation and example 6 (63fb22cc) · Commits · AGS_public / pythermogis

docs/images/npv.png

−800 B (213 KiB)

Loading image diff...

docs/usage/portfoliorun_analysis.md

+94 −86

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