Module curlew.synthetic

Generate synthetic datasets and models for testing purposes.

Functions

def anderson(shape=None, offset1=225, offset2=250, **kwargs)

Expand source code

def anderson( shape=None, offset1=225, offset2=250, **kwargs ):

    """
    Return a synthetic model with a slightly curved layer-cake stratigraphy cut
    by two intersecting normal faults.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    offset1 : tuple
        The offset of the first (older) normal fault.
    offet2 : tuple
        The offset of the second (younger) normal fault.
    
    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`
    
    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)
    vert = 1 if ndim == 2 else 2

    s0 = strati('s0', C=QuadraticField( 'f0', input_dim=ndim, curve=_g2((-0.00005, 0), ndim), origin=_p2((1000, 500), ndim)  ) )
    s1 = fault( 's1',
           C=LinearField( 'f1', input_dim=ndim, origin=_p2((950, 550), ndim), gradient=_g2((np.cos( np.deg2rad(35) ), np.sin( np.deg2rad(35) )), ndim)  ),
           offset=offset1, shortening=_g2((0, -1), ndim) )
    s2 = fault( 's2',
           C=LinearField( 'f2', input_dim=ndim, origin=_p2((1050, 500), ndim), gradient=_g2((-np.cos( np.deg2rad(35) ), np.sin( np.deg2rad(35) )), ndim)  ),
           offset=offset2, shortening=_g2((0, -1), ndim) )
    
    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
        
    M = GeoModel( [s0, s1, s2], grid=G, name='anderson' )
    s = M.predict(G)

    Cs = sample(s, M, pv='rgb', bind=False, **kwargs)
    kw_fault = dict(kwargs)
    kw_fault.pop('breaks', None)
    kw_fault.pop('pv', None)
    kw_fault['pval'] = kw_fault.get('pval', 1.0)
    Cf1 = sample(s1.predict(G), M, pv='rgb', breaks=[0.5], xstep=600, bind=False, **kw_fault)
    Cf2 = sample(s2.predict(G), M, pv='rgb', breaks=[0.5], bind=False, **kw_fault)

    ordered = _csets_in_order(Cs)
    strat = ordered[0]
    prop = ordered[-1] if 'property' in Cs else None
    f1_c = _csets_in_order(Cf1)[0]
    f2_c = _csets_in_order(Cf2)[0]

    gv = strat.gv
    if gv is not None:
        mask = (gv[:, vert] > 0.9) & (gv[:, vert] < 1.1)
        strat.gp = strat.gp[mask]
        strat.gv = strat.gv[mask]
        strat.gop = strat.gop[mask]
        strat.gov = strat.gov[mask]

    for fc in (f1_c, f2_c):
        if fc.vv is not None:
            fc.vv = fc.vv * 0

    M['s0'].field.bind(strat)
    M['s1'].field.bind(f1_c)
    M['s2'].field.bind(f2_c)
    if prop is not None:
        M.bind(prop)

    return M

Return a synthetic model with a slightly curved layer-cake stratigraphy cut by two intersecting normal faults.

Parameters

shape : tuple

The width and height of the generated data.

offset1 : tuple

The offset of the first (older) normal fault.

offet2 : tuple

The offset of the second (younger) normal fault.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def extract_constraints(M, events=None, include_property=False)

Expand source code

def extract_constraints(M, events=None, include_property=False):
    """
    Return bound constraint sets from a synthetic ``curlew.geology.geomodel.GeoModel`` for use when fitting neural fields.

    Parameters
    ----------
    M : curlew.geology.geomodel.GeoModel
        A synthetic model returned by one of the builders in this module (after ``sample`` has run).
    events : list of str, optional
        Event names to include. Defaults to all events in ``M.events``.
    include_property : bool
        If True, include global property constraints under the key ``'property'``.

    Returns
    -------
    dict
        Constraint sets keyed by event name (and ``'property'`` when requested).
    """
    if events is None:
        events = [e.name for e in M.events]
    out = {}
    for name in events:
        c = M[name].field.C
        if c is not None:
            out[name] = c.numpy()
    if include_property and M.C is not None:
        out["property"] = M.C.numpy()
    return out

Return bound constraint sets from a synthetic GeoModel for use when fitting neural fields.

Parameters

M : GeoModel

A synthetic model returned by one of the builders in this module (after sample() has run).

events : list of str, optional

Event names to include. Defaults to all events in M.events.

include_property : bool

If True, include global property constraints under the key 'property'.

Returns

dict

Constraint sets keyed by event name (and 'property' when requested).

def goguel(shape=None, zstep=300, dip=60.0, offset=0.2, n_faults=3, **kwargs)

Expand source code

def goguel(shape=None, zstep=300, dip=60.0, offset=0.2, n_faults=3, **kwargs):
    """
    Return a synthetic model with horizontal stratigraphy and several finite faults.

    Each fault has planar geometry (shared implicit field), constant slip scaled by a
    cosine-tapered ellipsoidal patch (finite fault extent).

    Parameters
    ----------
    shape : tuple, optional
        The width and height of the generated data. 
    zstep : float
        Vertical spacing between fault seed positions.
    dip : float
        Fault dip in degrees (from horizontal).
    offset : float
        Slip magnitude. Values in ``(0, 1)`` are interpreted as a fraction of the
        ellipsoid major-axis radius; values ``>= 1`` are absolute slip (same units as extent).
    n_faults : int
        Number of finite faults (default 3).

    Keywords
    --------
        Passed to :func:`curlew.synthetic.sample`.

    Returns
    -------
    M : curlew.geology.geomodel.GeoModel.
    """
    # Grid, extent, and model centre (fault/strat seeds are placed relative to these).
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)
    vert = 1 if ndim == 2 else 2
    origin = np.asarray(G.center, dtype=float)

    # Ellipsoid major axis and fault spacing scale with model height; slip from offset arg.
    patch_radius = 2.5 * extent[vert] / (n_faults + 1)
    xstep = 1.5 * extent[0] / (n_faults + 1)
    if isinstance(offset, (int, np.integer)) or (isinstance(offset, float) and offset >= 1):
        slip_mag = float(offset)
    else:
        slip_mag = offset * patch_radius
    slip = [slip_mag] * n_faults

    # Stagger fault centres in x and z, then re-centre the cluster on the grid midpoint.
    cx = np.array([xstep * (i + 1) for i in range(n_faults)])
    cz = np.array([extent[vert] - zstep * (i + 1) for i in range(n_faults)])
    cx = cx - np.mean(cx) + origin[0]
    cz = cz - np.mean(cz) + origin[vert]

    # Layer-cake basement: horizontal scalar field and regularly spaced isosurface seeds.
    s0 = strati(
        "s0",
        type=LinearField,
        origin=_g2(origin, ndim),
        gradient=_g2(np.array([0.0, 1.0]), ndim),
    )
    for seed, name in _isosurface_seeds(extent, ndim, 8):
        s0.addIsosurface(name, seed=seed)

    # One fault GeoEvent: shared plane dipping toward +x; vertical shortening gives normal (top-down) slip.
    fault_grad = np.array([np.sin(np.deg2rad(dip)), np.cos(np.deg2rad(dip))])
    contacts = [f"ff{i}" for i in range(n_faults)]
    modifiers = [f"fe{i}" for i in range(n_faults)]
    s1 = fault(
        "s1",
        type=LinearField,
        origin=_g2(origin, ndim),
        gradient=_g2(fault_grad, ndim),
        shortening=_g2((0, -1), ndim),
        offset=slip,
        contact=contacts,
        width=0,
        modifier=modifiers,
        n_steps=1,
    )

    # Ellipse frame: tangent along strike, normal across the fault plane; axes set patch size.
    fn = _g2(s1["s1"].grad.detach().cpu().numpy(), ndim)
    ft = _g2(np.array([-fn[1], fn[0]]), ndim)
    axes = _g2(np.array([patch_radius, 0.75 * patch_radius]), ndim)

    # Per fault: isosurface seed on the shared plane + cosine ellipsoid that tapers slip to zero.
    for i, (_x, _z) in enumerate(zip(cx, cz)):
        seed = _g2(np.array([_x, _z]), ndim)
        s1.addIsosurface(
            f"ff{i}",
            seed=seed.copy(),
            field="s1",
        )
        if ndim == 2:
            directions = np.array([ft, fn])
        else:
            directions = np.stack([ft, np.cross(ft, fn), fn])
        s1.addField(
            f"fe{i}",
            field=EllipsoidalField(
                f"fe{i}",
                input_dim=ndim,
                origin=seed.copy(),
                axes=axes,
                directions=directions,
                decay="cosine",
            ),
        )

    # Assemble model and sample constraints (bedding on s0, fault plane on s1, global property).
    M = GeoModel([s0, s1], grid=G, name="goguel")
    s = M.predict(G)

    Cs = sample(s, M, pv="rgb", bind=False, **kwargs)
    kw_fault = dict(kwargs)
    kw_fault.pop("breaks", None)
    kw_fault["pval"] = kw_fault.get("pval", 1.0)
    fault_breaks = [s1.getIsovalue(f"ff{i}") for i in range(n_faults)]
    Cf = sample(
        s1.predict(G),
        M,
        pv="rgb",
        breaks=fault_breaks,
        bind=False,
        **kw_fault,
    )

    ordered = _csets_in_order(Cs)
    strat = ordered[0]
    prop = ordered[-1] if "property" in Cs else None
    fault_c = _csets_in_order(Cf)[0]

    # Keep near-horizontal bedding orientations only (s0 gradient is vertical).
    gv = strat.gv
    if gv is not None:
        mask = (gv[:, vert] > 0.9) & (gv[:, vert] < 1.1)
        strat.gp = strat.gp[mask]
        strat.gv = strat.gv[mask]
        strat.gop = strat.gop[mask]
        strat.gov = strat.gov[mask]

    if fault_c.vv is not None:
        fault_c.vv = fault_c.vv * 0

    M["s0"].field.bind(strat)
    M["s1"]["s1"].bind(fault_c)
    if prop is not None:
        M.bind(prop)

    return M

Return a synthetic model with horizontal stratigraphy and several finite faults.

Each fault has planar geometry (shared implicit field), constant slip scaled by a cosine-tapered ellipsoidal patch (finite fault extent).

Parameters

shape : tuple, optional

The width and height of the generated data.

zstep : float

Vertical spacing between fault seed positions.

dip : float

Fault dip in degrees (from horizontal).

offset : float

Slip magnitude. Values in (0, 1) are interpreted as a fraction of the ellipsoid major-axis radius; values >= 1 are absolute slip (same units as extent).

n_faults : int

Number of finite faults (default 3).

Keywords

Passed to :func:<code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a></code>.

Returns

M : curlew.geology.geomodel.GeoModel.

def hutton(shape=None, **kwargs)

Expand source code

def hutton( shape=None, **kwargs ):
    """
    Return a synthetic model with a folded basement cut by an unconformity.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    
    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`

    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)

    s0 = strati('s0', C=QuadraticField( 'f1', input_dim=ndim, gradient=_g2((0, 1), ndim), 
                                       origin=_p2((0, 0), ndim), curve=_g2((-0.00002, 0), ndim) ))
    d1 = fold('d1', origin=_p2((0, 0), ndim),
                    extension=_g2((0, 1), ndim),
                    compression=_g2((1, 0.3), ndim),
                    wavelength=2000,
                    amplitude=250, sharpness=0.7)
    s1 = strati('s1', C=QuadraticField( 'f1', input_dim=ndim, gradient=_g2((0.1, 0.9), ndim), 
                                       origin=_p2((1000, 500), ndim), curve=_g2((-0.00002, 0), ndim) ), base=0)

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
    for seed, name in _isosurface_seeds(extent, ndim, 10):
        s1.addIsosurface(name, seed=seed)

    M  = GeoModel( [s0,d1,s1], grid=G, name='hutton' )
    s = M.predict(G)
    Cs = sample(s, M, pv='rgb', bind=False, **kwargs)
    ordered = _csets_in_order(Cs)
    M['s0'].field.bind(ordered[1])
    M['s1'].field.bind(ordered[0])
    if len(ordered) > 2:
        M.bind(ordered[2])

    return M

Return a synthetic model with a folded basement cut by an unconformity.

Parameters

shape : tuple

The width and height of the generated data.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def lehmann(shape=None, **kwargs)

Expand source code

def lehmann( shape=None, **kwargs ):
    """
    Return a synthetic model with a folded basement.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    
    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`

    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)

    if ndim == 3:
        s0 = strati('s0', C=PeriodicField('f0', input_dim=ndim, gradient=np.array([0., 0., 1.]), axialPlane=np.array([1., 0., 0.])))
    else:
        s0 = strati('s0', C=PeriodicField('f0', input_dim=ndim))
    
    M = GeoModel([s0], grid=G, name='lehmann' )

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
    
    s = M.predict(G)
    sample(s, M, pv='rgb', **kwargs)
    
    return M

Return a synthetic model with a folded basement.

Parameters

shape : tuple

The width and height of the generated data.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def michell(shape=None, offset=100, **kwargs)

Expand source code

def michell( shape=None, offset=100, **kwargs ):

    """
    Return a synthetic model with a slightly curved layer-cake stratigraphy cut
    by a thrust fault.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    offset : tuple
        The offset of the fault.

    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`
    
    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)
    o = _p2((1000, 500), ndim)

    s0 = strati('s0', C=QuadraticField( 'f0', input_dim=ndim, curve=_g2((-0.00005, 0), ndim), origin=o ) )
    s1 = fault( 's1',
           C=LinearField( 'f1', input_dim=ndim, origin=o, gradient=_g2((0.5, 0.5), ndim)  ),
           offset=offset, shortening=_g2((-1, 0), ndim) )

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)

    M = GeoModel( [s0,s1], grid=G, name='michell' )
    s = M.predict(G)
    Cs = sample(s, M, pv='rgb', bind=False, **kwargs)
    kwargs['pval'] = kwargs.get('pval', 1.0)
    kw_fault = dict(kwargs)
    kw_fault.pop('breaks', None)
    Cf = sample(s1.predict(G), M, pv='rgb', breaks=[0.5], bind=False, **kw_fault)

    ordered = _csets_in_order(Cs)
    strat = ordered[0]
    fault_c = _csets_in_order(Cf)[0]
    prop = ordered[-1] if 'property' in Cs else None

    gv = strat.gv
    if gv is not None:
        mask = gv[:, 0] > -0.68
        strat.gp = strat.gp[mask]
        strat.gv = strat.gv[mask]
        strat.gop = strat.gop[mask]
        strat.gov = strat.gov[mask]

    if fault_c.vv is not None:
        fault_c.vv = fault_c.vv * 0

    M['s0'].field.bind(strat)
    M['s1'].field.bind(fault_c)
    if prop is not None:
        M.bind(prop)

    return M

Return a synthetic model with a slightly curved layer-cake stratigraphy cut by a thrust fault.

Parameters

shape : tuple

The width and height of the generated data.

offset : tuple

The offset of the fault.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def playfair(shape=None, width=50, addFault=False, **kwargs)

Expand source code

def playfair( shape=None, width=50, addFault=False, **kwargs ):

    """
    Return a synthetic model with a layer-cake stratigraphy cut
    by a dyke.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    width : float
        The half-width of the added dyke.
    addFault : bool
        If True, add a fault crosscutting the dykes and stratigraphy. Default is False. 

    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`
    
    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)
    o = _p2((1000, 500), ndim)

    s0 = strati('s0', C=QuadraticField( 'f0', input_dim=ndim, curve=_g2((-0.00005, 0), ndim), origin=o ) )
    s1 = sheet( 's1',
           C=LinearField( 'f1', input_dim=ndim, origin=o, gradient=_g2((0.5, 0.5), ndim) ), contact=(-width, width) )

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
    
    if addFault:
        s2 = fault( 's2',
           C=LinearField( 'f2', input_dim=ndim, origin=_p2((1050, 500), ndim), gradient=_g2((-np.cos( np.deg2rad(35) ), np.sin( np.deg2rad(35) )), ndim)  ),
           offset=100, shortening=_g2((0, -1), ndim) )
        M = GeoModel( [s0, s1, s2], grid=G, name='newcastle' )
    else:
        M = GeoModel( [s0, s1], grid=G, name='playfair' )
    
    kwargs['pval'] = kwargs.get('pval', 1.0) # change default to sample all value constraints
    Cs = sample(M.predict(G), M, pv='rgb', bind=False, **kwargs)
    C1 = sample(s1.predict(G), M, pv='rgb', breaks=[-width, width], bind=False, **kwargs)
    ordered = _csets_in_order(Cs)
    dyke = _csets_in_order(C1)[0]
    M['s0'].field.bind(ordered[1])
    M['s1'].field.bind(dyke)
    if len(ordered) > 2:
        M.bind(ordered[2])

    return M

Return a synthetic model with a layer-cake stratigraphy cut by a dyke.

Parameters

shape : tuple

The width and height of the generated data.

width : float

The half-width of the added dyke.

addFault : bool

If True, add a fault crosscutting the dykes and stratigraphy. Default is False.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def sample(G, M, pv=None, breaks=19, init=100, xstep=300, pval=0.6, cmap='tab20', seed=42, bind=True)

Expand source code

def sample( G, M, pv=None, breaks=19, init=100, xstep=300, pval=0.6, cmap='tab20', seed=42, bind=True ):
    """
    Sample value, orientation, equality (contact trace), and property constraints from a scalar field and associated gradients.

    Parameters
    ----------
    G : curlew.core.Geode
        A Geode object containing the results from a GeoEvent or GeoModel.
    M : curlew.geology.geomodel.GeoModel
        The geomodel used to construct the geode G. 
    pv : np.ndarray or str
        A (N, d) array of n-dimensional property vectors (e.g., color). Can also be 'rgb' to create synthetic colors.
    breaks : list
        A set of scalar values at which contacts (changes in color) should be placed.
    init : int
        The index of the first "drillhole" in the x-direction.
    xstep : int
        The separation between "drillholes" in the x-direction.
    pval : float
        The probability that an observation is a scalar value (rather than just a gradient) constraint.
    cmap : str
        The name of the Matplotlib colormap to use for sampling colors that determine where geological contacts are.
        Must be a discrete colormap.
    seed : int
        Random seed to facilitate reproducible results.
    bind : bool
        If True, the sampled constraints will be added to their respective fields (and global constraints bound to the
        GeoModel).
        
    Notes
    -----
    For 3D grids, constraints are sampled on x-z sections at y = 25%, 50% and 75%
    of the model y extent (1500 m).

    Returns
    -------
    dict
        Sampled constraints keyed by structure ID (int) and ``'property'`` when present.
        Empty when ``bind=True`` (constraints are attached to ``M`` and its events).
    """

    np.random.seed(seed)

    sf = G.grid.reshape(curlew._numpy(G.scalar))
    sid = G.grid.reshape(curlew._numpy(G.structureID).astype(int))
    xy = G.grid.reshape(G.grid.coords())
    constraints = {
        int(k): {"vp": [], "vv": [], "gp": [], "gv": [], "gop": [], "gov": [], "eq": []}
        for k in np.unique(sid)
    }
    eq_accum = {int(k): {} for k in np.unique(sid)}

    if G.grid.ndim == 3:
        y_global = xy[0, :, 0, 1]
        for y_frac in (0.25, 0.5, 0.75):
            j = int(np.argmin(np.abs(y_global - G.grid.dims[1] * y_frac)))
            sf2 = sf[:, j, :]
            sid2 = sid[:, j, :]
            xy2 = xy[:, j, :, :]
            gx = np.diff(sf2, axis=0)
            gz = np.diff(sf2, axis=1)
            c = colour(sf2, breaks=breaks, cmap=cmap)
            breaks_arr = _breaks_array(sf2, breaks)
            contacts = np.sum(np.abs(np.diff(c, axis=1, append=0)), axis=-1) > 0
            domains = np.sum(
                [np.abs(np.diff(sid2, axis=1, append=0)), np.abs(np.diff(sid2, axis=0, append=0))], axis=0
            ) > 0
            pv2 = c if pv == "rgb" else pv
            _sample_section(
                sf2, sid2, xy2, contacts, domains, gx, gz, init, xstep, pval, pv2,
                constraints, eq_accum, breaks_arr,
            )
    else:
        gx = np.diff(sf, axis=0)
        gy = np.diff(sf, axis=1)
        c = colour(sf, breaks=breaks, cmap=cmap)
        breaks_arr = _breaks_array(sf, breaks)
        contacts = np.sum(np.abs(np.diff(c, axis=1, append=0)), axis=-1) > 0
        domains = np.sum(
            [np.abs(np.diff(sid, axis=1, append=0)), np.abs(np.diff(sid, axis=0, append=0))], axis=0
        ) > 0
        if pv == "rgb":
            pv = c
        _sample_section(
            sf, sid, xy, contacts, domains, gx, gy, init, xstep, pval, pv,
            constraints, eq_accum, breaks_arr,
        )

    for i, levels in eq_accum.items():
        traces = [np.array(pts) for pts in levels.values() if len(pts) >= 2]
        if traces:
            constraints[i]["eq"].extend(traces)

    out = _finalize_constraints(constraints, M, bind=bind)
    return {} if bind else out

Sample value, orientation, equality (contact trace), and property constraints from a scalar field and associated gradients.

Parameters

G : Geode

A Geode object containing the results from a GeoEvent or GeoModel.

M : GeoModel

The geomodel used to construct the geode G.

pv : np.ndarray or str

A (N, d) array of n-dimensional property vectors (e.g., color). Can also be 'rgb' to create synthetic colors.

breaks : list

A set of scalar values at which contacts (changes in color) should be placed.

init : int

The index of the first "drillhole" in the x-direction.

xstep : int

The separation between "drillholes" in the x-direction.

pval : float

The probability that an observation is a scalar value (rather than just a gradient) constraint.

cmap : str

The name of the Matplotlib colormap to use for sampling colors that determine where geological contacts are. Must be a discrete colormap.

seed : int

Random seed to facilitate reproducible results.

bind : bool

If True, the sampled constraints will be added to their respective fields (and global constraints bound to the GeoModel).

Notes

For 3D grids, constraints are sampled on x-z sections at y = 25%, 50% and 75% of the model y extent (1500 m).

Returns

dict

Sampled constraints keyed by structure ID (int) and 'property' when present. Empty when bind=True (constraints are attached to M and its events).

def seuss(shape=None, nlayers=6, **kwargs)

Expand source code

def seuss(shape=None, nlayers=6, **kwargs):
    """
    Return a synthetic model with layered stratigraphy, a dyke, an intrusion
    domain boundary, and two listric faults (one older, one younger), matching
    the "Seuss" GeoModel from the Building Analytical Models tutorial.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. Default (1500, 700) matches
        the tutorial notebook.
    nlayers : int
        Number of stratigraphic layers (isosurfaces) in each package. Default 6.

    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)` when
        generating constraints.

    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel with listric faults and domain boundaries (name='Seuss').
    """
    shape = shape or EXTENT_2D
    G = _make_grid(shape)
    dims = _physical_extent(shape)
    
    # First stratigraphic package (layer-cake)
    s0 = strati(
        "s0",
        C=LinearField("f0", input_dim=2, gradient=(0.1, 0.9)),
    )
    sy1 = np.linspace(0, dims[1], nlayers)
    sx1 = np.full_like(sy1, dims[0] / 10)
    for i, (x, y) in enumerate(zip(sx1, sy1)):
        s0.addIsosurface("S%d" % (i + 1), seed=np.array([x, y]))

    # Dyke
    s1 = sheet(
        "s1",
        C=LinearField(
            "f1",
            input_dim=2,
            origin=np.array([750.0, 300.0]),
            gradient=np.array([-0.5, 0.5]),
        ),
        contact=(-20, 20),
        aperture=2,
    )

    # Domain boundary: intrusion (constant -2) below, s0 and s1 above
    intrusion = GeoEvent("i0", field=-2, type=int) # quick and easy way to define a constant scalar field
    d1 = domainBoundary(
        "d1",
        C=QuadraticField(
            "d1f",
            input_dim=2,
            origin=np.array([400.0, 50.0]),
            gradient=np.array([0.4, 1.0]),
            curve=(0, 0.002),
        ),
        bound=0,
        lt=[intrusion],
        gt=[s0, s1],
    )

    # Older listric fault
    f1 = fault(
        "f1",
        C=ListricField(
            "f1f",
            input_dim=2,
            origin=np.array([650.0, 700.0]),
            fault_floor=0.0,
            fault_ceil=700.0,
            curvature_rate=0.006,
        ),
        width=1e-9,
        n_steps=3,
        offset=100,
        shortening=np.array([1, -1]),
    )

    # Second stratigraphic package (unconformable cover)
    s2 = strati(
        "s2",
        C=LinearField(
            "f2",
            input_dim=2,
            origin=np.array([0.0, 500.0]),
            gradient=np.array([0.0, 1.0]),
        ),
    )
    sy2 = np.linspace(350, dims[1], nlayers)
    sx2 = np.full_like(sy2, 9 * dims[0] / 10)
    for i, (x, y) in enumerate(zip(sx2, sy2)):
        s2.addIsosurface("S%d" % (i + 1), seed=np.array([x, y]))

    # Unconformity domain boundary: d1 and f1 below, s2 above
    d2 = domainBoundary(
        "d2",
        C=LinearField(
            "d2f",
            input_dim=2,
            origin=np.array([700.0, 500.0]),
            gradient=np.array([0.2, 0.9]),
        ),
        bound=0,
        lt=[d1, f1],
        gt=[s2],
    )
    
    f2 = fault('f2', type=ListricField,
           contact=0.0, 
           offset=100, # The displacement magnitude
           n_steps=3, # apply displacement in several steps, due to high curvature
           shortening=[1, -1],
           input_dim=2,
           origin=np.array([250.0, 700.0]),
           fault_ceil=700., # Defines the ceiling of the fault
           fault_floor=0., # Defines the floor of the fault
           width=1e-9,
           curvature_rate=0.006 # Defines the steepness of the decay
           
        )
    
    M = GeoModel([d2, f2], grid=G, name="Seuss")
    s = M.predict(G)
    sample(s, M, pv="rgb", **kwargs)
    return M

Return a synthetic model with layered stratigraphy, a dyke, an intrusion domain boundary, and two listric faults (one older, one younger), matching the "Seuss" GeoModel from the Building Analytical Models tutorial.

Parameters

shape : tuple

The width and height of the generated data. Default (1500, 700) matches the tutorial notebook.

nlayers : int

Number of stratigraphic layers (isosurfaces) in each package. Default 6.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code> when
generating constraints.

Returns

M : GeoModel

Geomodel with listric faults and domain boundaries (name='Seuss').

def steno(shape=None, **kwargs)

Expand source code

def steno( shape=None, **kwargs ):
    """
    Return a synthetic model with a slightly curved layer-cake stratigraphy.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    
    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`
    
    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)

    s0 = strati('s0', C=QuadraticField( 'f0', 
                                       input_dim=ndim, 
                                       gradient=_g2((0.00001, 1), ndim), 
                                       curve=_g2((-0.00005, 0), ndim), 
                                       origin=_p2((1000, 500), ndim) ) )

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
    
    M = GeoModel([s0], grid=G, name="steno")

    s = s0.predict(G)
    sample(s, M, pv='rgb', **kwargs)

    return M

Return a synthetic model with a slightly curved layer-cake stratigraphy.

Parameters

shape : tuple

The width and height of the generated data.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.

def walker(shape=None,
width=[60, 50, 40, 50, 40],
pos=[0, 100, 200, 400, 600],
addFault=True,
**kwargs)

Expand source code

def walker( shape=None, width=[60,50,40,50,40], pos=[0,100,200,400,600], addFault=True, **kwargs ):

    """
    Return a synthetic model with a layer-cake stratigraphy cut
    by several parallel dykes.

    Parameters
    ----------
    shape : tuple
        The width and height of the generated data. 
    width : float
        A list specifying the dyke widths.
    pos : float
        A list specifying the dyke positions in the x-axis. Must have the same length as width.
    addFault : bool
        If True, add a fault crosscutting the dykes and stratigraphy. Default is False.
    Keywords
    ---------
        All keywords are passed to `curlew.synthetic.sample(...)`
    
    Returns
    --------
    M : curlew.geology.geomodel.GeoModel
        Geomodel of the synthetic model with sampled constraints bound to events.
    """
    G = _make_grid(shape)
    extent = _physical_extent(shape)
    ndim = _model_ndim(shape)
    o = _p2((0, extent[1] / 2 if ndim == 2 else extent[2] / 2), ndim)
    s0 = strati('s0', C=QuadraticField( 'f0', input_dim=ndim, curve=_g2((-0.00005, 0), ndim), origin=o ) )
    contact = [(pos[i], pos[i]+width[i]) for i in range(len(pos))]
    s1 = sheet( 's1',
           C=LinearField( 'f1', input_dim=ndim, origin=o, gradient=_g2((0.5, 0.5), ndim) ), contact=contact )

    for seed, name in _isosurface_seeds(extent, ndim, 5):
        s0.addIsosurface(name, seed=seed)
    
    if addFault:
        s2 = fault( 's2',
           C=LinearField( 'f2', input_dim=ndim, origin=_p2((1050, 500), ndim), gradient=_g2((-np.cos( np.deg2rad(35) ), np.sin( np.deg2rad(35) )), ndim)  ),
           offset=100, shortening=_g2((0, -1), ndim) )
        M = GeoModel( [s0, s1, s2], grid=G, name='newcastle' )
    else:
        M = GeoModel( [s0, s1], grid=G, name='walker' )
    
    kwargs['pval'] = kwargs.get('pval', 1.0) # change default to sample all value constraints
    Cs = sample(M.predict(G), M, pv='rgb', bind=False, **kwargs)
    C1 = sample(s1.predict(G), M, pv='rgb', breaks=np.hstack(contact), bind=False, **kwargs)
    ordered = _csets_in_order(Cs)
    dyke = _csets_in_order(C1)[0]
    M['s0'].field.bind(ordered[1])
    M['s1'].field.bind(dyke)
    if len(ordered) > 2:
        M.bind(ordered[2])

    return M

Return a synthetic model with a layer-cake stratigraphy cut by several parallel dykes.

Parameters

shape : tuple

The width and height of the generated data.

width : float

A list specifying the dyke widths.

pos : float

A list specifying the dyke positions in the x-axis. Must have the same length as width.

addFault : bool

If True, add a fault crosscutting the dykes and stratigraphy. Default is False.

Keywords

All keywords are passed to <code><a title="curlew.synthetic.sample" href="#curlew.synthetic.sample">sample()</a>(...)</code>

Returns

M : GeoModel

Geomodel of the synthetic model with sampled constraints bound to events.