Improve scripts in RMS re. 4d_ampl_maps_dataio.py and obs_4d_ampl_maps_dataio.py

[jcrivenaes]

In general allow more flexibility wrt to non-symmetric windows for getting seismic attributes

• Extraction of seismic amplitude maps between two horizons with offsets / one horizon with asymmetric window
• Also possible to extract several windows / horizon offsets – to differentiate the maps, a tag specifying the window is added optionally
• Manual setting of interval name to give sensible names for superzones (e.g. all, upperx, lowerx) consistent with simulator map superzones

[jcrivenaes]

"""
Compute and store 4D attribute maps, and export results using fmu.dataio Attributes are
calculated with xtgeo. For a list of available attributes, see
https://xtgeo.readthedocs.io/en/latest/apireference.html#xtgeo.Cube.compute_attributes..

As long as xtgeo documentation is not updated, check source code and search for
'compute_attributes_in_window':
/prog/res/komodo/<komodo_version>/root/lib/python3.11/site-packages/xtgeo/cube/cube1.py

In this case, the attribute maps are calculated in a symmetric window around a single
horizon.

Read env variable to test for hist vs pred mode If undefined it will assume history run.
The variable is usually be set in ERT (setenv FLOWSIM_IS_PREDICTION True).

EZA/JRIV/RNYB
JRIV, may 2025, latest update
ABEC, june 2025, modified to allow for extraction windows around two horizons with offsets.
"""

from os import getenv

import numpy as np
import xtgeo
from fmu.dataio import ExportData
from fmu.tools.rms import import_localmodule

PRJ = project  # noqa

# import modules inside this RMS project; libgeneral.py and libseis.py
libgen = import_localmodule(PRJ, "libgeneral")
libseis = import_localmodule(PRJ, "libseis")

# User settings for which attribute maps to calculate
# Available calculations, see xtgeo (link in header)
ATTRNAME = "amp"  # Only used for the file name and metadata
CALCTYPES = ["rms"]  # ["max", "min", "rms"]
SCALE_FACTOR = 1.0  # Synthetic seismic cubes are multiplied by this factor
STACK_OFFSET = "far"  # Only used as metadata

# Extraction windows
# ATTR_HORIZONS_4D_WIN = [
#     {
#         "name": "all",
#         "top": "topgarn32",
#         "top_offs": -30,
#         "base": "toprorupper",
#         "base_offs": 30,
#         "addtag": False,
#     },
# ]

# Depth maps from share/results/maps of type
# DEPTHMAPS = "ds_gf_final_extracted"

# Standard settings (usually no changes below here)
# -------------------------------------------------

INITIAL_CFG = libgen.global_config("../../fmuconfig/output/global_variables.yml")

# Which horizons to calculate attribute maps for is defined in fmuconfig _seismic.yml
ATTR_HORIZONS_4D_WIN = INITIAL_CFG["global"]["seismic"]["ATTR_HORIZONS_4D_WIN"]
DEPTHMAPS = INITIAL_CFG["global"]["seismic"]["DEPTHMAP_TYPE"]

# Read environment variable FLOWSIM_IS_PREDICTION, normally ERT will set this
SIM2SEIS_PRED = getenv("FLOWSIM_IS_PREDICTION")

if SIM2SEIS_PRED:
    # i.e. if flowsim prediction mode; use other global_variables file
    print("Run in flowsim PREDICTION mode")
    CFG = libgen.global_config("../../fmuconfig/output/global_variables_pred.yml")
else:
    print("Run in flowsim HISTORICAL mode")
    CFG = INITIAL_CFG

# paths and standard naming
SEISMIC_PATH = "../../sim2seis/output/seismic_forward/"
RESULT_MAPS_PATH = "../../share/results/maps/"
SEISMIC_PREFIX = "syntseis_"
SEISMIC_POSTFIX = "_seismic_depth.segy"


def get_pred_or_hist_seisdiffdates():
    """Get correct diffdates as list of date pairs as [["YYYYMMDD","YYYYMMDD"], ...]."""

    if SIM2SEIS_PRED:
        usedates = "SEISMIC_PRED_DIFFDATES"
    else:
        usedates = "SEISMIC_HIST_DIFFDATES"

    return libseis.get_listed_seisdiffdates(CFG["global"]["dates"][usedates])


def slice_cube_4d_export():
    """Processing seismic 4D maps and export to share/result and SUMO."""

    for datepair in get_pred_or_hist_seisdiffdates():
        monitor_date, base_date = datepair
        print(f"\nDate: {datepair}")

        seismic_base_name = SEISMIC_PATH + SEISMIC_PREFIX + base_date + SEISMIC_POSTFIX
        seismic_monitor_name = (
            SEISMIC_PATH + SEISMIC_PREFIX + monitor_date + SEISMIC_POSTFIX
        )

        seismic_base = xtgeo.cube_from_file(seismic_base_name)
        seismic_monitor = xtgeo.cube_from_file(seismic_monitor_name)
        seismic_monitor.values -= seismic_base.values
        seismic_monitor.values *= SCALE_FACTOR

        for attr_dict in ATTR_HORIZONS_4D_WIN:
            name = attr_dict["name"]
            topname = attr_dict["top"]
            top_offs = int(attr_dict["top_offs"])
            basename = attr_dict["base"]
            base_offs = int(attr_dict["base_offs"])
            addtag = attr_dict["addtag"]

            # Create the window specification string used in filename and metadata
            if name == topname and name == basename:
                window_s = f"{str(top_offs)}_{str(base_offs)}"
            else:
                if top_offs == 0 and base_offs == 0:
                    window_s = f"{topname}_{basename}"
                elif top_offs == 0:
                    window_s = f"{topname}_{basename}_{str(base_offs)}"
                elif base_offs == 0:
                    window_s = f"{topname}_{str(top_offs)}_{basename}"
                else:
                    window_s = f"{topname}_{str(top_offs)}_{basename}_{str(base_offs)}"
            print(f"  Window: {window_s}")

            # Load depth map from file since this is a post modelling run!
            topsurf_file = f"{RESULT_MAPS_PATH}{topname.lower()}--{DEPTHMAPS}.gri"
            topsurf = xtgeo.surface_from_file(topsurf_file)
            if topname == basename:
                basesurf = topsurf.copy()
            else:
                basesurf_file = f"{RESULT_MAPS_PATH}{basename.lower()}--{DEPTHMAPS}.gri"
                basesurf = xtgeo.surface_from_file(basesurf_file)

            # Trying to define a 'meaningful' zrange metadata attribute
            z_med_top = np.ma.median(topsurf.values)
            z_med_base = np.ma.median(basesurf.values)
            zrange = int(z_med_base - z_med_top) + (base_offs - top_offs)
            print(f"  zrange: {z_med_top=:.1f} {z_med_base=:.1f} {top_offs=} {base_offs=} {zrange=}")

            # This new method calculates all available attributes
            attr_maps = seismic_monitor.compute_attributes_in_window(
                topsurf + top_offs, basesurf + base_offs
            )

            for calctype in CALCTYPES:
                outsurf = attr_maps[calctype]

                tagname = f"{ATTRNAME}_{calctype}"
                if addtag:
                    tagname += f"_{window_s}"

                outname = ExportData(
                    config=CFG,
                    content="seismic",
                    content_metadata={
                        "attribute": ATTRNAME,
                        "calculation": calctype,
                        "zrange": zrange,
                        "filter_size": -999.0,
                        "scaling_factor": SCALE_FACTOR,
                        "stacking_offset": STACK_OFFSET,
                    },
                    unit=None,
                    timedata=[[monitor_date, "monitor"], [base_date, "base"]],
                    workflow="rms 4D workflow",
                    name=name,
                    tagname=tagname,
                ).export(outsurf)

                print("Saved:", outname)

    print("Finished making 4d maps. Written to share folder with metadata")


def main():
    """Main is the entry point."""
    slice_cube_4d_export()


if __name__ == "__main__":
    main()

[jcrivenaes]

"""
Compute and store 4D attribute maps, and export results using fmu.dataio Attributes are
calculated with xtgeo. For a list of available attributes, see
https://xtgeo.readthedocs.io/en/latest/apireference.html#xtgeo.Cube_compute_attributes...

In this case, the attribute maps are calculated in a symmetric window around a single
horizon.

Read env variable to test for hist vs pred mode If undefined it will assume history run.
The variable is usually be set in ERT (setenv FLOWSIM_IS_PREDICTION True).

EZA/JRIV/RNYB
JRIV, may 2025, latest update
ABEC, june 2025, modified to allow for extraction windows around two horizons with offsets.
"""

import numpy as np
import xtgeo
from fmu.dataio import ExportData
from fmu.tools.rms import import_localmodule

PRJ = project  # noqa

# import modules inside this RMS project; libgeneral.py and libseis.py
libgen = import_localmodule(PRJ, "libgeneral")
libseis = import_localmodule(PRJ, "libseis")

# User settings for which attribute maps to calculate
# Available calculations, see xtgeo (link in header)
ATTRNAME = "amp"  # Only used for the file name and metadata
CALCTYPES = ["rms"]  # ["max", "min", "rms"]
SCALE_FACTOR = 1.0  # Synthetic seismic cubes are multiplied by this factor
STACK_OFFSET = "far"  # Only used as metadata

# Extraction windows
# ATTR_HORIZONS_4D_WIN = [
#     {
#         "name": "all",
#         "top": "topgarn32",
#         "top_offs": -30,
#         "base": "toprorupper",
#         "base_offs": 30,
#         "addtag": False,
#     },
# ]

# Depth maps from share/results/maps of type
# DEPTHMAPS = "ds_gf_final_extracted"

# Standard settings (usually no changes below here)
# -------------------------------------------------

CFG = libgen.global_config("../../fmuconfig/output/global_variables.yml")

# Which horizons to calculate attribute maps for is defined in fmuconfig _seismic.yml
ATTR_HORIZONS_4D_WIN = CFG["global"]["seismic"]["ATTR_HORIZONS_4D_WIN"]
DEPTHMAPS = CFG["global"]["seismic"]["DEPTHMAP_TYPE"]

SEIS_DIFFDATES = CFG["global"]["dates"]["SEISMIC_HIST_DIFFDATES"]

# paths and standard naming
SEISMIC_PATH = "../../share/preprocessed/seismic/"
RESULT_MAPS_PATH = "../../share/results/maps/"
SEISMIC_PREFIX = "seismic--amp_ts_depth--"
SEISMIC_POSTFIX = ".segy"


QUICKPLOT = False


def slice_cube_and_plot():
    """For each date and top, get attribute in cube window sampled to surface."""
    for datepair in libseis.get_listed_seisdiffdates(SEIS_DIFFDATES):
        monitor_date, base_date = datepair
        date = monitor_date + "_" + base_date

        seismic_name = SEISMIC_PATH + SEISMIC_PREFIX + date + SEISMIC_POSTFIX
        print("Reading cube:", seismic_name)
        seismic_cube = xtgeo.cube_from_file(seismic_name)

        for attr_dict in ATTR_HORIZONS_4D_WIN:
            name = attr_dict["name"]
            topname = attr_dict["top"]
            top_offs = int(attr_dict["top_offs"])
            basename = attr_dict["base"]
            base_offs = int(attr_dict["base_offs"])
            addtag = attr_dict["addtag"]

            # Create the window specification string used in filename and metadata
            if name == topname and name == basename:
                window_s = f"{str(top_offs)}_{str(base_offs)}"
            else:
                if top_offs == 0 and base_offs == 0:
                    window_s = f"{topname}_{basename}"
                elif top_offs == 0:
                    window_s = f"{topname}_{basename}_{str(base_offs)}"
                elif base_offs == 0:
                    window_s = f"{topname}_{str(top_offs)}_{basename}"
                else:
                    window_s = f"{topname}_{str(top_offs)}_{basename}_{str(base_offs)}"
            print(f"  Window: {window_s}")

            # Load depth map from file
            topsurf_file = f"{RESULT_MAPS_PATH}{topname.lower()}--{DEPTHMAPS}.gri"
            topsurf = xtgeo.surface_from_file(topsurf_file)
            if topname == basename:
                basesurf = topsurf.copy()
            else:
                basesurf_file = f"{RESULT_MAPS_PATH}{basename.lower()}--{DEPTHMAPS}.gri"
                basesurf = xtgeo.surface_from_file(basesurf_file)

            # Trying to define a 'meaningful' zrange metadata attribute
            z_med_top = np.ma.median(topsurf.values)
            z_med_base = np.ma.median(basesurf.values)
            zrange = int(z_med_base - z_med_top) + (base_offs - top_offs)
            print(f"  zrange: {z_med_top=:.1f} {z_med_base=:.1f} {top_offs=} {base_offs=} {zrange=}")

            # This new method calculates all available attributes
            attr_maps = seismic_cube.compute_attributes_in_window(
                topsurf + top_offs, basesurf + base_offs
            )

            for calctype in CALCTYPES:
                print(f"Calculation for '{ATTRNAME}': type={calctype}, window={window_s}")

                outsurf = attr_maps[calctype]

                tagname = f"{ATTRNAME}_{calctype}"
                if addtag:
                    tagname += f"_{window_s}"

                # export data, with rich metadata
                outname = ExportData(
                    config=CFG,
                    name=name,
                    tagname=tagname,
                    content="seismic",
                    rep_include=False,
                    content_metadata={
                        "attribute": ATTRNAME,
                        "calculation": calctype,
                        "zrange": zrange,
                        "stacking_offset": STACK_OFFSET,
                    },
                    is_observation=True,
                    is_prediction=False,
                    unit=None,
                    timedata=[[monitor_date, "monitor"], [base_date, "base"]],
                    preprocessed=True,
                ).export(outsurf)

                print("Saved:", outname)

                if QUICKPLOT:
                    out_img = outname.replace(".gri", ".png")
                    outsurf.quickplot(
                        out_img,
                        colormap="RdBu",
                        minmax=(-0.5, 0.5),
                        title=f"{name}--{tagname}--{date}",
                        infotext=calctype + "_" + window_s,
                    )

    print("Finished making 4d maps. Written to share folder")


def main():
    """Main routine, entry point for script."""
    slice_cube_and_plot()


if __name__ == "__main__":
    main()