Skip to content

Improve accuracy in calculation of differences between posterior and prior ensemble mean and standard deviation #855

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
oddvarlia wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Improve accuracy in calculation of differences between posterior and prior ensemble mean and standard deviation #855

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 85 additions & 36 deletions src/subscript/field_statistics/field_statistics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1339,6 +1339,11 @@ def calc_stats(
if not use_geogrid_fields:
return

# Get list of active realization (Must be active for all iterations in iter_list)
active_real, number_of_skipped = get_active_real(
iter_list, ens_path, nreal, geogrid_name
)

ensemble_path = ens_path

logger.info(f"Number of realizations: {nreal}")
Expand All @@ -1356,8 +1361,8 @@ def calc_stats(
(ertbox_size[0], ertbox_size[1], ertbox_size[2], nreal),
dtype=np.float32,
)
number_of_skipped = 0
for real_number in range(nreal):

for real_number in active_real:
grid_dimensions, subgrids, property_param = (
read_ensemble_realization(
ensemble_path,
Expand All @@ -1368,11 +1373,7 @@ def calc_stats(
geogrid_name,
)
)
if grid_dimensions is None:
txt = f" Skip non-existing realization: {real_number}"
logger.info(txt)
number_of_skipped += 1
continue
assert grid_dimensions is not None
ertbox_prop_values = get_values_in_ertbox(
grid_dimensions,
subgrids,
Expand Down Expand Up @@ -1445,8 +1446,8 @@ def calc_stats(
(ertbox_size[0], ertbox_size[1], ertbox_size[2], nreal),
dtype=np.int32,
)
number_of_skipped = 0
for real_number in range(nreal):

for real_number in active_real:
grid_dimensions, subgrids, property_param = (
read_ensemble_realization(
ensemble_path,
Expand All @@ -1457,11 +1458,7 @@ def calc_stats(
geogrid_name,
)
)

if grid_dimensions is None:
logger.info(f" Skip realization: {real_number}")
number_of_skipped += 1
continue
assert grid_dimensions is not None
ertbox_prop_values = get_values_in_ertbox(
grid_dimensions,
subgrids,
Expand Down Expand Up @@ -1576,6 +1573,8 @@ def calc_temporary_field_stats(
# Check if any need to continue to calculation
if not use_temporary_fields:
return
# Get list of active realization (Must be active for all iterations in iter_list)
active_real, number_of_skipped = get_active_real(iter_list, ens_path, nreal)

# Import realizations of temporary field parameters
for param_name in param_list:
Expand All @@ -1591,20 +1590,15 @@ def calc_temporary_field_stats(
dtype=np.float32,
)

number_of_skipped = 0
for real_number in range(nreal):
for real_number in active_real:
filepath = (
ens_path
/ Path(
"realization-" + str(real_number) + "/iter-" + str(iteration)
)
/ Path(full_param_filename)
)
if not filepath.exists():
txt = f" Skip non-existing realization: {real_number}"
logger.info(txt)
number_of_skipped += 1
continue
assert filepath.exists()
property = xtgeo.gridproperty_from_file(filepath, fformat="roff")
values = property.values
all_values[:, :, :, real_number] = values
Expand Down Expand Up @@ -1657,6 +1651,39 @@ def calc_temporary_field_stats(
xtgeo_ertbox_stdev.to_file(result_stdev_file_path, fformat="roff")


def get_active_real(
iter_list: list, ens_path: Path, nreal: int, geogrid_name: str = "NotUsed"
):
Comment on lines +1654 to +1656
Copy link
Collaborator

@alifbe alifbe Mar 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Propose to change default geogrid_name to empty string "" or None rather than arbitrary string

"""Get a list of active realizations"""
number_of_skipped = 0
active_real = []
for real_number in range(nreal):
real_exist = True
for iteration in iter_list:
ensemble_path = ens_path / Path(
"realization-" + str(real_number) + "/iter-" + str(iteration)
)
if geogrid_name != "NotUsed":
grid_path = Path("share/results/grids/" + geogrid_name + ".roff")
file_path_grid = ensemble_path / grid_path
txt = f" Skip non-existing realization of grid: {real_number}"
else:
file_path_grid = ensemble_path
txt = f" Skip non-existing realization: {real_number}"

if not file_path_grid.exists():
logger.info(txt)
number_of_skipped += 1
# No need to check other iterations since active_real should
# only be those realizations that exists for all specified
# iterations in iter_list
real_exist = False
continue
if real_exist:
active_real.append(real_number)
return active_real, number_of_skipped


def generate_script(
rms_load_script, ert_config_path, result_path, field_stat_config_file
):
Expand Down Expand Up @@ -1853,25 +1880,47 @@ def main():
if init_path and param_names:
for param_name in param_names:
for iteration in iter_list:
new_name = "mean_" + param_name + "_" + str(iteration)
param_file_name = Path(result_path) / Path(new_name + ".roff")
prop_param = xtgeo.gridproperty_from_file(
param_file_name, fformat="roff"
mean_name = "mean_" + param_name + "_" + str(iteration)
std_name = "stdev_" + param_name + "_" + str(iteration)

mean_param_file_name = Path(result_path) / Path(mean_name + ".roff")
mean_prop_param = xtgeo.gridproperty_from_file(
mean_param_file_name, fformat="roff"
)
print(f"Read: {{new_name}} into {{GRIDNAME}}")
if label:
new_name = new_name + "_" + label
prop_param.to_roxar(PRJ, GRIDNAME, new_name)
print(f"Read: {{mean_name}} into {{GRIDNAME}}")

new_name = "stdev_" + param_name + "_" + str(iteration)
param_file_name = Path(result_path) / Path(new_name + ".roff")
prop_param = xtgeo.gridproperty_from_file(
param_file_name, fformat="roff"
std_param_file_name = Path(result_path) / Path(std_name + ".roff")
std_prop_param = xtgeo.gridproperty_from_file(
std_param_file_name, fformat="roff"
)
print(f"Read: {{new_name}} into {{GRIDNAME}}")
print(f"Read: {{std_name}} into {{GRIDNAME}}")

if label:
new_name = new_name + "_" + label
prop_param.to_roxar(PRJ, GRIDNAME, new_name)
new_mean_name = mean_name + "_" + label
new_std_name = std_name + "_" + label
difference_mean_name = "diff_" + new_mean_name
difference_std_name = "diff_" + new_std_name

mean_prop_param.to_roxar(PRJ, GRIDNAME, new_mean_name)
std_prop_param.to_roxar(PRJ, GRIDNAME, new_std_name)

if iteration == iter_list[0]:
# Init
mean_prop_param_init = mean_prop_param
std_prop_param_init = std_prop_param
elif iteration == iter_list[-1]:
mean_prop_param_upd = mean_prop_param
std_prop_param_upd = std_prop_param
diff_mean_prop_param = mean_prop_param_upd.copy()
diff_std_prop_param = std_prop_param_upd.copy()
diff_mean_prop_param.values = \
mean_prop_param_upd.values - mean_prop_param_init.values
diff_std_prop_param.values = \
std_prop_param_upd.values - std_prop_param_init.values
diff_mean_prop_param.to_roxar(PRJ,
GRIDNAME, difference_mean_name)
diff_std_prop_param.to_roxar(PRJ,
GRIDNAME, difference_std_name)

if __name__ == "__main__":
main()
Expand Down