FieldStatistics: Export and copy result files to realization-0 (#744)

* Add script to calculate mean, stdev for continuous 3D field parameters and estimated facies probabilities for discrete 3D parameter

* Generated test data is written to tmp directory, added subscript doc minimum version

* Updated pyproject.toml and removed unnecessary files

* Fixed a missing copy of files

Some doc updates and some result files name changes

src/subscript/field_statistics/field_statistics.py

@@ -53,6 +53,7 @@
 
 The assumption behind this method (using ERTBOX grid as a fixed common grid for
 all realizations) is:
+
     * The lateral extension of the geogrid is close to a regular grid with same
       orientation and grid resolution as the ERTBOX grid.
     * The ERTBOX grid should be the same as used in ERT when field parameters
@@ -162,23 +163,26 @@
 
 -- Installation of the ERT workflow:
 DEFINE <FIELD_STAT_CONFIG_FILE>  ../input/config/field_param_stat.yml
-LOAD_WORKFLOW_JOB  ../../bin/jobs/WF_FIELD_PARAM_STATISTICS
+-- The workflow job FIELD_STATISTICS is generated automatically by ERT
+-- The workflow file wf_field_param_statistics run FIELD_STATISTICS
 LOAD_WORKFLOW           ../../bin/workflows/wf_field_param_statistics
 
--- The workflow file to be located under ert/bin/workflows:
+-- The workflow file to be located under ert/bin/workflows run FIELD_STATISTICS:
+-- Example of a workflow file can be
 FIELD_STATISTICS -c <FIELD_STAT_CONFIG_FILE>
                  -p <CONFIG_PATH>
                  -e <SCRATCH>/<USER>/<CASE_DIR>
                  -r <RESULT_PATH>
--- Workflow job for ERT to calculate:
---   Mean and standard deviatons of specified continuous 3D parameters.
---   Estimate of facies probabilities from discrete 3D parameter for facies.
--- The input realizations are found under:
---   <ensemble_path>/realization-*/iter-*/share/results/grids/geogrid--<name>.roff
--- The output mean and standard deviations and facies probability estimates are saved
--- under a directory specified by the user.
--- The first three command line arguments are required, the last one (<RESULT_PATH>)
--- has default 'share/grid_statistics' under <ensemble_path>.
+-- where <FIELD_STAT_CONFIG_FILE> is the usre specification for this script,
+-- and where <CONFIG_PATH> is the ERT <CONFIG_PATH> for the ERT project,
+-- and where ensemble directory is specified by  the '-e' option and
+-- where the result directory is specified by the <RESULT_PATH>. This is optional
+-- since share/grid_statistics is used as default.
+-- The results from iter-0 is also copied to 'realization-0/iter-0/share/results/grids
+-- and results from iter-3 is copied to 'realization-0/iter-3/share/results/grids.
+-- Workflow job for ERT to calculate field statistics is automatically
+-- generated by ERT from the subscript repository, but when setting it 
+-- up manually, it look like this:
 INTERNAL   False
 EXECUTABLE  ../scripts/field_statistics.py
 
@@ -233,6 +237,11 @@ def field_stat(args):
     if args.resultpath:
         relative_result_path = Path(args.resultpath)
     result_path = ens_path / relative_result_path
+
+    rms_load_script = None
+    if args.generate_rms_load_script:
+        rms_load_script = args.generate_rms_load_script
+
     glob_var_config_path = ert_config_path / Path(GLOBAL_VARIABLES_FILE)
     cfg_global = utils.yaml_load(glob_var_config_path)["global"]
     keyword = "FACIES_ZONE"
@@ -252,6 +261,12 @@ def field_stat(args):
     calc_stats(
         field_stat, ens_path, facies_per_zone, result_path, ert_config_path, ertbox_size
     )
+    ertbox_path = ert_config_path / ERTBOX_GRID_PATH
+    copy_ertbox_grid_to_result_path(ertbox_path, result_path)
+
+    if rms_load_script:
+        generate_script(rms_load_script, ert_config_path, result_path, config_file)
+
     logger.info(
         "Finished running workflow to calculate statistics "
         "for ensemble of field parameters"
@@ -308,6 +323,14 @@ def get_parser() -> argparse.ArgumentParser:
         action="version",
         version="%(prog)s (subscript version " + subscript.__version__ + ")",
     )
+
+    parser.add_argument(
+        "-z",
+        "--generate_rms_load_script",
+        type=str,
+        default="tmp_import_ensemble_field_statistics.py",
+        help=("Name of script for loading results into RMS for visualization. "),
+    )
     return parser
 
 
@@ -389,7 +412,6 @@ def get_values_in_ertbox(
         ertbox_prop_values[:, :, start_layer_ertbox:] = prop_values[
             :, :, start_layer:end_layer
         ]
-
     return ertbox_prop_values
 
 
@@ -601,6 +623,31 @@ def get_ertbox_size(ertbox_path):
     return ertbox_grid.dimensions
 
 
+def copy_ertbox_grid_to_result_path(ertbox_path, result_path):
+    if not Path(ertbox_path).exists():
+        raise IOError(f"The ertbox file does not exist in:  {ertbox_path}")
+    ertbox_grid = xtgeo.grid_from_file(ertbox_path, fformat="egrid")
+    grid_file_name = result_path / Path("ertbox.roff")
+    print(f"Copy ertbox grid file from {ertbox_path} to {grid_file_name}")
+    ertbox_grid.to_file(grid_file_name, fformat="roff")
+
+
+def copy_to_real0_dirs(field_stat, result_path, ens_path):
+    import glob
+    import shutil
+
+    iteration_list = field_stat["iterations"]
+    for iter in iteration_list:
+        source_files = result_path / Path(f"ertbox--*_{iter}.roff")
+        target_dir = ens_path / Path(f"realization-0/iter-{iter}/share/results/grids")
+        print(f"Source_files:  {source_files}")
+        print(f"Target dir: {target_dir}")
+        for f in glob.glob(source_files.as_posix()):
+            shutil.copy(f, target_dir.as_posix())
+        source_file = result_path / Path("ertbox.roff")
+        shutil.copy(source_file.as_posix(), target_dir.as_posix())
+
+
 def check_zone_conformity(zone_code_names, zone_names_used, zone_conformity):
     for zone_name, conformity in zone_conformity.items():
         if zone_name not in list(zone_code_names.values()):
@@ -887,6 +934,128 @@ def calc_stats(
                         logger.info(txt)
 
 
+def generate_script(
+    rms_load_script, ert_config_path, result_path, field_stat_config_file
+):
+    template_string = """#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from  pathlib import Path
+import xtgeo
+import yaml
+import fmu.config.utilities as utils
+
+PRJ = project
+
+GRIDNAME = "ERTBOX"
+
+ERT_CONFIG_PATH = "{ert_config_path}"
+
+GLOBAL_VARIABLES_FILE  = \
+    ERT_CONFIG_PATH / Path("../../fmuconfig/output/global_variables.yml")
+
+FIELD_STAT_CONFIG_FILE = "{field_stat_config_file}"
+
+RESULT_PATH = "{result_path}"
+
+
+LABEL = "drogon"
+
+def read_field_stat_config(config_file_name):
+    print(f"Read file: {{config_file_name}}")
+    with open(config_file_name, encoding="utf-8") as yml_file:
+        return yaml.safe_load(yml_file)
+
+def get_facies_per_zone(glob_var_file):
+    cfg_global = utils.yaml_load(glob_var_file)["global"]
+    keyword = "FACIES_ZONE"
+    if keyword in cfg_global:
+        facies_per_zone = cfg_global[keyword]
+    else:
+        raise KeyError(f"Missing keyword: {{keyword}} in {{GLOBAL_VARIABLES_FILE}}")
+    return facies_per_zone
+
+def main():
+    config_dict = read_field_stat_config(FIELD_STAT_CONFIG_FILE)
+    field_stat = config_dict["field_stat"]
+    zone_code_names = field_stat["zone_code_names"]
+    facies_per_zone = get_facies_per_zone(GLOBAL_VARIABLES_FILE)
+    result_path = RESULT_PATH
+    zone_list= list(zone_code_names.values())
+    stat_list= ["mean", "stdev"]
+    iter_list = field_stat["iterations"]
+
+    cont_prop_dict = field_stat["continuous_property_param_per_zone"]
+
+    discrete_prop_dict = field_stat["discrete_property_param_per_zone"]
+
+    label = LABEL
+    for stat in stat_list:
+        for zone in zone_list:
+            for iteration in iter_list:
+                if cont_prop_dict:
+                    if zone in cont_prop_dict:
+                        for prop_name in cont_prop_dict[zone]:
+                            name = \
+                            "ertbox--" + stat + "_" + zone + "_" + prop_name \
+                            + "_" + str(iteration)
+                            print(f"Read: {{name}} into {{GRIDNAME}}")
+                            filename = Path(result_path) / Path(name + ".roff")
+                            prop_param = xtgeo.gridproperty_from_file(
+                                filename,
+                                fformat="roff")
+                            prop_param.to_roxar(PRJ, GRIDNAME, name)
+                            new_name = name
+                            if label:
+                                new_name = name + "_" + label
+                            prop_param.name = new_name
+                        name = "ertbox--nactive_" + zone + "_" + str(iteration)
+                        print(f"Read: {{name}} into {{GRIDNAME}}")
+                        filename =  Path(result_path) / Path(name + ".roff")
+                        prop_param = xtgeo.gridproperty_from_file(
+                            filename,
+                            fformat="roff"
+                        )
+                        new_name = name
+                        if label:
+                            new_name = name + "_" + label
+                        prop_param.name = new_name
+                        prop_param.to_roxar(PRJ, GRIDNAME, new_name)
+                if discrete_prop_dict:
+                    code_names_per_zone = facies_per_zone[zone]
+                    for _, fname in code_names_per_zone.items():
+                        name = \
+                        "ertbox--prob_" + zone + "_" + fname + "_" + str(iteration)
+                        print(f"Read: {{name}} into {{GRIDNAME}}")
+                        filename = Path(result_path) / Path(name + ".roff")
+                        prop_param = \
+                            xtgeo.gridproperty_from_file(filename, fformat="roff")
+                        new_name = name
+                        if label:
+                            new_name = name + "_" + label
+                        prop_param.name = new_name
+                        prop_param.to_roxar(PRJ, GRIDNAME, new_name)
+                    name = "ertbox--nactive_" + zone + "_" + str(iteration)
+                    print(f"Read: {{name}} into {{GRIDNAME}}")
+                    filename =  Path(result_path) / Path(name + ".roff")
+                    prop_param = xtgeo.gridproperty_from_file(filename, fformat="roff")
+                    prop_param.to_roxar(PRJ, GRIDNAME, name)
+
+if __name__ == "__main__":
+    main()
+"""
+    print(f"Write file: {rms_load_script}")
+    with open(rms_load_script, "w") as file:
+        file.write(
+            template_string.format(
+                ert_config_path=ert_config_path,
+                field_stat_config_file=field_stat_config_file,
+                result_path=result_path,
+            )
+        )
+        file.write("\n")
+
+
 class FieldStatistics(ErtScript):
     """This class defines the ERT workflow hook.
 

tests/test_field_statistics.py

@@ -28,7 +28,7 @@
 ERT_CONFIG_PATH = Path("ert/model")
 DATADIR = Path(__file__).absolute().parent / TESTDATA
 GLOBAL_VARIABLES_FILE = Path("../../fmuconfig/output/global_variables.yml")
-
+RMS_LOAD_SCRIPT_NAME = "tmp_import_field_stat_results.py"
 
 CONFIG_DICT = {
     "nreal": 10,
@@ -64,7 +64,7 @@ def make_box_grid(dimensions, grid_name, result_path):
     filename_egrid = result_path / Path(grid_name.upper() + ".EGRID")
 
     grid = xtgeo.create_box_grid(dimensions)
-    grid.name = grid_name
+    grid.name = grid_name.lower()
     print(f"Grid name:  {grid.name}")
     print(f"Grid dimensions: {grid.dimensions}")
     print(f"Write grid to file:  {filename}")
@@ -232,6 +232,7 @@ def make_ensemble_test_data(
                         xtgeo_geogrid.set_actnum(xtgeo_active)
                         set_subgrid_names(xtgeo_geogrid, new_subgrids=subgrid_dict)
                         xtgeo_geogrid.to_file(filename_grid, fformat="roff")
+
             if print_info:
                 print(
                     "Testdata for ensemble for zone "
@@ -431,7 +432,6 @@ def compare_with_referencedata(ens_path, result_path, print_check=False):
         if words[0] in ["mean", "stdev", "prob"]:
             fullfilename = result_path / Path("ertbox--" + name)
             reference_filename = result_path / Path("referencedata") / Path(name)
-
             grid_property = xtgeo.gridproperty_from_file(fullfilename, fformat="roff")
             grid_property_reference = xtgeo.gridproperty_from_file(
                 reference_filename, fformat="roff"
@@ -891,6 +891,7 @@ def test_main(tmp_path, config_file, config_dict, print_info=True):
     ens_path = tmp_testdata_path / ENSEMBLE
     result_path = ens_path / RESULT_PATH
 
+    rms_load_script = result_path / RMS_LOAD_SCRIPT_NAME
     # Run the main script as a subprocess
     script_name = Path(__file__).absolute().parent.parent / Path(
         "src/subscript/field_statistics/field_statistics.py"
@@ -911,6 +912,8 @@ def test_main(tmp_path, config_file, config_dict, print_info=True):
             ens_path.as_posix(),
             "-r",
             result_path.as_posix(),
+            "-z",
+            rms_load_script.as_posix(),
         ]
     )
     # For this test not to fail, the CONFIG_DICT and the specified