Examples inside RMS

RMS is a licensed proprietary modeling software developed by AspenTech. From version 10 it has its own python engine integrated, and XTGeo is designed to work inside this environment. The integration will be continuously improved.

Hence XTGeo can read most datatypes that are exposed in RMS’ API (called ROXAPI), and then all native methods in XTGeo can be applied on those data. For example, if you want to write a surface from RMS to a format that ROXAPI does not support, but XTGeo supports, then it is quite easy. XTGeo can also read data from external files and store the data in the RMS data tree.

Note that all these script examples are assumed to be ran inside a python job within RMS.

Get and set data

In general, data are imported into XTGeo by a from_roxar() or by a xtgeo.xxx_from_roxar() (where xx is “surface”, “grid”, etc). Then the altered instance can be stored in roxar/RMS by a to_roxar() method.

The to_roxar() method will not do a project save when inside RMS or when inside a virtual project setting. However, if a project is applied as a file path, then to_roxar will save implicitly. Examples:

Inside RMS GUI

import xtgeo

surf = xtgeo.surface_from_roxar(project, "TopReek", "DS_extracted")
surf.values += 100
surf.to_roxar(project)

# Note: project save needs to be done by user (GUI action)

Outside RMS, direct access

import xtgeo

myproject = "/some/file/path/reek.rms11.1.1"

surf = xtgeo.surface_from_roxar(myproject, "TopReek", "DS_extracted")
surf.values += 100
surf.to_roxar(myproject)

# Note: project save is done automatic

Outside RMS, project mode

import xtgeo

myproject = "/some/file/path/reek.rms11.1.1"
rox = xtgeo.RoxUtils(myproject)

surf = xtgeo.surface_from_roxar(rox.project, "TopReek", "DS_extracted")
surf.values += 100
surf.to_roxar(rox.project)

# Note: project save is not done automatic, you need to:

rox.project.save()
rox.project.close()

Surface data

Here are some simple examples on how to use XTGeo to interact with RMS data, and e.g. do quick exports to files.

Export a surface in RMS to irap binary format

import xtgeo

# import (transfer) data from RMS to XTGeo and export
surf = xtgeo.surface_from_roxar(project, "TopReek", "DS_extracted")

surf.to_file("topreek.gri")

# modify surface, add 1000 to all map nodes
surf.values += 1000

# store in RMS (category must exist)
surf.to_roxar(project, "TopReek", "DS_whatever")

Export a surface in RMS to zmap ascii format

Note here that an automatic resampling to a nonrotated regular grid will be done in case the RMS map has a rotation.

import xtgeo as xt

# surface names
hnames = ["TopReek", "MiddleReek", "LowerReek"]

# loop over stratigraphy
for name in hnames:
    surf = xt.surface_from_roxar(project, name, "DS_extracted")
    fname = name.lower()  # lower case file name
    surf.to_file(fname + ".zmap", fformat="zmap_ascii")

print("Export done")

Take a surface in RMS and multiply values with 2:

import xtgeo

surf = xtgeo.surface_from_roxar(project, "TopReek", "DS_tmp")

surf.values *= 2  # values is the masked 2D numpy array property

# store the surface back to RMS
surf.to_roxar(project, "TopReek", "DS_tmp")

Do operations on surfaces, also inside polygons:

Find the diff maps in time domain, of the main surfaces. Also make a a version where cut by polygons where surfaces has interp (minimum common multiplum)

import xtgeo
from fmu.config import utilities as ut

CFG = ut.yaml_load("../../fmuconfig/output/global_variables.yml")["rms"]

# ========= SETTINGS ===================================================================

PRJ = project  # noqa

# input
TSCAT1 = "TS_interp_raw_ow"
PCAT = "TL_interp_raw_approx_outline"


# output
ISCAT1 = "IS_twt_main_interp_raw_ow"
ISCAT2 = "IS_twt_main_interp_raw_ow_cut"

# ========= END SETTINGS ===============================================================


def main():

    topmainzones = CFG["horizons"]["TOP_MAINRES"]
    mainzones = CFG["zones"]["MAIN_ZONES"]
    for znum, mzone in enumerate(mainzones):

        surf1 = xtgeo.surface_from_roxar(PRJ, topmainzones[znum], TSCAT1)
        surf2 = xtgeo.surface_from_roxar(PRJ, topmainzones[znum + 1], TSCAT1)

        diff = surf2.copy()
        diff.values -= surf1.values
        diff.to_roxar(PRJ, mzone, ISCAT1, stype="zones")
        print("Store {} at {}".format(mzone, ISCAT1))

        # extract differences inside a polygon and compute min/max values:

        poly = xtgeo.polygons_from_roxar(PRJ, topmainzones[znum], PCAT)
        surf1.eli_outside(poly)
        surf2.eli_outside(poly)
        diff2 = surf2.copy()
        diff2.values -= surf1.values
        print(
           "Min and max values inside polygons {} : {} (negative OK) for {}".format(
                 diff2.values.min(), diff2.values.max(), mzone
                 )
             )
        diff2.to_roxar(PRJ, mzone, ISCAT2, stype="zones")
        print("Store cut surface {} at {}".format(mzone, ISCAT2))


 if __name__ == "__main__":
     main()
     print("Done, see <{}> and <{}>".format(ISCAT1, ISCAT2))

3D grid data

Exporting geometry to ROFF file

import xtgeo

# import (transfer) data from RMS to XTGeo and export
mygrid = xtgeo.grid_from_roxar(project, "Geomodel")

mygrid.to_file("topreek.roff")  # roff binary is default format

Edit a porosity in a 3D grid

import xtgeo

# import (transfer) data from RMS to XTGeo
myporo = xtgeo.gridproperty_from_roxar(project, "Geomodel", "Por")

# now I want to limit porosity to 0.35 for values above 0.35:

myporo.values[myporo.values > 0.35] = 0.35

# store to another icon
poro.to_roxar(project, "Geomodel", "PorNew")

Edit a permeability given a porosity cutoff

import numpy as np
import xtgeo

myporo = xtgeo.gridproperty_from_roxar(project, "Geomodel", "Por")
myperm = xtgeo.gridproperty_from_roxar(project, "Geomodel", "Perm")

# if poro < 0.01 then perm is 0.001, otherwise keep as is, illustrated with np.where()
myperm.values = np.where(myporo.values < 0.1, 0.001, myperm.values)

# store to another icon
myperm.to_roxar(project, "Geomodel", "PermEdit")

Edit a 3D grid porosity inside polygons

# Example where I want to read a 3D grid porosity, and set value
# to 99 inside polygons

import xtgeo

mygrid = xtgeo.grid_from_roxar(project, "Reek_sim")
myprop = xtgeo.gridproperty_from_roxar(project, "Reek_sim", "PORO")

# read polygon(s), from Horizons, Faults, Zones or Clipboard
mypoly = xtgeo.polygons_from_roxar(project, "TopUpperReek", "DL_test")

# need to connect property to grid geometry when using polygons
myprop.geometry = mygrid

myprop.set_inside(mypoly, 99)

# Save in RMS as a new icon
myprop.to_roxar(project, "Reek_sim", "NEWPORO_setinside")

Create region polygons from the grid

import numpy as np
import xtgeo

GNAME = "Simgrid"
REGNAME = "Regions"
ZONENAME = "Zone"

CB_FOLDER = "Region_polygons"

ZONE_FILTER = [2, 3]

# factor that controls the precision of the polygons
# higher value gives smoother more convex polygons.
ALPHA_FACTOR = 1

def create_region_polygons():
    """Create region polygons and store them on the clipboard"""
    grid = xtgeo.grid_from_roxar(project, GNAME)
    reg = xtgeo.gridproperty_from_roxar(project, GNAME, REGNAME)
    zone = xtgeo.gridproperty_from_roxar(project, GNAME, ZONENAME)

    for regnum, regname in reg.codes.items():
        print(f"Creating boundary polygon for region {regname}")

        # create a filter array to extract boundaries for the region
        # zone was used as filter to minimice overlap of the final polygons
        # Note: a layer filter could have been applied instead
        filter_array = (reg.values==regnum) & (np.isin(zone.values, [ZONE_FILTER]))

        pol = grid.get_boundary_polygons(ALPHA_FACTOR, filter_array=filter_array)

        # in case of several polygons, keep only the largest (first)
        pol.filter_byid([0])

        # store polygon to the clipboard
        pol.to_roxar(project, regname, CB_FOLDER, stype="clipboard")

     print(f"Complete, region polygons are stored under clipboard folder {CB_FOLDER}")

 if __name__ == "__main__":
     create_region_polygons()

Make a hybrid grid

XTGeo can convert a conventional grid to a so-called hybrid-grid where a certain depth interval has horizontal layers.

import xtgeo

PRJ = project  # noqa
GNAME_INPUT = "Mothergrid"
GNAME_HYBRID = "Simgrid"
REGNAME = "Region"
HREGNAME = "Hregion"

NHDIV = 22
REGNO = 1
TOP = 1536
BASE = 1580


def hregion():
    """Make a custom region property for hybrid grid"""
    tgrid = xtgeo.grid_from_roxar(PRJ, GNAME_INPUT)
    reg = xtgeo.gridproperty_from_roxar(PRJ, GNAME_INPUT, REGNAME)

    reg.values[:, :, :] = 1
    reg.values[:, 193:, :] = 0  # remember 0 base in NP arrays

    reg.to_roxar(PRJ, GNAME_INPUT, HREGNAME)  # store for info/check

    return tgrid, reg


def make_hybrid(grd, reg):
    """Convert to hybrid and store in RMS project"""
    grd.convert_to_hybrid(nhdiv=NHDIV, toplevel=TOP, bottomlevel=BASE, region=reg,
                          region_number=1)

    grd.inactivate_by_dz(0.001)
    grd.to_roxar(PRJ, GNAME_HYBRID)


if __name__ == "__main__":

    print("Make hybrid...")
    grd, reg = hregion()
    make_hybrid(grd, reg)
    print("Make hybrid... done!")

Cube data

Slicing a surface in a cube

Examples to come…

Well data

Get average properties per zone

import xtgeo

PRJ = project  # noqa
WELLNAME = "DC1-1V4_ref"
TRAJNAME = "Imported trajectory"
ZONELOGNAME = "ZONELOG"
ZNAMES = {0: "UPPER", 1: "MIDDLE", 2: "LOWER"}


def get_well():
    """Get XTGeo Well() object"""
    wll = xtgeo.well_from_roxar(PRJ, WELLNAME, trajectory=TRAJNAME)
    return wll


def compute_avg_per_zone(wll):
    """Compute avg per zone without any other criteria"""

    df = wll.get_dataframe()
    df_avgs = df.groupby(ZONELOGNAME).mean()
    df_avgs.rename(index=ZNAMES, inplace=True)  # rename zonelog numbers with true name

    print("Average properties per zone")
        print(df_avgs)

        # e.g. get avg PORO for MIDDLE, rounded to 3 decimals:
        print("\nAVG poro for MIDDLE is {:2.3f}\n".format(df_avgs.loc["MIDDLE", "PORO"]))


def compute_avg_per_zone_smarter(wll):
    """Compute avg per zone by looking only at intervals that increase"""

    wll.zonelogname = ZONELOGNAME
    wll.make_zone_qual_log("QUAL")

    # This quality log will be 1 if zonelog is truly increasing, or 2 if truly
    # decreasing, so here I will only here filter on increasing (downward)
    # cf: https://xtgeo.readthedocs.io/en/latest/apiref/xtgeo.well.well1.html#
    # xtgeo.well.well1.Well.make_zone_qual_log

    df = wll.get_dataframe()
    df = df[df.QUAL == 1]  # only get the increasing part
    df_avgs = df.groupby(ZONELOGNAME).mean()
    df_avgs.rename(index=ZNAMES, inplace=True)  # rename zonelog numbers with name

    print("\n\nAverage properties per zone where penetrating zone downwards")
    print(df_avgs)


def main():
    mywell = get_well()
    compute_avg_per_zone(mywell)
    compute_avg_per_zone_smarter(mywell)


if __name__ == "__main__":
    main()

Filter logs on facies/zone boundaries

Petrophysical property modelling can be more precise if so-called shoulder effects are filtered. Here is a small example on how to do this:

import xtgeo

PRJ = project

TRAJNAME = "Drilled trajectory"
LRUNNAME = "log"
ZONELOGNAME = "Zone"
FACIESLOGNAME = "Facies"
INLOGS = [ZONELOGNAME, FACIESLOGNAME]
PETROLOGS = {"KLOGH": "KLOGH_orig", "PHIT": "PHIT_orig", "Sw": "Sw_orig"}
FILTER: {"tvd": 1.5}  # filter 1.5m below and above boundary in TVD


def filter_shoulder():
    """Filter shoulder bed data."""
    for rms_well in PRJ.wells:
        wll = xtgeo.well_from_roxar(
            PRJ, rms_well.name, trajectory=TRAJNAME, logrun=LRUNNAME
        )

        # skip wells without facies
        if FACIESLOGNAME not in wll.get_dataframe() or not rms_well.name.startswith("55"):
            continue

        print("Use: ", rms_well.name)

        # keep the original logs and work on copy:
        well_df = wll.get_dataframe()
        for target, orig in PETROLOGS.items():
            if target in well_df.columns:
                if orig not in well_df.columns:
                    # first time; create an "_orig" column
                    print("Create", orig)
                    wll.create_log(orig)
                    dataframe = wll.get_dataframe()
                    dataframe[orig] = dataframe[target].copy()
                    wll.set_dataframe(dataframe)

        uselogs = list(PETROLOGS.keys())

        wll.mask_shoulderbeds(inputlogs=INLOGS, targetlogs=uselogs, nsamples=2)
        wll.to_roxar(PRJ, rms_well.name, trajectory=TRAJNAME, logrun=LRUNNAME)


if __name__ == "__main__":
    filter_shoulder()

Blocked well data

Remember that RMS define blocked wells as a special grid property while XTGeo treats blocked wells as a subclass of Well() data.

Make new blocked logs from facies

In the following example, the blocked facies is used to make new logs that will be input to Equinor’s APS module.

import numpy as np
import xtgeo

PRJ = project
GNAME = "Geogrid_Valysar"
BWNAME = "BW"
FACIES = "Facies"

APS_FACIES = {0: "Floodplain", 1: "Channel", 2: "Crevasse", 5: "Coal"}
PREFIX = "aps_"

# note it is possible to "play with" probabilities that are not just 0 or 1
MINPROB = 0.0
MAXPROB = 1.0

def main():
    """Main work, looping wells and make APS relevant logs"""

    for well in PRJ.wells:

        blw = xtgeo.blockedwell_from_roxar(
            PRJ, GNAME, BWNAME, well.name, lognames=[FACIES]
        )
        dfr = blw.get_dataframe()
        for code, faciesname in APS_FACIES.items():
            newname = PREFIX + faciesname
            dfr[newname] = MINPROB
            dfr[newname][dfr[FACIES] == code] = MAXPROB

            # if facies is undefined, also probability shall be undefined
            dfr[newname][np.isnan(dfr[FACIES])] = np.nan

        blw.set_dataframe(dfr)
        blw.to_roxar(PRJ, GNAME, BWNAME, well.name)


if __name__ == "__main__":
    main()

Line point data

Add to or remove points inside or outside polygons

In the following example, remove or add to points being inside or outside polygons on clipboard.

import xtgeo

PRJ = project

POLYGONS = ["mypolygons", "myfolder"]  # mypolygons in folder myfolder on clipboard
POINTSET1 = ["points1", "myfolder"]
POINTSET2 = ["points2", "myfolder"]

POINTSET1_UPDATED = ["points1_edit", "myfolder"]
POINTSET2_UPDATED = ["points2_edit", "myfolder"]

def main():
    """Operations on points inside or outside polygons."""

    poly = xtgeo.polygons_from_roxar(PRJ, *POLYGONS, stype="clipboard")
    po1 = xtgeo.points_from_roxar(PRJ, *POINTSET1, stype="clipboard")
    po2 = xtgeo.points_from_roxar(PRJ, *POINTSET2, stype="clipboard")

    po1.eli_inside_polygons(poly)
    po1.to_roxar(PRJ, *POINTSET1_UPDATED, stype="clipboard")  # store

    # now add 100 inside polugons for POINTSET2, and then remove all points outside
    po2.add_inside_polygons(poly, 100)
    po2.eli_outside_polygons(poly)
    po2.to_roxar(PRJ, *POINTSET2_UPDATED, stype="clipboard")  # store


if __name__ == "__main__":
    main()