docs examples to new usage (23e40f76) · Commits · AGS_public / pythermogis

pythermogis/cli.py

+2 −2

Original line number	Original line	Diff line number	Diff line
	@@ -3,8 +3,8 @@ from pathlib import Path

	import typer		import typer

	from pytg3.aquifer import Aquifer		from pythermogis.pytg3.aquifer import Aquifer
	from pytg3.doublet import ThermoGISDoublet		from pythermogis.pytg3.doublet import ThermoGISDoublet
	from pythermogis.potential import calculate_potential		from pythermogis.potential import calculate_potential
	from pythermogis.properties import instantiate_thermogis_parameters		from pythermogis.properties import instantiate_thermogis_parameters

pythermogis/pytg3/aquifer.py

+11 −0

Original line number	Original line	Diff line number	Diff line
	@@ -48,12 +48,23 @@ class BaseAquifer(BaseModel):
	) -> None:		) -> None:
	raise NotImplementedError		raise NotImplementedError

			def to_dataset(self) -> xr.Dataset:
			return xr.Dataset(
			{
			field_name: value
			for field_name, value in self
			if isinstance(value, xr.DataArray)
			}
			)



	class Aquifer(BaseAquifer):		class Aquifer(BaseAquifer):
	thickness: FloatOrArray		thickness: FloatOrArray
	permeability: FloatOrArray \| None = None		permeability: FloatOrArray \| None = None
	transmissivity: FloatOrArray \| None = None		transmissivity: FloatOrArray \| None = None


	@model_validator(mode="after")		@model_validator(mode="after")
	def check_permeability_or_transmissivity(self) -> Aquifer:		def check_permeability_or_transmissivity(self) -> Aquifer:
	has_perm = self.permeability is not None		has_perm = self.permeability is not None

tests/test_doc_examples.py

+120 −107

Original line number	Original line	Diff line number	Diff line
	@@ -6,29 +6,23 @@ import xarray as xr
	from matplotlib import pyplot as plt		from matplotlib import pyplot as plt
	from pygridsio import plot_grid, read_grid, resample_grid		from pygridsio import plot_grid, read_grid, resample_grid

	from pythermogis.doublet import (		from pythermogis.pytg3.aquifer import Aquifer, StochasticAquifer
	calculate_doublet_performance,		from pythermogis.pytg3.settings import Settings
	calculate_doublet_performance_stochastic,		from pythermogis.pytg3.doublet import ThermoGISDoublet
	)
	from pythermogis.plot import plot_exceedance		from pythermogis.plot import plot_exceedance
	from pythermogis.postprocess import calculate_pos		from pythermogis.postprocess import calculate_pos
	from pythermogis.properties import instantiate_thermogis_parameters


	def test_example_1():		def test_example_1():
	import xarray as xr		aquifer = Aquifer(
			thickness=300,
	input_data = xr.Dataset(		ntg=0.5,
	{		porosity=0.2,
	"thickness": 300,		depth=2000,
	"ntg": 0.5,		permeability=300
	"porosity": 0.2,
	"depth": 2000,
	"permeability": 300,
	}
	)		)

	results = calculate_doublet_performance(input_data)		results = ThermoGISDoublet(aquifer=aquifer).simulate().to_dataset()
	print(results)		print(results)


	@@ -44,42 +38,38 @@ def test_example_2():
	coords={"x": [0, 1], "y": [10, 20]},		coords={"x": [0, 1], "y": [10, 20]},
	)		)

	results = calculate_doublet_performance(input_data)		aquifer = Aquifer(
			thickness=input_data["thickness"],
			ntg=input_data["ntg"],
			porosity=input_data["porosity"],
			depth=input_data["depth"],
			permeability=input_data["permeability"],
			)


			results = ThermoGISDoublet(aquifer=aquifer).simulate().to_dataset()
	print(results)		print(results)


	def test_example_3():		def test_example_3():
	input_data = xr.Dataset(		aquifer = Aquifer(
	{		thickness=300, ntg=0.5, porosity=0.2, depth=2000, permeability=300
	"thickness": 300,		)
	"ntg": 0.5,
	"porosity": 0.2,		settings = Settings()
	"depth": 2000,		settings.set_use_heatpump(True)
	"permeability": 300,
	}		results = ThermoGISDoublet(aquifer=aquifer, settings=settings).simulate().to_dataset()
	)

	tg_properties = instantiate_thermogis_parameters()
	tg_properties.setUseHeatpump(True)
	utc_properties = tg_properties.setupUTCParameters()
	results = calculate_doublet_performance(input_data, utc_properties=utc_properties)
	print(results)		print(results)


	def test_example_4():		def test_example_4():
	input_data = xr.Dataset(		aquifer = StochasticAquifer(
	{		thickness_mean=300, thickness_sd=50, ntg=0.5, porosity=0.2, depth=2000, ln_permeability_mean=4, ln_permeability_sd=0.5
	"thickness_mean": 300,
	"thickness_sd": 50,
	"ntg": 0.5,
	"porosity": 0.2,
	"depth": 2000,
	"ln_permeability_mean": 4,
	"ln_permeability_sd": 0.5,
	}
	)		)
			results = ThermoGISDoublet(aquifer=aquifer).simulate(p_values=[50]).to_dataset()

	results = calculate_doublet_performance_stochastic(input_data)
	print(results)		print(results)


	@@ -89,25 +79,24 @@ def test_example_5():
	)		)
	output_data_path.mkdir(parents=True, exist_ok=True)		output_data_path.mkdir(parents=True, exist_ok=True)

	input_data = xr.Dataset(		aquifer = StochasticAquifer(
	{		thickness_mean=300,
	"thickness_mean": 300,		thickness_sd=50,
	"thickness_sd": 50,		ntg=0.5,
	"ntg": 0.5,		porosity=0.2,
	"porosity": 0.2,		depth=2000,
	"depth": 2000,		ln_permeability_mean=5,
	"ln_permeability_mean": 5,		ln_permeability_sd=0.5,
	"ln_permeability_sd": 0.5,
	}
	)
	results = calculate_doublet_performance_stochastic(
	input_data, p_values=[1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99]
	)		)

			results = ThermoGISDoublet(aquifer=aquifer).simulate(p_values=[1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99]).to_dataset()


	fig, axe = plt.subplots(figsize=(10, 5))		fig, axe = plt.subplots(figsize=(10, 5))
	kh = results.transmissivity_with_ntg * 1.0		kh = results.transmissivity_with_ntg * 1.0
	kh.plot(y="p_value", ax=axe)		kh.plot(y="p_value", ax=axe)
	axe.set_title("Aquifer Net transmissivity\n kH")		axe.set_title("Aquifer Net transmissivity\n kH")
	temp = input_data.temperature		temp = aquifer.temperature
	inj_temp = results.inj_temp.sel(p_value=50, method="nearest")		inj_temp = results.inj_temp.sel(p_value=50, method="nearest")
	prd_temp = results.prd_temp.sel(p_value=50, method="nearest")		prd_temp = results.prd_temp.sel(p_value=50, method="nearest")
	axe.axhline(		axe.axhline(
	@@ -129,12 +118,12 @@ def test_example_6():
	# the location of the input data:		# the location of the input data:
	# the data can be found in the resources/example_data directory of the repo		# the data can be found in the resources/example_data directory of the repo
	input_data_path = (		input_data_path = (
	Path(__file__).parent.parent / "resources" / "test_input" / "example_data"		Path(__file__).parent / "resources" / "test_input" / "example_data"
	)		)

	# create a directory to write the output files to		# create a directory to write the output files to
	output_data_path = (		output_data_path = (
	Path(__file__).parent.parent / "resources" / "test_output" / "example6"		Path(__file__).parent / "resources" / "test_output" / "example6"
	)		)
	output_data_path.mkdir(parents=True, exist_ok=True)		output_data_path.mkdir(parents=True, exist_ok=True)

	@@ -145,37 +134,47 @@ def test_example_6():
	# grids can be in .nc, .asc, .zmap or .tif format: see pygridsio package		# 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;		new_cellsize = 5000 # in m, this sets the resolution of the model;
	# to speed up calcualtions you can increase the cellsize		# to speed up calcualtions you can increase the cellsize
	reservoir_properties = resample_grid(		thickness_mean = resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize		read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize
	).to_dataset(name="thickness_mean")		)
	reservoir_properties["thickness_sd"] = resample_grid(		thickness_sd = resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"),		read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"),
	new_cellsize=new_cellsize,		new_cellsize=new_cellsize,
	)		)
	reservoir_properties["ntg"] = resample_grid(		ntg = resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize		read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize
	)		)
	reservoir_properties["porosity"] = (		porosity = (
	resample_grid(		resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"),		read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"),
	new_cellsize=new_cellsize,		new_cellsize=new_cellsize,
	)		)
	/ 100		/ 100
	)		)
	reservoir_properties["depth"] = resample_grid(		depth = resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize		read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize
	)		)
	reservoir_properties["ln_permeability_mean"] = np.log(		ln_permeability_mean = np.log(
	resample_grid(		resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"),		read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"),
	new_cellsize=new_cellsize,		new_cellsize=new_cellsize,
	)		)
	)		)
	reservoir_properties["ln_permeability_sd"] = resample_grid(		ln_permeability_sd = resample_grid(
	read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"),		read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"),
	new_cellsize=new_cellsize,		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,
			)

	# output inputs		# output inputs
	variables_to_plot = ["depth", "thickness_mean", "thickness_sd"]		variables_to_plot = ["depth", "thickness_mean", "thickness_sd"]
	fig, axes = plt.subplots(		fig, axes = plt.subplots(
	@@ -187,7 +186,7 @@ def test_example_6():
	)		)
	for i, variable in enumerate(variables_to_plot):		for i, variable in enumerate(variables_to_plot):
	plot_grid(		plot_grid(
	reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True		getattr(aquifer, variable), axes=axes[i], add_netherlands_shapefile=True
	) # See documentation on plot_grid in pygridsio,		) # See documentation on plot_grid in pygridsio,
	# you can also provide your own shapefile		# you can also provide your own shapefile
	plt.tight_layout() # ensure there is enough spacing		plt.tight_layout() # ensure there is enough spacing
	@@ -203,7 +202,7 @@ def test_example_6():
	)		)
	for i, variable in enumerate(variables_to_plot):		for i, variable in enumerate(variables_to_plot):
	plot_grid(		plot_grid(
	reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True		getattr(aquifer, variable), axes=axes[i], add_netherlands_shapefile=True
	) # See documentation on plot_grid in pygridsio,		) # See documentation on plot_grid in pygridsio,
	# you can also provide your own shapefile		# you can also provide your own shapefile
	plt.tight_layout() # ensure there is enough spacing		plt.tight_layout() # ensure there is enough spacing
	@@ -215,9 +214,8 @@ def test_example_6():
	if (output_data_path / "output_results.nc").exists and not run_simulation:		if (output_data_path / "output_results.nc").exists and not run_simulation:
	results = xr.load_dataset(output_data_path / "output_results.nc")		results = xr.load_dataset(output_data_path / "output_results.nc")
	else:		else:
	results = calculate_doublet_performance_stochastic(
	reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]		results = ThermoGISDoublet(aquifer=aquifer).simulate(p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]).to_dataset()
	)
	results.to_netcdf(		results.to_netcdf(
	output_data_path / "output_results.nc"		output_data_path / "output_results.nc"
	) # write doublet simulation results to file		) # write doublet simulation results to file
	@@ -284,22 +282,28 @@ def test_example_7():
	output_data_path.mkdir(parents=True, exist_ok=True)		output_data_path.mkdir(parents=True, exist_ok=True)

	# read in reservoir property maps		# read in reservoir property maps
	reservoir_properties = read_grid(		thickness_mean = read_grid(
	input_data_path / "ROSL_ROSLU__thick.zmap"		input_data_path / "ROSL_ROSLU__thick.zmap"
	).to_dataset(name="thickness_mean")		)
	reservoir_properties["thickness_sd"] = read_grid(		thickness_sd = read_grid(
	input_data_path / "ROSL_ROSLU__thick_sd.zmap"		input_data_path / "ROSL_ROSLU__thick_sd.zmap"
	)		)
	reservoir_properties["ntg"] = read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap")		ntg = read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap")
	reservoir_properties["porosity"] = 0.01 * read_grid(		porosity = 0.01 * read_grid(
	input_data_path / "ROSL_ROSLU__poro.zmap"		input_data_path / "ROSL_ROSLU__poro.zmap"
	)		)
	reservoir_properties["depth"] = read_grid(input_data_path / "ROSL_ROSLU__top.zmap")		depth = read_grid(input_data_path / "ROSL_ROSLU__top.zmap")
	reservoir_properties["ln_permeability_mean"] = np.log(		ln_permeability_mean = np.log(read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"))
	read_grid(input_data_path / "ROSL_ROSLU__perm.zmap")		ln_permeability_sd = read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap")
	)
	reservoir_properties["ln_permeability_sd"] = read_grid(		aquifer = StochasticAquifer(
	input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"		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,
	)		)

	# Select locations of interest from the reservoir		# Select locations of interest from the reservoir
	@@ -313,16 +317,18 @@ def test_example_7():
	]		]

	# Collect selected locations, create a new Dataset with the coordinate 'location'		# Collect selected locations, create a new Dataset with the coordinate 'location'
			aquifer_ds = aquifer.to_dataset()

	locations = []		locations = []
	for i, (x, y) in enumerate(portfolio_coords):		for i, (x, y) in enumerate(portfolio_coords):
	point = reservoir_properties.sel(x=x, y=y, method="nearest")		point = aquifer_ds.sel(x=x, y=y, method="nearest")
	point = point.expand_dims(location=[i]) # Add new dim for stacking		point = point.expand_dims(location=[i])
	locations.append(point)		locations.append(point)
	portfolio_reservoir_properties = xr.concat(locations, dim="location")		portfolio_reservoir_properties = xr.concat(locations, dim="location")

	results_portfolio = calculate_doublet_performance_stochastic(		portfolio_aquifer = StochasticAquifer(**portfolio_reservoir_properties.data_vars)
	portfolio_reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]
	)		results_portfolio = ThermoGISDoublet(aquifer=portfolio_aquifer).simulate(p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]).to_dataset()

	# post process: clip the npv and calculate the probability of success		# post process: clip the npv and calculate the probability of success
	AEC = -3.5		AEC = -3.5
	@@ -357,7 +363,7 @@ def test_example_7():
	kh = results_loc.transmissivity_with_ntg		kh = results_loc.transmissivity_with_ntg
	kh.plot(y="p_value", ax=ax, color=colors[i_loc])		kh.plot(y="p_value", ax=ax, color=colors[i_loc])
	ax.set_title("Aquifer: ROSL_ROSLU\n kH")		ax.set_title("Aquifer: ROSL_ROSLU\n kH")
	temp = results_loc.temperature.sel(p_value=50)		temp = results_loc.temperature
	inj_temp = results_loc.inj_temp.sel(p_value=50)		inj_temp = results_loc.inj_temp.sel(p_value=50)
	prd_temp = results_loc.prd_temp.sel(p_value=50)		prd_temp = results_loc.prd_temp.sel(p_value=50)
	ax.axhline(		ax.axhline(
	@@ -409,7 +415,7 @@ def test_example_7():
	def test_example_8():		def test_example_8():
	# output directory to write output to		# output directory to write output to
	output_data_path = (		output_data_path = (
	Path(__file__).parent.parent / "resources" / "test_output" / "example8"		Path(__file__).parent / "resources" / "test_output" / "example8"
	)		)
	output_data_path.mkdir(parents=True, exist_ok=True)		output_data_path.mkdir(parents=True, exist_ok=True)

	@@ -432,9 +438,14 @@ def test_example_8():
	)		)

	# For every sample, run a doublet simulation store the output values		# For every sample, run a doublet simulation store the output values
	simulations = calculate_doublet_performance(		aquifer = Aquifer(
	reservoir_properties, print_execution_duration=True		thickness=reservoir_properties["thickness"],
			porosity=reservoir_properties["porosity"],
			ntg=reservoir_properties["ntg"],
			depth=reservoir_properties["depth"],
			permeability=reservoir_properties["permeability"],
	)		)
			simulations = ThermoGISDoublet(aquifer=aquifer).simulate().to_dataset()

	# post process the samples to calculate the percentiles of each variable		# post process the samples to calculate the percentiles of each variable
	percentiles = np.arange(1, 99)		percentiles = np.arange(1, 99)
	@@ -484,7 +495,7 @@ def test_example_8():
	def test_example_9():		def test_example_9():
	# output directory to write output to		# output directory to write output to
	output_data_path = (		output_data_path = (
	Path(__file__).parent.parent / "resources" / "test_output" / "example9"		Path(__file__).parent / "resources" / "test_output" / "example9"
	)		)
	output_data_path.mkdir(parents=True, exist_ok=True)		output_data_path.mkdir(parents=True, exist_ok=True)

	@@ -522,25 +533,27 @@ def test_example_9():
	# clip to minimum transmissivity of 1Dm		# clip to minimum transmissivity of 1Dm
	permeability = np.clip(permeability, 1000 / thickness, 1e10)		permeability = np.clip(permeability, 1000 / thickness, 1e10)

	reservoir_properties = xr.Dataset(		porosity_da = xr.DataArray(
	{		porosity / 100, dims=["depth", "samples"],
	"thickness": thickness,		coords={"depth": depths, "samples": np.arange(n_samples)}
	"porosity": (["depth", "samples"], porosity / 100),		)
	"ntg": ntg,		permeability_da = xr.DataArray(
	"permeability": (["depth", "samples"], permeability),		permeability, dims=["depth", "samples"],
	},		coords={"depth": depths, "samples": np.arange(n_samples)}
	coords={
	"samples": np.arange(n_samples),
	"depth": depths,
	},
	)		)
			depth_da = xr.DataArray(depths, dims=["depth"], coords={"depth": depths})

	# For every sample, run a doublet simulation		aquifer = Aquifer(
	simulations = calculate_doublet_performance(		thickness=thickness,
	reservoir_properties,		porosity=porosity_da,
	chunk_size=100,		ntg=ntg,
			depth=depth_da,
			permeability=permeability_da,
	)		)

			# For every sample, run a doublet simulation
			simulations = ThermoGISDoublet(aquifer=aquifer).simulate(chunk_size=100).to_dataset()

	simulations_stoch = simulations.quantile([0.1, 0.5, 0.9], dim="samples")		simulations_stoch = simulations.quantile([0.1, 0.5, 0.9], dim="samples")

	# make a plot showing there is a sweet spot between Doublet power, Cost of energy		# make a plot showing there is a sweet spot between Doublet power, Cost of energy