Full sim2seis workflow¶

In [1]:

            
                Copied!
                
%load_ext autoreload
%autoreload 2
%load_ext autoreload
%autoreload 2

In [2]:

            
                Copied!
                
import matplotlib.pyplot as plt
from matplotlib import rc
import numpy as np
import pandas as pd
import xarray as xr

rc("font", size=15)
import matplotlib.pyplot as plt
from matplotlib import rc
import numpy as np
import pandas as pd
import xarray as xr

rc("font", size=15)

Fetch the simulation data¶

In [3]:

            
                Copied!
                
from sim2x.datasets import fetch_t1a_example_data
from sim2x.datasets import fetch_t1a_example_data

In [4]:

            
                Copied!
                
t1a = fetch_t1a_example_data()["t1a.zip"]
t1a = fetch_t1a_example_data()["t1a.zip"]

Extract the Simulation Data¶

In [5]:

            
                Copied!
                
from eclx import EclDeck

sim = EclDeck()
sim.load_grid(t1a / "TUT1A.EGRID")
sim.load_init(t1a / "TUT1A.INIT")
sim.set_rst(t1a / "TUT1A.UNRST")
sim.load_rst(reports=[0, 5, 10], keys=["PRESSURE", "SWAT", "BO"])
from eclx import EclDeck

sim = EclDeck()
sim.load_grid(t1a / "TUT1A.EGRID")
sim.load_init(t1a / "TUT1A.INIT")
sim.set_rst(t1a / "TUT1A.UNRST")
sim.load_rst(reports=[0, 5, 10], keys=["PRESSURE", "SWAT", "BO"])

Loading KW: FIPNUM: 100%|██████████| 23/23 [00:00<00:00, 571.55it/s]
Loading RST Report 0:   0%|          | 0/3 [00:00<?, ?it/s]
  0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: SWAT:   0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: PRESSURE: 100%|██████████| 2/2 [00:00<00:00, 254.22it/s]
Loading RST Report 5:  33%|███▎      | 1/3 [00:00<00:00, 33.94it/s]
  0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: SWAT:   0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: PRESSURE: 100%|██████████| 2/2 [00:00<00:00, 283.74it/s]
Loading RST Report 10:  67%|██████▋   | 2/3 [00:00<00:00, 34.02it/s]
  0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: SWAT:   0%|          | 0/2 [00:00<?, ?it/s]
Loading KW: PRESSURE: 100%|██████████| 2/2 [00:00<00:00, 294.44it/s]
Loading RST Report 10: 100%|██████████| 3/3 [00:00<00:00, 34.97it/s]

Run sim2imp¶

In [6]:

            
                Copied!
                
                    
                    
                
                

        
from digirock import (
    WaterECL,
    WoodsFluid,
    Mineral,
    VRHAvg,
    NurCriticalPoroAdjust,
    GassmannRock,
    DeadOil,
    Transform,
)
from digirock.utils.ecl import EclStandardConditions
from digirock.fluids.bw92 import wat_salinity_brine
from digirock.typing import NDArrayOrFloat, NDArrayOrInt

sal = wat_salinity_brine(
    EclStandardConditions.TEMP.value, EclStandardConditions.PRES.value, 1.00916
)
wat = WaterECL(31.02641, 1.02, 3e-6, 0.8, sal)
oil = DeadOil(std_density=0.784905)
oil.set_pvt([1.25, 1.2, 1.15], [2.06843, 5.51581, 41.36854])
fl = WoodsFluid(["swat", "soil"], [wat, oil])

sand = Mineral(2.75, 32, 14)
clay = Mineral(2.55, 25, 8)
vrh = VRHAvg(["vsand", "vclay"], [sand, clay])
ncp = NurCriticalPoroAdjust(["poro"], vrh, 0.39)

grock = GassmannRock(ncp, vrh, fl)

# we need to create a VSH transform to turn NTG into VSH
class VShaleMult(Transform):

    _methods = ["vp", "vs", "density", "bulk_modulus", "shear_modulus"]

    def __init__(self, element, mult):
        super().__init__(["NTG"], element, self._methods)
        self._mult = mult

    def density(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.density(props, **kwargs)

    def vp(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.vp(props, **kwargs)

    def vs(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.vs(props, **kwargs)

    def bulk_modulus(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.bulk_modulus(props, **kwargs)

    def shear_modulus(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.shear_modulus(props, **kwargs)


rock = VShaleMult(grock, 0.5)
from digirock import (
    WaterECL,
    WoodsFluid,
    Mineral,
    VRHAvg,
    NurCriticalPoroAdjust,
    GassmannRock,
    DeadOil,
    Transform,
)
from digirock.utils.ecl import EclStandardConditions
from digirock.fluids.bw92 import wat_salinity_brine
from digirock.typing import NDArrayOrFloat, NDArrayOrInt

sal = wat_salinity_brine(
    EclStandardConditions.TEMP.value, EclStandardConditions.PRES.value, 1.00916
)
wat = WaterECL(31.02641, 1.02, 3e-6, 0.8, sal)
oil = DeadOil(std_density=0.784905)
oil.set_pvt([1.25, 1.2, 1.15], [2.06843, 5.51581, 41.36854])
fl = WoodsFluid(["swat", "soil"], [wat, oil])

sand = Mineral(2.75, 32, 14)
clay = Mineral(2.55, 25, 8)
vrh = VRHAvg(["vsand", "vclay"], [sand, clay])
ncp = NurCriticalPoroAdjust(["poro"], vrh, 0.39)

grock = GassmannRock(ncp, vrh, fl)

# we need to create a VSH transform to turn NTG into VSH
class VShaleMult(Transform):

    _methods = ["vp", "vs", "density", "bulk_modulus", "shear_modulus"]

    def __init__(self, element, mult):
        super().__init__(["NTG"], element, self._methods)
        self._mult = mult

    def density(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.density(props, **kwargs)

    def vp(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.vp(props, **kwargs)

    def vs(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.vs(props, **kwargs)

    def bulk_modulus(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.bulk_modulus(props, **kwargs)

    def shear_modulus(self, props, **kwargs):
        props["vclay"] = props["NTG"] * self._mult
        return self.element.shear_modulus(props, **kwargs)


rock = VShaleMult(grock, 0.5)

In [7]:

            
                Copied!
                
from sim2x import transform_units, sim2imp
from sim2x import transform_units, sim2imp

In [8]:

            
                Copied!
                
sim.data
sim.data

Out[8]:

	i	j	k	active	actnum	centerx	centery	centerz	DY	MINPVV	...	PORV	MULTY	EQLNUM	FIPNUM	SWAT_0	PRESSURE_0	SWAT_5	PRESSURE_5	SWAT_10	PRESSURE_10
0	0	0	0	0	1	250.0	250.0	8025.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4485.393555	0.256655	5227.921387	0.329152	5289.645508
1	1	0	0	1	1	750.0	250.0	8025.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4485.393555	0.261818	5378.076172	0.336713	5427.173828
2	2	0	0	2	1	1250.0	250.0	8025.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4485.393555	0.264957	5469.159180	0.336434	5507.694336
3	3	0	0	3	1	1750.0	250.0	8025.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4485.393555	0.269414	5524.484375	0.337670	5555.153809
4	4	0	0	4	1	2250.0	250.0	8025.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4485.393555	0.275811	5551.415527	0.342052	5577.749512
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
70	0	4	2	70	1	250.0	2250.0	8125.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4514.608398	0.515489	5633.163574	0.712018	5655.128418
71	1	4	2	71	1	750.0	2250.0	8125.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4514.608398	0.619617	5657.175781	0.756979	5673.387695
72	2	4	2	72	1	1250.0	2250.0	8125.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4514.608398	0.683569	5699.605469	0.778044	5709.502930
73	3	4	2	73	1	1750.0	2250.0	8125.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4514.608398	0.739840	5762.413086	0.792490	5768.013672
74	4	4	2	74	1	2250.0	2250.0	8125.0	500.0	0.000001	...	445269.0	1.0	1	1	0.25	4514.608398	0.784475	5861.002441	0.799081	5864.878418

75 rows × 37 columns

In [9]:

            
                Copied!
                
                    
                    
                
                

        
# convert sim untis from field to sim2x compatible
transform_units(sim)
imps = sim2imp(
    sim,  # the EclDeck
    rock,  # the full rock model
    {
        "PRESSURE": "pres",
        "SWAT": "swat",
        "PORO": "poro",
    },  # simulator props in sim.data mapped to rock model props
    sim.loaded_reports,  # the reports to run sim2imp on
    extra_props={"temp": 100},  # any extra props needed by `rock`
)
imps.keys()
# convert sim untis from field to sim2x compatible
transform_units(sim)
imps = sim2imp(
    sim,  # the EclDeck
    rock,  # the full rock model
    {
        "PRESSURE": "pres",
        "SWAT": "swat",
        "PORO": "poro",
    },  # simulator props in sim.data mapped to rock model props
    sim.loaded_reports,  # the reports to run sim2imp on
    extra_props={"temp": 100},  # any extra props needed by `rock`
)
imps.keys()

Out[9]:

dict_keys([0, 5, 10])

Extract a regular grid GI index¶

In [10]:

            
                Copied!
                
from sim2x import cpgrid_to_rg
from segysak import create2d_dataset, create3d_dataset
from sim2x import cpgrid_to_rg
from segysak import create2d_dataset, create3d_dataset

Create a 2D data volume with cdp geometry we want to model with sim2seis

In [11]:

            
                Copied!
                
                    
                    
                
                

        
arb_line = create2d_dataset((25, 1))
arb_line["cdp_x"] = (("cdp",), np.linspace(-50, 600, 25))
arb_line["cdp_y"] = (("cdp",), np.linspace(600, -50, 25))
arb_line = arb_line.set_coords(
    ("cdp_x", "cdp_y")
)  # make coordinates so they get inheritted by child volumes (e.g. gi_vol)
# arb_line["depth"] = (("depth", ), np.linspace(7980, 8220, 200))
# arb_line.expand_dims({"depth":np.linspace(7980, 8220, 200)})
arb_line = create2d_dataset((25, 1))
arb_line["cdp_x"] = (("cdp",), np.linspace(-50, 600, 25))
arb_line["cdp_y"] = (("cdp",), np.linspace(600, -50, 25))
arb_line = arb_line.set_coords(
    ("cdp_x", "cdp_y")
)  # make coordinates so they get inheritted by child volumes (e.g. gi_vol)
# arb_line["depth"] = (("depth", ), np.linspace(7980, 8220, 200))
# arb_line.expand_dims({"depth":np.linspace(7980, 8220, 200)})

In [14]:

            
                Copied!
                
plt.scatter(vol_3d.cdp_x, vol_3d.cdp_y, label="vol3d")
plt.scatter(arb_line.cdp_x, arb_line.cdp_y, label="arb_line")
plt.scatter(sim.xyzcorn["x0"], sim.xyzcorn["y0"], label="cp")
plt.legend()
plt.scatter(vol_3d.cdp_x, vol_3d.cdp_y, label="vol3d")
plt.scatter(arb_line.cdp_x, arb_line.cdp_y, label="arb_line")
plt.scatter(sim.xyzcorn["x0"], sim.xyzcorn["y0"], label="cp")
plt.legend()

Out[14]:

<matplotlib.legend.Legend at 0x7f7205072b20>

Use the arb-line and the corner point geometry to great a new volume at the arb-line locations (cdp_x, and cdp_y) that maps the sim cell index to a regular grid sample.

Convert sim2imp output to regular grids¶

In [18]:

            
                Copied!
                
from sim2x import map_to_gi_vol

report_rgs = {key: map_to_gi_vol(val, gi_vol) for key, val in imps.items()}
from sim2x import map_to_gi_vol

report_rgs = {key: map_to_gi_vol(val, gi_vol) for key, val in imps.items()}

In [19]:

            
                Copied!
                
report_rgs3d = {key: map_to_gi_vol(val, gi_vol3d) for key, val in imps.items()}
report_rgs3d = {key: map_to_gi_vol(val, gi_vol3d) for key, val in imps.items()}

In [21]:

            
                Copied!
                
                    
                    
                
                

        
# we will continue from here with just report 5
rem = report_rgs[5].copy()

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()
# we will continue from here with just report 5
rem = report_rgs[5].copy()

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()

Filling nan values¶

In [22]:

            
                Copied!
                
                    
                    
                
                

        
for key, fill_value in [
    ("vp", 3000),
    ("vs", 1550),
    ("density", 2.2),
    ("bulk_modulus", 14.5),
    ("shear_modulus", 5),
]:
    rem[key] = rem[key].where(np.invert(rem[key].isnull()), fill_value)

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()
for key, fill_value in [
    ("vp", 3000),
    ("vs", 1550),
    ("density", 2.2),
    ("bulk_modulus", 14.5),
    ("shear_modulus", 5),
]:
    rem[key] = rem[key].where(np.invert(rem[key].isnull()), fill_value)

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()

In [23]:

            
                Copied!
                
                    
                    
                
                

        
rem3d = report_rgs3d[5].copy()

for key, fill_value in [
    ("vp", 3000),
    ("vs", 1550),
    ("density", 2.2),
    ("bulk_modulus", 14.5),
    ("shear_modulus", 5),
]:
    rem3d[key] = rem3d[key].where(np.invert(rem3d[key].isnull()), fill_value)
rem3d = report_rgs3d[5].copy()

for key, fill_value in [
    ("vp", 3000),
    ("vs", 1550),
    ("density", 2.2),
    ("bulk_modulus", 14.5),
    ("shear_modulus", 5),
]:
    rem3d[key] = rem3d[key].where(np.invert(rem3d[key].isnull()), fill_value)

Smoothing the Earth Model¶

In [24]:

            
                Copied!
                
from sim2x.utils.dask import dask_gaussian_filter
from sim2x.utils.dask import dask_gaussian_filter

In [25]:

            
                Copied!
                
                    
                    
                
                

        
dask_gaussian_filter(rem, 2.0)

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()
dask_gaussian_filter(rem, 2.0)

fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
rem.vp.plot(ax=axs[0], yincrease=False)
rem.vs.plot(ax=axs[1], yincrease=False)
rem.density.plot(ax=axs[2], yincrease=False)
rem.bulk_modulus.plot(ax=axs[3], yincrease=False)
rem.shear_modulus.plot(ax=axs[4], yincrease=False)
fig.tight_layout()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:100, in _optional_import_(module, name, package)
     99 try:
--> 100     module = importlib.import_module(module)
    101     return module if name is None else importlib.import_module(name, module)

File /usr/lib/python3.8/importlib/__init__.py:127, in import_module(name, package)
    126         level += 1
--> 127 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1014, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:991, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:961, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:219, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1014, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:991, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:973, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'dask_image'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Input In [25], in <cell line: 1>()
----> 1 dask_gaussian_filter(rem, 2.0)
      3 fig, axs = plt.subplots(ncols=5, figsize=(20, 4))
      4 rem.vp.plot(ax=axs[0], yincrease=False)

File ~/.local/lib/python3.8/site-packages/sim2x/utils/dask.py:14, in dask_gaussian_filter(ds, sigma)
      7 def dask_gaussian_filter(ds: xr.Dataset, sigma: float = 0.5) -> None:
      8     """Applies filer to variables in place.
      9 
     10     Args:
     11         ds:
     12         sigma:
     13     """
---> 14     dndf = _optional_import_("dask_image.ndfilters", package="dask_image")
     15     for var in ds.data_vars:
     16         _ta = da.array(ds[var])

File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:111, in _optional_import_(module, name, package)
    108 def _failed_import(*args):
    109     raise ValueError(msg) from import_error
--> 111 _failed_import()

File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:109, in _optional_import_.<locals>._failed_import(*args)
    108 def _failed_import(*args):
--> 109     raise ValueError(msg) from import_error

ValueError: install the 'dask_image' package to make use of this feature

In [26]:

            
                Copied!
                
dask_gaussian_filter(rem3d, 1.5)
dask_gaussian_filter(rem3d, 1.5)

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:100, in _optional_import_(module, name, package)
     99 try:
--> 100     module = importlib.import_module(module)
    101     return module if name is None else importlib.import_module(name, module)

File /usr/lib/python3.8/importlib/__init__.py:127, in import_module(name, package)
    126         level += 1
--> 127 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1014, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:991, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:961, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:219, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1014, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:991, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:973, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'dask_image'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Input In [26], in <cell line: 1>()
----> 1 dask_gaussian_filter(rem3d, 1.5)

File ~/.local/lib/python3.8/site-packages/sim2x/utils/dask.py:14, in dask_gaussian_filter(ds, sigma)
      7 def dask_gaussian_filter(ds: xr.Dataset, sigma: float = 0.5) -> None:
      8     """Applies filer to variables in place.
      9 
     10     Args:
     11         ds:
     12         sigma:
     13     """
---> 14     dndf = _optional_import_("dask_image.ndfilters", package="dask_image")
     15     for var in ds.data_vars:
     16         _ta = da.array(ds[var])

File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:111, in _optional_import_(module, name, package)
    108 def _failed_import(*args):
    109     raise ValueError(msg) from import_error
--> 111 _failed_import()

File ~/.local/lib/python3.8/site-packages/sim2x/utils/tools.py:109, in _optional_import_.<locals>._failed_import(*args)
    108 def _failed_import(*args):
--> 109     raise ValueError(msg) from import_error

ValueError: install the 'dask_image' package to make use of this feature

Extracting surface from the simulation model¶

In [27]:

            
                Copied!
                
from eclx._grid import get_sim_surface

b = get_sim_surface(sim.xyzcorn, 1, face="top", slice_dir="k")
surf = arb_line.seis.surface_from_points(b[:, :2], b[:, 2])

surf.data.plot(yincrease=False)
from eclx._grid import get_sim_surface

b = get_sim_surface(sim.xyzcorn, 1, face="top", slice_dir="k")
surf = arb_line.seis.surface_from_points(b[:, :2], b[:, 2])

surf.data.plot(yincrease=False)

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Input In [27], in <cell line: 1>()
----> 1 from eclx._grid import get_sim_surface
      3 b = get_sim_surface(sim.xyzcorn, 1, face="top", slice_dir="k")
      4 surf = arb_line.seis.surface_from_points(b[:, :2], b[:, 2])

ImportError: cannot import name 'get_sim_surface' from 'eclx._grid' (/home/runner/.local/lib/python3.8/site-packages/eclx/_grid.py)

In [28]:

            
                Copied!
                
surf = vol_3d.seis.surface_from_points(b[:, :2], b[:, 2])
surf.data.plot()
surf = vol_3d.seis.surface_from_points(b[:, :2], b[:, 2])
surf.data.plot()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [28], in <cell line: 1>()
----> 1 surf = vol_3d.seis.surface_from_points(b[:, :2], b[:, 2])
      2 surf.data.plot()

NameError: name 'b' is not defined

Create the TWT/DEPTH MODEL¶

In [29]:

            
                Copied!
                
from sim2x import twt_conv

twt_int = twt_conv.time_integral(rem.vp)
twt_int.plot(yincrease=False)
# twt
from sim2x import twt_conv

twt_int = twt_conv.time_integral(rem.vp)
twt_int.plot(yincrease=False)
# twt

Out[29]:

<matplotlib.collections.QuadMesh at 0x7f7204ba83a0>

In [30]:

            
                Copied!
                
time_int3d = twt_conv.time_integral(rem3d.vp)
time_int3d.sel(iline=10).plot(yincrease=False)
time_int3d = twt_conv.time_integral(rem3d.vp)
time_int3d.sel(iline=10).plot(yincrease=False)

Out[30]:

<matplotlib.collections.QuadMesh at 0x7f7204cc63d0>

Create a 3d variation we can model

In [31]:

            
                Copied!
                
                    
                    
                
                

        
depth_samples = xr.DataArray(
    np.linspace(2430, 2490, 25), dims="cdp", coords={"cdp": twt_int.cdp}, name="depth"
)
dsx_, dsy_ = np.meshgrid(np.linspace(2430, 2490, 25), np.linspace(2430, 2490, 25))
dsm_ = (dsx_ + dsy_) / 2
depth_samples3d = xr.DataArray(
    dsm_,
    dims=("iline", "xline"),
    coords={"iline": time_int3d.iline, "xline": time_int3d.xline},
    name="depth",
)
depth_samples = xr.DataArray(
    np.linspace(2430, 2490, 25), dims="cdp", coords={"cdp": twt_int.cdp}, name="depth"
)
dsx_, dsy_ = np.meshgrid(np.linspace(2430, 2490, 25), np.linspace(2430, 2490, 25))
dsm_ = (dsx_ + dsy_) / 2
depth_samples3d = xr.DataArray(
    dsm_,
    dims=("iline", "xline"),
    coords={"iline": time_int3d.iline, "xline": time_int3d.xline},
    name="depth",
)

In [33]:

            
                Copied!
                
                    
                    
                
                

        
twt_2500mps = depth_samples * 2 / 4000
peg_twt2500 = twt_conv.peg_time_interval(
    twt_int, depth_samples, twt_2500mps, mapping_dims=("cdp",)
)

# converting to twt
rem_twt = twt_conv.domain_convert_vol(
    peg_twt2500, rem, mapping_dims=("cdp",), output_samples=np.arange(1.2, 1.27, 0.001)
)

twt3d_2500mps = depth_samples3d * 2 / 2500
peg3d_twt2500 = twt_conv.peg_time_interval(time_int3d, depth_samples3d, twt3d_2500mps)

rem3d_twt = twt_conv.domain_convert_vol(peg3d_twt2500, rem3d)
twt_2500mps = depth_samples * 2 / 4000
peg_twt2500 = twt_conv.peg_time_interval(
    twt_int, depth_samples, twt_2500mps, mapping_dims=("cdp",)
)

# converting to twt
rem_twt = twt_conv.domain_convert_vol(
    peg_twt2500, rem, mapping_dims=("cdp",), output_samples=np.arange(1.2, 1.27, 0.001)
)

twt3d_2500mps = depth_samples3d * 2 / 2500
peg3d_twt2500 = twt_conv.peg_time_interval(time_int3d, depth_samples3d, twt3d_2500mps)

rem3d_twt = twt_conv.domain_convert_vol(peg3d_twt2500, rem3d)

In [34]:

            
                Copied!
                
                    
                    
                
                

        
rem_twt_z = twt_conv.domain_convert_vol(
    peg_twt2500,
    rem_twt,
    mapping_dims=("cdp",),
    domain_vol_direction="reverse",
    from_dim="twt",
    to_dim="depth",
)
rem_twt_z = twt_conv.domain_convert_vol(
    peg_twt2500,
    rem_twt,
    mapping_dims=("cdp",),
    domain_vol_direction="reverse",
    from_dim="twt",
    to_dim="depth",
)

In [35]:

            
                Copied!
                
                    
                    
                
                

        
fig, axs = plt.subplots(ncols=4, figsize=(20, 3))
rem.vp.plot(ax=axs[0])
rem_twt.vp.plot(ax=axs[1])
rem_twt_z.vp.plot(ax=axs[2])
((rem_twt_z.vp - rem.vp) / rem.vp * 100).plot(ax=axs[3])
fig.tight_layout()
fig, axs = plt.subplots(ncols=4, figsize=(20, 3))
rem.vp.plot(ax=axs[0])
rem_twt.vp.plot(ax=axs[1])
rem_twt_z.vp.plot(ax=axs[2])
((rem_twt_z.vp - rem.vp) / rem.vp * 100).plot(ax=axs[3])
fig.tight_layout()

In [36]:

            
                Copied!
                
                    
                    
                
                

        
# check that time conversion matches in both directions
fig, axs = plt.subplots()
twt_conv.time_convert_surface(peg_twt2500, depth_samples, mapping_dims=("cdp",)).plot(
    ax=axs
)
twt_2500mps.plot(ax=axs)
# check that time conversion matches in both directions
fig, axs = plt.subplots()
twt_conv.time_convert_surface(peg_twt2500, depth_samples, mapping_dims=("cdp",)).plot(
    ax=axs
)
twt_2500mps.plot(ax=axs)

Out[36]:

[<matplotlib.lines.Line2D at 0x7f7204b28cd0>]

In [37]:

            
                Copied!
                
peg3d_z2500 = twt_conv.domain_convert_vol(
    peg3d_twt2500, xr.Dataset({"z": peg3d_twt2500.depth.broadcast_like(peg3d_twt2500)})
)
peg3d_z2500 = twt_conv.domain_convert_vol(
    peg3d_twt2500, xr.Dataset({"z": peg3d_twt2500.depth.broadcast_like(peg3d_twt2500)})
)

In [38]:

            
                Copied!
                
peg3d_z2500.sel(iline=10).z.plot()
peg3d_z2500.sel(iline=10).z.plot()

Out[38]:

<matplotlib.collections.QuadMesh at 0x7f7204a29970>

In [39]:

            
                Copied!
                
                    
                    
                
                

        
peg3d_z2500_rtwt = twt_conv.domain_convert_vol(
    peg3d_twt2500,
    peg3d_z2500,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
)
peg3d_z2500_rtwt = twt_conv.domain_convert_vol(
    peg3d_twt2500,
    peg3d_z2500,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
)

In [40]:

            
                Copied!
                
peg3d_twt2500.sel(xline=10).plot(yincrease=False)
peg3d_twt2500.sel(xline=10).plot(yincrease=False)

Out[40]:

<matplotlib.collections.QuadMesh at 0x7f71fd6a8760>

calculate reflectivity volume¶

In [41]:

            
                Copied!
                
from sim2x import seis_synth, seis_interface
from sim2x import seis_synth, seis_interface

In [42]:

            
                Copied!
                
                    
                    
                
                

        
a = seis_synth.reflectivity(
    5,
    rem_twt.sel(cdp=10).vp.values,
    rem_twt.sel(cdp=10).vs.values,
    rem_twt.sel(cdp=10).density.values,
    method="zoep_pud",
)
a = seis_synth.reflectivity(
    5,
    rem_twt.sel(cdp=10).vp.values,
    rem_twt.sel(cdp=10).vs.values,
    rem_twt.sel(cdp=10).density.values,
    method="zoep_pud",
)

In [43]:

            
                Copied!
                
plt.plot(a["Rp"])
plt.plot(a["Rp"])

Out[43]:

[<matplotlib.lines.Line2D at 0x7f71fcd1f310>]

In [44]:

            
                Copied!
                
rf = seis_synth.reflectivity_vol(rem_twt, 5, mapping_dims=("cdp",))

rf = seis_synth.reflectivity_vol(rem_twt, 5, mapping_dims=("cdp",))

In [45]:

            
                Copied!
                
rf.refl_5.T.plot(yincrease=False)
rf.refl_5.T.plot(yincrease=False)

Out[45]:

<matplotlib.collections.QuadMesh at 0x7f71fd2f31c0>

In [46]:

            
                Copied!
                
rf.refl_5.T.plot(yincrease=False)
rf.refl_5.T.plot(yincrease=False)

Out[46]:

<matplotlib.collections.QuadMesh at 0x7f71fcd82610>

In [47]:

            
                Copied!
                
from sim2x import wavelets as wv

ricker = wv.RickerWavelet(
    25,
)
ricker.resample(dt=0.001)
ricker.as_seconds()
from sim2x import wavelets as wv

ricker = wv.RickerWavelet(
    25,
)
ricker.resample(dt=0.001)
ricker.as_seconds()

In [48]:

            
                Copied!
                
seis = seis_synth.convolution1d_vol(rf, "refl_5", ricker.amp, mapping_dims=("cdp",))
seis = seis_synth.convolution1d_vol(rf, "refl_5", ricker.amp, mapping_dims=("cdp",))

In [49]:

            
                Copied!
                
seis["amp"].T.plot(yincrease=False)
seis["amp"].T.plot(yincrease=False)

Out[49]:

<matplotlib.collections.QuadMesh at 0x7f71fd903730>

In [50]:

            
                Copied!
                
                    
                    
                
                

        
seis_depth = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="linear",
)
seis_depth_cubic = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="cubic",
)
seis_depth_sinc = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="sinc",
)

fig, axs = plt.subplots(ncols=3, nrows=3, figsize=(15, 8))

seis_depth["amp"].plot(yincrease=False, ax=axs[0, 0])
seis_depth_cubic["amp"].plot(yincrease=False, ax=axs[1, 1])
seis_depth_sinc["amp"].plot(yincrease=False, ax=axs[2, 2])

(seis_depth["amp"] - seis_depth_cubic["amp"]).plot(yincrease=False, ax=axs[0, 1])
(seis_depth["amp"] - seis_depth_sinc["amp"]).plot(yincrease=False, ax=axs[0, 2])

(seis_depth_cubic["amp"] - seis_depth_sinc["amp"]).plot(yincrease=False, ax=axs[1, 2])

cdp = 12
seis_depth["amp"].sel(cdp=cdp).plot(ax=axs[1, 0], label="linear")
seis_depth_cubic["amp"].sel(cdp=cdp).plot(ax=axs[1, 0], label="cubic")

seis_depth["amp"].sel(cdp=cdp).plot(ax=axs[2, 0], label="linear")
seis_depth_sinc["amp"].sel(cdp=cdp).plot(ax=axs[2, 0], label="sinc")

seis_depth_cubic["amp"].sel(cdp=cdp).plot(ax=axs[2, 1], label="cubic")
seis_depth_sinc["amp"].sel(cdp=cdp).plot(ax=axs[2, 1], label="sinc")

fig.tight_layout()
seis_depth = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="linear",
)
seis_depth_cubic = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="cubic",
)
seis_depth_sinc = twt_conv.domain_convert_vol(
    peg_twt2500,
    seis,
    from_dim="twt",
    to_dim="depth",
    domain_vol_direction="reverse",
    mapping_dims=("cdp",),
    interpolation_kind="sinc",
)

fig, axs = plt.subplots(ncols=3, nrows=3, figsize=(15, 8))

seis_depth["amp"].plot(yincrease=False, ax=axs[0, 0])
seis_depth_cubic["amp"].plot(yincrease=False, ax=axs[1, 1])
seis_depth_sinc["amp"].plot(yincrease=False, ax=axs[2, 2])

(seis_depth["amp"] - seis_depth_cubic["amp"]).plot(yincrease=False, ax=axs[0, 1])
(seis_depth["amp"] - seis_depth_sinc["amp"]).plot(yincrease=False, ax=axs[0, 2])

(seis_depth_cubic["amp"] - seis_depth_sinc["amp"]).plot(yincrease=False, ax=axs[1, 2])

cdp = 12
seis_depth["amp"].sel(cdp=cdp).plot(ax=axs[1, 0], label="linear")
seis_depth_cubic["amp"].sel(cdp=cdp).plot(ax=axs[1, 0], label="cubic")

seis_depth["amp"].sel(cdp=cdp).plot(ax=axs[2, 0], label="linear")
seis_depth_sinc["amp"].sel(cdp=cdp).plot(ax=axs[2, 0], label="sinc")

seis_depth_cubic["amp"].sel(cdp=cdp).plot(ax=axs[2, 1], label="cubic")
seis_depth_sinc["amp"].sel(cdp=cdp).plot(ax=axs[2, 1], label="sinc")

fig.tight_layout()

In [ ]: