update

src/fmu/tools/rms/generate_bw_per_facies.py

@@ -1,83 +1,95 @@
-import xtgeo
 import copy
+
 import numpy as np
-from xtgeo.common.constants import UNDEF
+import xtgeo
 
 DEBUG_PRINT = False
-PRJ =  project
+PRJ = project  # noqa
 GRID_MODEL = "Geogrid_Valysar"
 BW_SET = "BW_copy"
 ORIGINAL_PETRO_LOGS = ["PHIT", "KLOGH"]
-FACIES_LOG ="Facies"
+FACIES_LOG = "Facies"
 FACIES_CODE_NAMES = {
-        0: "Floodplain",
-        1: "Channel",
-        2: "Crevasse",
-        5: "Coal",
-    }
+    0: "Floodplain",
+    1: "Channel",
+    2: "Crevasse",
+    5: "Coal",
+}
+
 
 def main():
+
     """Script to be used as RMS python job to create new petrophysical
     logs from original logs, but with one log per facies.
     All grid blocks for the blocked wells not belonging to the facies
     is set to undefined.
 
-    Purpose: 
+    Purpose
+
         Create the blocked well logs to be used to condition
-        petrophysical realizations where all grid cells are 
+        petrophysical realizations where all grid cells are
         assumed to belong to only one facies. This script will
         not modify any of the original logs, only create new logs where only
         petrophysical log values for one facies is selected
         and all other are set ot undefined.
-    
-    Input:
+
+    Input
+
         grid model name,
         blocked well set name,
         list of original log names to use,
         facies log name,
         facies code with facies name dictionary
 
-    Output:
+    Output
+
         One new petro log per facies per petro variables in the
         input list of original log names. The output will be saved in
         the given blocked well set specified in the input.
+
     """
+
     create_bw_per_facies(
-        PRJ, GRID_MODEL, BW_SET,
+        PRJ,
+        GRID_MODEL,
+        BW_SET,
         ORIGINAL_PETRO_LOGS,
         FACIES_LOG,
         FACIES_CODE_NAMES,
-        debug_print= DEBUG_PRINT)
+        debug_print=DEBUG_PRINT,
+    )
 
 
 def create_bw_per_facies(
-        project,
-        grid_name,
-        bw_name,
-        original_petro_log_names,
-        facies_log_name,
-        facies_code_names,
-        debug_print=False
-    ):
+    project,
+    grid_name,
+    bw_name,
+    original_petro_log_names,
+    facies_log_name,
+    facies_code_names,
+    debug_print=False,
+):
     if debug_print:
         print(f"Petro lognames to use:  {original_petro_log_names}")
         print(f"Facies logname:  {facies_log_name}")
     original_log_names = copy.copy(original_petro_log_names)
     original_log_names.append(facies_log_name)
-    bw = xtgeo.blockedwells_from_roxar(project, grid_name, bw_name, lognames=original_log_names)
+    bw = xtgeo.blockedwells_from_roxar(
+        project, grid_name, bw_name, lognames=original_log_names
+    )
 
     for well in bw.wells:
         if debug_print:
             print(f"Wellname: {well.name}")
 
-        # Update the new logs by only keeping petro variables belonging to the current facies
+        # Update the new logs by only keeping petro variables
+        # belonging to the current facies
         df = well.get_dataframe()
         new_log_names = []
         for facies_code, fname in facies_code_names.items():
-            filtered_rows = (df[facies_log_name] != int(facies_code))
+            filtered_rows = df[facies_log_name] != int(facies_code)
             for petro_name in original_petro_log_names:
                 if petro_name in well.lognames:
-
                     new_log_name = fname + "_" + petro_name
                     well.create_log(new_log_name)
 
@@ -91,26 +103,23 @@ def create_bw_per_facies(
         if debug_print:
             print(f"Well:  {well.name}")
             print(f"All logs:  {well.lognames_all}")
-            print(f"Dataframe for facies log and new logs:")
+            print("Dataframe for facies log and new logs:")
             df_updated = well.get_dataframe()
-            selected_log_names =[]
+            selected_log_names = []
             selected_log_names.append(facies_log_name)
             selected_log_names.extend(new_log_names)
             print(f"{df_updated[selected_log_names]}")
-            print("-"*100)
+            print("-" * 100)
         print(f"Create new logs for well {well.name}: {new_log_names}")
-        well.to_roxar(project, grid_name,bw_name, well.name,
-            lognames=new_log_names, update_option= "overwrite")
-
-def main():
-    create_bw_per_facies(
-        PRJ, GRID_MODEL, BW_SET,
-        ORIGINAL_PETRO_LOGS,
-        FACIES_LOG,
-        FACIES_CODE_NAMES,
-        debug_print= DEBUG_PRINT)
+        well.to_roxar(
+            project,
+            grid_name,
+            bw_name,
+            well.name,
+            lognames=new_log_names,
+            update_option="overwrite",
+        )
 
 
 if __name__ == "__main__":
-        main()
-
+    main()

tests/rms/test_generate_bw_per_facies.py

@@ -0,0 +1,204 @@
+"""Run tests in RMS using rmsapi to new blocked well logs by using
+generate_bw_per_facies.
+
+Creates a tmp RMS project in given version.
+
+This requires a ROXAPI license, and to be ran in a "roxenvbash" environment
+
+"""
+
+import contextlib
+import shutil
+from os.path import isdir
+from pathlib import Path
+
+import numpy as np
+import pytest
+import xtgeo
+
+import roxarimport roxar.jobs
+
+with contextlib.suppress(ImportError):
+    import roxar
+
+from fmu.tools.rms.generate_bw_per_facies import create_bw_per_facies
+
+# ======================================================================================
+# settings to create RMS project!
+
+REMOVE_RMS_PROJECT_AFTER_TEST = True
+DEBUG_PRINT = True
+TMPD = Path("TMP")
+TMPD.mkdir(parents=True, exist_ok=True)
+PROJNAME = "tmp_project_generate_bw_logs.rmsxxx"
+PRJ = str(TMPD / PROJNAME)
+RESULTDIR = TMPD / "bw"
+RESULTDIR.mkdir(parents=True, exist_ok=True)
+REFERENCE_DIR = Path("tests/rms/generate_bw_testdata")
+REFERENCE_FILE = REFERENCE_DIR / "BW.txt"
+
+GRIDDATA = REFERENCE_DIR / "grid.roff"
+GRIDNAME = "GridModel"
+NWELLS = 3
+BLOCKED_WELL_SET = "BW"
+BLOCKED_WELL_JOB = "BW_job"
+ORIGINAL_PETRO_LOGS = ["PORO", "PERM"]
+FACIES_LOG = "FACIES"
+FACIES_CODE_NAMES = {
+    0: "F1",
+    1: "F2",
+    2: "F3",
+}
+
+def create_wells(project, nwells):
+    well_list = []
+    well_names = []
+    for i in range(nwells):
+        # well head and name
+        well_name = f"W{i}"
+        well_names.append(well_name)
+        well = project.wells.create(well_name)
+        well.rkb = 100
+        east = 1000 + 200*i
+        north = 2000 + 100*i
+        well.well_head = (east, north)
+
+        # trajectory
+        trajectories = well.wellbore.trajectories
+        drilled_trajectory = trajectories.create("Drilled trajectory")
+        surveypoints = drilled_trajectory.survey_point_series
+        array_with_points = surveypoints.generate_survey_points(15)
+        array_with_points[:,0] = np.arange(1000, 1200, 10)  #From MD depth 1000 to MD depth 1200
+        array_with_points[:,3].fill(east)   # vertical well (constant x,y position for all points)
+        array_with_points[:,4].fill(north)  
+        array_with_points[:,5] = np.arange(1500, 1700, 10)  #TVD
+        surveypoints.set_survey_points(array_with_points)
+
+        # log curves
+        log_run = drilled_trajectory.log_runs.create("Log run 1")
+        measured_depths = np.arange(1499.0,1701,1.0) # MD points for log values
+        log_run.set_measured_depths(measured_depths)
+        nval = len(measured_depths)
+
+        poro_log_curve = log_run.log_curves.create("PORO")
+        values = np.zeros(nval, dtype=np.float32)
+        for i in range(nval):
+            angle = i*10*np.pi/nval
+            values[i] = 0.05 * np.sin(angle) + 0.15
+        poro_log_curve.set_values(values)
+
+        perm_log_curve = log_run.log_curves.create("PERM")
+        for i in range(nval):
+            angle = i*10*np.pi/nval
+            values[i] = 99 * np.sin(angle) + 100
+        perm_log_curve.set_values(values)
+
+        facies_log_curve = log_run.log_curves.create_discrete("FACIES")
+        values = facies_log_curve.generate_values()
+        for i in range(nval):
+            angle = i*10*np.pi/nval
+            values[i] = i % 3  # (values 0,1 or 2, must be consistent with facies code names)
+        facies_log_curve.set_values(values)
+        code_names = FACIES_CODE_NAMES
+        facies_log_curve.set_code_names(code_names)
+        well_list.append(well)
+
+    return well_names
+
+def create_bw_job(owner_strings,job_type, job_name, well_names):
+    bw_job = roxar.jobs.Job.create(owner=owner_strings, type=job_type, name=job_name)
+
+    params = {
+        "BlockedWellsName": BLOCKED_WELL_SET,
+        "Continuous Blocked Log": [
+            {
+                "Name": "PORO",
+                "Interpolate": True,
+                "CellLayerAveraging": True,
+             },
+            {
+                "Name": "PERM",
+                "Interpolate": True,
+                "CellLayerAveraging": True,
+             },
+        ],
+        "Discrete Blocked Log": [
+            {
+                "Name": "FACIES",
+                "CellLayerAveraging": True, 
+            },
+        ],
+        "Wells": [["Wells", well.name] for well in well_names],
+        "Zone Blocked Log": [
+            {
+                "Name": "ZONELOG",
+                "CellLayerAveraging": True,
+                "ZoneLogArray": [1],
+            }
+        ],
+    }
+    check, err_list, warn_list = bw_job.check(params)
+    if not check:
+        print(f"Error when creating blocked well job:")
+        for i in range(len(err_list)):
+            print(f"  {err_list[i]}")
+        print(f"Warnings when creating blocked well job:")
+        for i in range(len(warn_list)):
+            print(f"  {warn_list[i]}")
+
+    bw_job.set_arguments(params)
+    bw_job.save()
+
+    return bw_job
+
+
+def create_project():
+    """Create a tmp RMS project for testing, populate with basic data.
+
+    """
+
+    prj1 = str(PRJ)
+
+    print("\n******** Setup RMS project!\n")
+    if isdir(prj1):
+        print("Remove existing project! (1)")
+        shutil.rmtree(prj1)
+
+    project = roxar.Project.create()
+
+    rox = xtgeo.RoxUtils(project)
+    print("Roxar version is", rox.roxversion)
+    print("RMS version is", rox.rmsversion(rox.roxversion))
+    assert "1." in rox.roxversion
+
+
+    # populate with grid and props
+    grd = xtgeo.grid_from_file(GRIDDATA, fformat="roff")
+    grd.to_roxar(project, GRIDNAME)
+
+    well_names = create_wells(project, NWELLS)
+    owner_strings = ["Grid Models", GRIDNAME, "Grid"]
+    bw_job = create_bw_job(owner_strings,"Block Wells", BLOCKED_WELL_JOB, well_names)
+    bw_job.execute()
+    project.save_as(prj1)
+    return project
+
+
+
+
+@pytest.mark.skipunlessroxar
+def test_generate_bw():
+    """Test generate_new blocked well logs"""
+
+    project = create_project()
+    create_bw_per_facies(
+        PRJ,
+        GRIDNAME,
+        BLOCKED_WELL_SET,
+        ORIGINAL_PETRO_LOGS,
+        FACIES_LOG,
+        FACIES_CODE_NAMES,
+        debug_print=DEBUG_PRINT,
+
+    # Export bw set før og etter og sammenlign med referanse
+