curlew documentation|Stable (main)Development (dev)

Module curlew.geology

Functions and objects used for creating and manipulating geological structures in Curlew.

Sub-modules

curlew.geology.geoevent

A class defining a structural scalar (implicit) field. This is where a lot of the magic happens, including the chaining of multiple neural fields and …

curlew.geology.geomodel

A class representing a time-aware geological model and facilitating interactions with the underlying linked-list of GeoEvent instances (that represent …

curlew.geology.interactions

Functions defining how different scalar fields interact with each other to create generative, kinematic and hybrid events …

curlew.geology.property

Implementation of various types of "geological" forward models that convert scalar field values to measurable rock properties (density, magnetic …

Functions

def domainBoundary(name,
*,
C,
bound=0,
gt=0,
lt=1,
mode='below',
lithoSharpness=1.0,
structureSharpness=1.0,
**kwargs)
Expand source code
def domainBoundary( name, *, C, bound = 0, gt = 0, lt = 1, mode="below",
                    lithoSharpness=1.0, structureSharpness=1.0, **kwargs ):
    """
    Create a GeoEvent representing a domain boundary, in which different sub-models are modelled on either side of the boundary. 
    This can be very useful for modelling e.g., sedimentary basins, where the basin fill is modelled on one side of the boundary (onlap relations), 
    or for modelling domain boundary faults (in which there is no known or meaningful relationship between the fault's hangingwall and footwall).

    Parameters
    ------------
    name : str
        A name for the created domain boundary.
    C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
        Either a pre-constructed neural field, explicit field or a set of 
        constraints to use to construct a new GeoEvent representing this domain boundary.
    bound : float, str
        A float specifying the value (isosurface value or name) of the interpolated scalar field that represents the domain boundary.
    gt : float | curlew.geology.geoevent.GeoEvent | list
        A float or a list of floats or `curlew.geology.geoevent.GeoEvent` instances that define the scalar field(s) used to populate the hangingwall (val > bound) of this domain boundary. In essence these
        define the geological sub-model used on the hangingwall side of the domain boundary. If a float is provided, it will be populated with a constant value.
    lt : float | curlew.geology.geoevent.GeoEvent | list
        A float or a list of floats or `curlew.geology.geoevent.GeoEvent` instances that define the scalar field(s) used to populate the footwall (val < bound) of this domain boundary. In essence these
        define the geological sub-model used on the footwall side of the domain boundary. If a float is provided, it will be populated with a constant value.
    mode : str
        The overprinting mode. Options are:
            - `"above"`: replace all regions greater than the provided threshold). Useful for e.g., unconformities.
            - `"below"`: replace all regions less than than the provided threshold). Useful for e.g., intrusions.
    lithoSharpness : float, optional
        Sigmoid sharpness for soft lithology at the domain boundary (see
        :class:`~curlew.geology.interactions.Overprint`). Default is 1.0.
    structureSharpness : float, optional
        Sigmoid sharpness for ``softStructureID`` at the domain boundary. Default is 1.0.
    
    Keywords
    ----------
    All keywords are passed to `curlew.geology.geoevent.GeoEvent.__init__(...)`, many of which are then used to initialise the 
    underlying analytical or neural field. See `curlew.geology.geoevent.GeoEvent.__init__(...)` for further details.

    Returns
    ---------
    A `curlew.geology.geoevent.GeoEvent` instance for the created domain boundary. This can then be included in a `curlew.geology.geomodel.GeoModel` class. Note that the submodels (i.e. `gt` and `lt`) are now
    associated with this `curlew.geology.geoevent.GeoEvent` instance, so do not need to be (directly) passed to the `curlew.geology.geomodel.GeoModel` constructor.
    """

    # create field representing domain boundary
    o = Overprint(
        threshold=bound, mode=mode,
        lithoSharpness=lithoSharpness, structureSharpness=structureSharpness,
    )
    f = _initF( name, C=C, overprint=o, **kwargs)

    # define parent / child relationships for domain boundary field
    if not (isinstance(gt, list) or isinstance(gt, list)):
        gt = [gt]
    if not (isinstance(lt, list) or isinstance(lt, list)):
        lt = [lt]
    if isinstance(gt[-1], GeoEvent):
        gt[-1].child = f
    if isinstance(lt[-1], GeoEvent):
        lt[-1].child = f
    f.parent = gt[-1]
    f.parent2 = lt[-1]

    # build linked list for gt and lt domains
    _linkE( gt + [f])
    _linkE( lt + [f])
    
    return f

Create a GeoEvent representing a domain boundary, in which different sub-models are modelled on either side of the boundary. This can be very useful for modelling e.g., sedimentary basins, where the basin fill is modelled on one side of the boundary (onlap relations), or for modelling domain boundary faults (in which there is no known or meaningful relationship between the fault's hangingwall and footwall).

Parameters

name : str
A name for the created domain boundary.
C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
Either a pre-constructed neural field, explicit field or a set of constraints to use to construct a new GeoEvent representing this domain boundary.
bound : float, str
A float specifying the value (isosurface value or name) of the interpolated scalar field that represents the domain boundary.
gt : float | curlew.geology.geoevent.GeoEvent | list
A float or a list of floats or GeoEvent instances that define the scalar field(s) used to populate the hangingwall (val > bound) of this domain boundary. In essence these define the geological sub-model used on the hangingwall side of the domain boundary. If a float is provided, it will be populated with a constant value.
lt : float | curlew.geology.geoevent.GeoEvent | list
A float or a list of floats or GeoEvent instances that define the scalar field(s) used to populate the footwall (val < bound) of this domain boundary. In essence these define the geological sub-model used on the footwall side of the domain boundary. If a float is provided, it will be populated with a constant value.
mode : str
The overprinting mode. Options are: - "above": replace all regions greater than the provided threshold). Useful for e.g., unconformities. - "below": replace all regions less than than the provided threshold). Useful for e.g., intrusions.
lithoSharpness : float, optional
Sigmoid sharpness for soft lithology at the domain boundary (see :class:~curlew.geology.interactions.Overprint). Default is 1.0.
structureSharpness : float, optional
Sigmoid sharpness for softStructureID at the domain boundary. Default is 1.0.

Keywords

All keywords are passed to GeoEvent(…), many of which are then used to initialise the underlying analytical or neural field. See GeoEvent(…) for further details.

Returns

A GeoEvent instance for the created domain boundary. This can then be included in a GeoModel class. Note that the submodels (i.e. gt and lt) are now associated with this GeoEvent instance, so do not need to be (directly) passed to the GeoModel constructor.

def fault(name,
*,
C=None,
shortening=None,
learn_sigma=False,
contact=0,
offset=0,
width=0,
modifier=None,
n_steps=2,
dt=-1.0,
**kwargs)
Expand source code
def fault(name, *, C=None, shortening=None, learn_sigma=False, contact=0, offset=0, width=0, modifier=None, n_steps=2, dt=-1.0, **kwargs):
    """
    Create a GeoEvent representing a fault, shear zone or (optionally) dilatant shear vein.

    Parameters
    -----------
    name : str
        A name for the created stratigraphic series (and GeoEvent that represents it).
    C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
        The constraints or predefined field used to constrain this geological structure.
    shortening : np.ndarray
        A numpy array of shape (n,) defining the principal compressive stress vector. This
        is used to resolve the slip direction, by projection onto the tangent of the fault's
        interpolated scalar field. Defaults to vertical.
    learn_sigma : bool
        True if shortening should be converted to a learnable parameter. Default is False.
    contact : float | str | list
        The value (or name) defining the fault (iso)surface. Default is zero. Multi-faults (i.e. parallel
        faults with different displacements and locations but with geometry defined by a single implicit
        function) can be constructed by making `contact`, `offset` and `width` a list with one entry for
        each fault surface.
    offset : float | tuple | list
        The mode-2 slip on this (shear) fault. If a float is passed, 
        the offset will be fixed. If a tuple containing (initial, minimum, maximum) is
        passed then the offset will be learned, by constrained to the specified range.
    width : float | tuple | list
        The width of ductile deformation associated with this fault. If a float is passed
        then this will specify the width of ductile deformation (using a sigmoid function),
        such that large widths can be used for shear zones and small values of width
        used for brittle faults. The width is converted to scale factors for a sigmoid fuction using the 
        following formula: `displacementM = slipM x sigmoid(s x (4/width)) x 0.5`.
        
        Optionally, a tuple can also be passed to combine two widths, one for the fault "core"
        and one for broader ductile deformation (e.g., drag folds). This tuple should contain:
        `(width_core, width_ductile, proportion)`, where `width_core` defines the width of the 
        fault core, `width_ductile` defines the (larger) width of surrounding ductile deformation,
        and `proportion` defines the partioning of the total offset between these two deformation
        types.
    modifier : str | list, optional
        Name of a field on the fault GeoEvent (e.g. from ``addField``) evaluated at each
        point and used to scale fault slip (e.g. an ellipsoidal patch field). For
        multi-faults, pass a list with one modifier name per contact (or ``None`` entries where no scaling is needed).
    n_steps : int
        Number of explicit Euler substeps for integrating fault slip (see :class:`~curlew.geology.interactions.FaultOffset`). Default is 2.
    dt : float
        Time step per Euler substep. Default is -1.0 (i.e. reconstruct from modern to paleo-coords).

    Keywords
    ----------
    All keywords are passed to `curlew.geology.geoevent.GeoEvent.__init__(...)`, many of which are then used to initialise the 
    underlying analytical or neural field. See `curlew.geology.geoevent.GeoEvent.__init__(...)` for further details.

    Returns
    ---------
    A `curlew.geology.geoevent.GeoEvent` instance for the created structure.
    """
    # shortening
    if shortening is None: 
        shortening = np.zeros( curlew.default_dim )
        shortening[-1] = -1 # default is vertical vector
    shortening = _tensor( shortening, dev=curlew.device, dt=curlew.dtype) 
    
    if isinstance(contact, list) or isinstance(contact, tuple) or isinstance(contact, np.ndarray):
        if not (isinstance(offset, list) or isinstance(offset, tuple) or isinstance(offset, np.ndarray)):
            offset = [offset for i in contact] # listify
        if not (isinstance(width, list) or isinstance(width, tuple) or isinstance(width, np.ndarray)):
            width = [width for i in contact] # listify
        if modifier is not None and not (
            isinstance(modifier, list) or isinstance(modifier, tuple) or isinstance(modifier, np.ndarray)
        ):
            modifier = [modifier for i in contact]
        assert len(offset) == len(contact)
        assert len(width) == len(contact)
        if modifier is not None:
            assert len(modifier) == len(contact)
    else:
        contact = [contact] # listify
        offset = [offset]
        width = [width]
        if modifier is not None:
            modifier = [modifier]

    if modifier is None:
        modifier = [None] * len(contact)

    # build offset object(s)
    O = []
    for _c, _o, _w, _m in zip(contact, offset, width, modifier):
        offs = FaultOffset(
            shortening=shortening,
            offset=_o,
            contact=_c,
            width=_w,
            modifier=_m,
            n_steps=n_steps,
            dt=dt,
        )

        # handle constant or learnable offsets and/or slip direction
        init=False
        if isinstance( _o, tuple):
            _o, smin, smax = _o # shear (mode II) offset
            offs.offset = nn.Parameter( _tensor( _o, dev=curlew.device, dt=curlew.dtype) ) # will now be changed by optimiser!
            offs.offsetRange = (smin, smax) # use min and max to constrain the range of values allowed
            init=True
        else:
            offs.offset = _tensor(_o, dev=curlew.device, dt=curlew.dtype)
        if learn_sigma:
            offs.shortening = nn.Parameter( offs.shortening )
            init=True
        
        if init: # initialise the optimiser to include the new parameters
            offs.init_optim(lr=kwargs.get('learning_rate', 1e-1))
        O.append(offs)
    
    if len(O) == 1:
        O = O[0] # un-listify for single fault surfaces
    
    # build field
    f = _initF( name, C=C, deformation=O, **kwargs)

    return f

Create a GeoEvent representing a fault, shear zone or (optionally) dilatant shear vein.

Parameters

name : str
A name for the created stratigraphic series (and GeoEvent that represents it).
C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
The constraints or predefined field used to constrain this geological structure.
shortening : np.ndarray
A numpy array of shape (n,) defining the principal compressive stress vector. This is used to resolve the slip direction, by projection onto the tangent of the fault's interpolated scalar field. Defaults to vertical.
learn_sigma : bool
True if shortening should be converted to a learnable parameter. Default is False.
contact : float | str | list
The value (or name) defining the fault (iso)surface. Default is zero. Multi-faults (i.e. parallel faults with different displacements and locations but with geometry defined by a single implicit function) can be constructed by making contact, offset and width a list with one entry for each fault surface.
offset : float | tuple | list
The mode-2 slip on this (shear) fault. If a float is passed, the offset will be fixed. If a tuple containing (initial, minimum, maximum) is passed then the offset will be learned, by constrained to the specified range.
width : float | tuple | list

The width of ductile deformation associated with this fault. If a float is passed then this will specify the width of ductile deformation (using a sigmoid function), such that large widths can be used for shear zones and small values of width used for brittle faults. The width is converted to scale factors for a sigmoid fuction using the following formula: displacementM = slipM x sigmoid(s x (4/width)) x 0.5.

Optionally, a tuple can also be passed to combine two widths, one for the fault "core" and one for broader ductile deformation (e.g., drag folds). This tuple should contain: (width_core, width_ductile, proportion), where width_core defines the width of the fault core, width_ductile defines the (larger) width of surrounding ductile deformation, and proportion defines the partioning of the total offset between these two deformation types.

modifier : str | list, optional
Name of a field on the fault GeoEvent (e.g. from addField) evaluated at each point and used to scale fault slip (e.g. an ellipsoidal patch field). For multi-faults, pass a list with one modifier name per contact (or None entries where no scaling is needed).
n_steps : int
Number of explicit Euler substeps for integrating fault slip (see :class:~curlew.geology.interactions.FaultOffset). Default is 2.
dt : float
Time step per Euler substep. Default is -1.0 (i.e. reconstruct from modern to paleo-coords).

Keywords

All keywords are passed to GeoEvent(…), many of which are then used to initialise the underlying analytical or neural field. See GeoEvent(…) for further details.

Returns

A GeoEvent instance for the created structure.

def fold(name,
*,
origin,
compression,
extension,
wavelength,
amplitude=1.0,
sharpness=1.0,
estStrain=False)
Expand source code
def fold( name, *, origin, compression, extension, wavelength, amplitude=1.0, sharpness=1.0, estStrain=False):
    """
    Create a GeoEvent representing a fold structure.

    This constructs an explicitly defined scalar field aligned with the
    fold-shortening axis, then computes displacement vectors that reproduce a
    fold geometry with the specified wavelength, amplitude, and sharpness.
    Optionally, an estimated finite strain required to “undo” the folding can
    be computed from the waveform geometry (and applied when transforming from
    model to paleo-coordinates).

    Parameters
    ------------
    name : str
        A name for the created stratigraphic series (and GeoEvent that represents it).
    origin : array-like
        A point through which the fold scalar field passes; used as the field’s origin.
    compression : array-like
        Vector describing the principal compression direction. This is
        normalized and scaled internally to represent the fold wavelength. This vector
        should be perpendicular to the fold's axial foliation.
    extension : array-like
        Vector describing the principal extension direction, such that the fold axis
        is the cross product between the extension and compression vectors.
    wavelength : float
        The fold wavelength. Used to scale the compression vector into a
        periodic scalar field.
    amplitude : float, default=1.0
        Amplitude of the fold waveform.
    sharpness : float, default=1.0
        Sharpness of the fold waveform, controlling the peakedness of the
        blended wave.
    estStrain : bool, default=False
        If True, numerically estimate the finite strain required to restore the
        folded layer to its unfolded length using a line integral of the
        waveform.

    Keywords
    ----------
    All keywords are passed to `curlew.geology.geoevent.GeoEvent.__init__(...)`, many of which are then used to initialise the 
    underlying analytical or neural field. See `curlew.geology.geoevent.GeoEvent.__init__(...)` for further details.

    Returns
    ---------
    A `curlew.geology.geoevent.GeoEvent` instance representing the fold structure, with an
    associated analytical scalar field and deformation function.
    """
    # create a fold field and associated deformation function
    compression = compression / np.linalg.norm(compression) # direction of principal compression
    compression *= (2 * np.pi) / wavelength  # scale normal vector to give appropriate wavelength 
    fa = LinearField( name, input_dim=len(compression), 
              origin=origin, gradient=compression ) # TODO; allow also interpolated fields here
    
    # evaluate fold strain
    if estStrain:
        x = _tensor( np.linspace(0,2,1000), dev=curlew.device, dt=curlew.dtype )
        y = blended_wave(x, f=sharpness, A=amplitude, T=2) # evaluate one waveform
        dx = torch.mean( torch.diff(x) )
        dy = torch.diff( y )
        l0 = torch.sum( torch.sqrt( dx**2 + dy**2 ) ) # line-integral gives initial length
        l1 = x[-1] - x[0] # current length is known
        strain = ((l0-l1) / l1).item() # hence get strain needed to undo folding
    else:
        strain = 0 # ignore shortening
    
    # create a lambda function for evaluating folds from scalar value
    f = lambda x: blended_wave( x, f=sharpness, A=amplitude, T=2)

    # creat fold object
    extension = extension / np.linalg.norm(extension) # principal stretching direction
    defo = FoldOffset( thicker=extension, shorter=compression, shortening=strain, periodic=f )

    # create and return our GeoEvent instance
    return GeoEvent( name, None, field=fa, # use pre-existing (analytical) field rather than creating a new one
                    deformation=defo )

Create a GeoEvent representing a fold structure.

This constructs an explicitly defined scalar field aligned with the fold-shortening axis, then computes displacement vectors that reproduce a fold geometry with the specified wavelength, amplitude, and sharpness. Optionally, an estimated finite strain required to “undo” the folding can be computed from the waveform geometry (and applied when transforming from model to paleo-coordinates).

Parameters

name : str
A name for the created stratigraphic series (and GeoEvent that represents it).
origin : array-like
A point through which the fold scalar field passes; used as the field’s origin.
compression : array-like
Vector describing the principal compression direction. This is normalized and scaled internally to represent the fold wavelength. This vector should be perpendicular to the fold's axial foliation.
extension : array-like
Vector describing the principal extension direction, such that the fold axis is the cross product between the extension and compression vectors.
wavelength : float
The fold wavelength. Used to scale the compression vector into a periodic scalar field.
amplitude : float, default=1.0
Amplitude of the fold waveform.
sharpness : float, default=1.0
Sharpness of the fold waveform, controlling the peakedness of the blended wave.
estStrain : bool, default=False
If True, numerically estimate the finite strain required to restore the folded layer to its unfolded length using a line integral of the waveform.

Keywords

All keywords are passed to GeoEvent(…), many of which are then used to initialise the underlying analytical or neural field. See GeoEvent(…) for further details.

Returns

A GeoEvent instance representing the fold structure, with an associated analytical scalar field and deformation function.

def restore(name,
*,
velocity,
C=None,
H=None,
base=-inf,
mode='above',
onlap=False,
lithoSharpness=1.0,
structureSharpness=1.0,
constraints=None,
fault_lift=None,
**kwargs)
Expand source code
def restore(name, *, velocity, C=None, H=None, base=-np.inf, mode="above", onlap=False,
            lithoSharpness=1.0, structureSharpness=1.0, constraints=None,
            fault_lift=None, **kwargs):
    """
    Create a GeoEvent for diffeomorphic restoration (fold / structure inversion).

    This is a :func:`strati`-style **generative** event — it gets the same
    :class:`~curlew.geology.interactions.Overprint` (``base``/``mode``/``onlap``)
    and thus the same isosurface/lithology behaviour as any other stratigraphic
    package — just parameterised by a different kind of interpolator: a
    :class:`~curlew.fields.restoration.RestorationField`, whose two (2-D: one,
    3-D: two) Clebsch potentials produce a divergence-free *velocity* field.
    That velocity is integrated into a diffeomorphism Φ; the event's scalar
    value is the restored depth ``φ(x) = (Φ⁻¹x)[depth_axis]``, and losses
    compare bedding orientation/values *after* integrating through Φ (Nanson's
    relation) to observed orientations/values, rather than fitting the scalar
    field directly (see :meth:`~curlew.fields.restoration.RestorationField.loss`).

    The event's ``deformation`` is set to the field's own
    :class:`~curlew.geology.interactions.FlowOffset` integrator
    (``field.integrator``) rather than a separately constructed one, so younger
    events are undeformed through the *same* Φ (and parameters) used for scalar
    evaluation and lithology.

    Parameters
    ----------
    name : str
        Event name (also used for the :class:`RestorationField`).
    velocity : :class:`~curlew.fields.clebsch.ClebschVelocity` or compatible module
        Pre-built divergence-free velocity field (e.g. from
        :func:`~curlew.fields.clebsch.clebsch_velocity` or
        :func:`~curlew.fields.clebsch.temporal_clebsch_velocity`).  Must expose
        ``dim``, ``forward``, and ``forward_and_jacobian``.  When using fault
        lifts, build with ``lift_dim`` matching :attr:`~curlew.fields.lift.FaultLift.lift_dim`.
    C : :class:`~curlew.core.CSet`, optional
        Constraints in **modern** coordinates.  Interpreted as:

        - ``gp``/``gv`` → bedding normals (Nanson; weight ``H.grad_loss``)
        - ``vp``/``vv`` → restored depth values (weight ``H.value_loss``)
        - ``eq`` → horizon traces, flatness in restored depth (weight ``H.eq_loss``)
        - ``iq`` → stratigraphic ordering on **restored depth** ``(Φ⁻¹ x)[depth_axis]``
          (weight ``H.iq_loss``; same ``CSet.iq`` format as scalar fields)
        - ``gop``/``gov`` → sign-invariant normals (weight ``H.ori_loss``)
    constraints : sequence of :class:`~curlew.fields.restoration_constraints.RestorationConstraint`, optional
        Pluggable priors (e.g. layer thickness, orthogonal thickness) attached
        to the :class:`RestorationField` — not stored on ``CSet``.
    fault_lift : :class:`~curlew.fields.lift.FaultLift`, optional
        Pre-built GWN sheet lift for fault-aware restoration.
    H : :class:`~curlew.core.HSet`, optional
        Loss weights.  Defaults enable gradient, value, and equality terms;
        ``iq_loss`` and scalar-field terms (monotonicity / thickness) are off
        by default.
    base, mode, onlap, lithoSharpness, structureSharpness
        Same as :func:`strati` — control the basal :class:`Overprint` for this
        (folded) stratigraphic package, and isosurface/lithology sharpness.

    Keywords
    --------
    Passed to :class:`~curlew.fields.restoration.RestorationField` via
    :class:`~curlew.geology.geoevent.GeoEvent` (e.g. ``n_steps``,
    ``depth_axis``, ``signed_normals``, ``sigma_floor``).  ``constraints`` may
    also be passed here instead of as a keyword to :func:`restore`.

    Returns
    -------
    :class:`~curlew.geology.geoevent.GeoEvent`
    """
    from curlew.core import HSet
    from curlew.fields.restoration import RestorationField

    if fault_lift is not None:
        kwargs["fault_lift"] = fault_lift
    kwargs["velocity"] = velocity

    if H is None:
        H = HSet(
            value_loss="1.0",
            grad_loss="1.0",
            eq_loss="1.0",
            mono_loss=0,
            thick_loss=0,
            flat_loss=0,
            ori_loss=0,
            iq_loss=0,
        )
    if constraints is not None:
        kwargs["constraints"] = list(constraints)
    kwargs.setdefault("type", RestorationField)
    kwargs["H"] = H
    event = strati(
        name, C=C, base=base, mode=mode, onlap=onlap,
        lithoSharpness=lithoSharpness, structureSharpness=structureSharpness,
        **kwargs,
    )
    event.deformation = event.field.integrator  # reuse the field's own Phi; no duplicate integrator
    return event

Create a GeoEvent for diffeomorphic restoration (fold / structure inversion).

This is a :func:strati()-style generative event — it gets the same :class:~curlew.geology.interactions.Overprint (base/mode/onlap) and thus the same isosurface/lithology behaviour as any other stratigraphic package — just parameterised by a different kind of interpolator: a :class:~curlew.fields.restoration.RestorationField, whose two (2-D: one, 3-D: two) Clebsch potentials produce a divergence-free velocity field. That velocity is integrated into a diffeomorphism Φ; the event's scalar value is the restored depth φ(x) = (Φ⁻¹x)[depth_axis], and losses compare bedding orientation/values after integrating through Φ (Nanson's relation) to observed orientations/values, rather than fitting the scalar field directly (see :meth:~curlew.fields.restoration.RestorationField.loss).

The event's deformation is set to the field's own :class:~curlew.geology.interactions.FlowOffset integrator (field.integrator) rather than a separately constructed one, so younger events are undeformed through the same Φ (and parameters) used for scalar evaluation and lithology.

Parameters

name : str
Event name (also used for the :class:RestorationField).
velocity : :class:~curlew.fields.clebsch.ClebschVelocityor </code>compatible module
Pre-built divergence-free velocity field (e.g. from :func:~curlew.fields.clebsch.clebsch_velocity or :func:~curlew.fields.clebsch.temporal_clebsch_velocity). Must expose dim, forward, and forward_and_jacobian. When using fault lifts, build with lift_dim matching :attr:~curlew.fields.lift.FaultLift.lift_dim.
C : :class:~curlew.core.CSet``, optional

Constraints in modern coordinates. Interpreted as:

  • gp/gv → bedding normals (Nanson; weight H.grad_loss)
  • vp/vv → restored depth values (weight H.value_loss)
  • eq → horizon traces, flatness in restored depth (weight H.eq_loss)
  • iq → stratigraphic ordering on restored depth (Φ⁻¹ x)[depth_axis] (weight H.iq_loss; same CSet.iq format as scalar fields)
  • gop/gov → sign-invariant normals (weight H.ori_loss)
constraints : sequence of :class:~curlew.fields.restoration_constraints.RestorationConstraint``, optional
Pluggable priors (e.g. layer thickness, orthogonal thickness) attached to the :class:RestorationField — not stored on CSet.
fault_lift : :class:~curlew.fields.lift.FaultLift``, optional
Pre-built GWN sheet lift for fault-aware restoration.
H : :class:~curlew.core.HSet``, optional
Loss weights. Defaults enable gradient, value, and equality terms; iq_loss and scalar-field terms (monotonicity / thickness) are off by default.
base, mode, onlap, lithoSharpness, structureSharpness
Same as :func:strati() — control the basal :class:Overprint for this (folded) stratigraphic package, and isosurface/lithology sharpness.

Keywords

Passed to :class:~curlew.fields.restoration.RestorationField via :class:~curlew.geology.geoevent.GeoEvent (e.g. n_steps, depth_axis, signed_normals, sigma_floor). constraints may also be passed here instead of as a keyword to :func:restore().

Returns

:class:~curlew.geology.geoevent.GeoEvent

def sheet(name,
*,
C=None,
contact=(-1, 1),
aperture=2,
n_steps=1,
dt=-1.0,
lithoSharpness=1.0,
structureSharpness=1.0,
**kwargs)
Expand source code
def sheet(name, *, C=None, contact=(-1,1), aperture=2, n_steps=1, dt=-1.0,
          lithoSharpness=1.0, structureSharpness=1.0, **kwargs):
    """
    Create a GeoEvent representing a sheet intrusion (dyke, sill or vein).

    Parameters
    -----------
    name : str
        A name for the created stratigraphic series (and GeoEvent that represents it).
    C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
        The constraints or predefined field used to constrain this geological structure.
    contact : tuple | list [ tuple ]
        A tuple of floats specifying the scalar values for the (upper, lower) sides of 
        the intrusion, or a list of (upper, lower) tuples if multiple dykes are defined.
    aperture : float | list
        The aperture (Mode I opening) of the dyke, or a list thereof. Used to displace surrounding rocks.
    n_steps : int
        Number of Euler substeps for integrating sheet-parallel displacement (see :class:`~curlew.geology.interactions.SheetOffset`). Default is 1.
    dt : float
        Time step per Euler substep. Default is -1.0 (i.e. reconstruct from modern to paleo-coords).
    lithoSharpness : float, optional
        Sigmoid sharpness passed to each interval :class:`~curlew.geology.interactions.Overprint`
        (and used for soft lithology on this event). Default is 1.0.
    structureSharpness : float, optional
        Sigmoid sharpness for ``softStructureID`` on each interval overprint. Default is 1.0.

    Keywords
    ----------
    All keywords are passed to `curlew.geology.geoevent.GeoEvent.__init__(...)`, many of which are then used to initialise the 
    underlying analytical or neural field. See `curlew.geology.geoevent.GeoEvent.__init__(...)` for further details.

    Returns
    ---------
    A `curlew.geology.geoevent.GeoEvent` instance for the created structure.
    """
    if isinstance(contact, float) or isinstance(contact, int):
        contact = (-contact, contact) # define as upper and lower surface (assuming symmetry)
    if not isinstance(contact, list):
        contact = [contact]
    if not isinstance(aperture, list):
        aperture = [aperture for c in contact]
    assert len(contact) == len(aperture)
    
    offset = []
    overprint = []
    for _c, _a in zip(contact, aperture):
        offset.append(
            SheetOffset(contact=_c, aperture=_a, n_steps=n_steps, dt=dt)
        )
        overprint.append(
            Overprint(
                threshold=_c, mode='in',
                lithoSharpness=lithoSharpness, structureSharpness=structureSharpness,
            )
        )
    if len(offset) == 1:
        offset = offset[0]
        overprint = overprint[0]

    return _initF( name, C=C, deformation=offset, overprint=overprint, **kwargs)

Create a GeoEvent representing a sheet intrusion (dyke, sill or vein).

Parameters

name : str
A name for the created stratigraphic series (and GeoEvent that represents it).
C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
The constraints or predefined field used to constrain this geological structure.
contact : tuple | list [ tuple ]
A tuple of floats specifying the scalar values for the (upper, lower) sides of the intrusion, or a list of (upper, lower) tuples if multiple dykes are defined.
aperture : float | list
The aperture (Mode I opening) of the dyke, or a list thereof. Used to displace surrounding rocks.
n_steps : int
Number of Euler substeps for integrating sheet-parallel displacement (see :class:~curlew.geology.interactions.SheetOffset). Default is 1.
dt : float
Time step per Euler substep. Default is -1.0 (i.e. reconstruct from modern to paleo-coords).
lithoSharpness : float, optional
Sigmoid sharpness passed to each interval :class:~curlew.geology.interactions.Overprint (and used for soft lithology on this event). Default is 1.0.
structureSharpness : float, optional
Sigmoid sharpness for softStructureID on each interval overprint. Default is 1.0.

Keywords

All keywords are passed to GeoEvent(…), many of which are then used to initialise the underlying analytical or neural field. See GeoEvent(…) for further details.

Returns

A GeoEvent instance for the created structure.

def stock(name, *, C, H, contact=0, **kwargs)
Expand source code
def stock(name, *, C, H, contact=0, **kwargs):
    """
    Create a GeoEvent representing a stock, pluton or batholith. 
    """
    pass

Create a GeoEvent representing a stock, pluton or batholith.

def strati(name,
*,
C=None,
base=-inf,
mode='above',
onlap=False,
lithoSharpness=1.0,
structureSharpness=1.0,
**kwargs)
Expand source code
def strati( name, *, C=None, base = -np.inf, mode="above", onlap=False,
            lithoSharpness=1.0, structureSharpness=1.0, **kwargs):
    """
    Create a GeoEvent representing a stratigraphic series (base stratigraphy or unconformity).

    Parameters
    ------------
    name : str
        A name for the created stratigraphic series (and GeoEvent that represents it).
    C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
        Either a pre-constructed neural field, explicit field or a set of 
        constraints to use to construct a new GeoEvent representing this stratigraphic series.
    base : float | str
        Scalar value or isosurface name representing the basement contact of this stratigraphic series. This
        is important for unconformities, as only values greater than `base` are considered
        part of this event. Default is -infinity. If a string is provided, an isosurface with the corresponding
        name *must* be added to the returned GeoEvent.
    mode : str
        The overprinting mode. Options are:
            - `"above"`: replace all regions greater than the provided threshold). Useful for e.g., unconformities.
            - `"below"`: replace all regions less than than the provided threshold). Useful for e.g., intrusions.
    onlap : bool
        True if the geometry of the basal unconformity for the created stratigraphic package should be defined
        by the parent (older) field rather than the created one. If True, bedding generated by this event will
        onlap onto the basal unconformity.
    lithoSharpness : float, optional
        Sigmoid sharpness for ``softLithoID`` at the basal overprint and for isosurface lithology assignment
        on this event (see :class:`~curlew.geology.interactions.Overprint`). Default is 1.0.
    structureSharpness : float, optional
        Sigmoid sharpness for ``softStructureID`` at the basal overprint. Default is 1.0.
        
    Keywords
    ----------
    All keywords are passed to `curlew.geology.geoevent.GeoEvent.__init__(...)`, many of which are then used to initialise the 
    underlying analytical or neural field. See `curlew.geology.geoevent.GeoEvent.__init__(...)` for further details.

    Returns
    ---------
    A `curlew.geology.geoevent.GeoEvent` instance for the created structure.
    """
    o = Overprint(
        threshold=base, mode=mode,
        lithoSharpness=lithoSharpness, structureSharpness=structureSharpness,
    )
    if onlap:
        o.defaultDomain = 'parent' # older unit defines unconformity geometry
    return _initF( name, C=C, overprint=o, **kwargs)

Create a GeoEvent representing a stratigraphic series (base stratigraphy or unconformity).

Parameters

name : str
A name for the created stratigraphic series (and GeoEvent that represents it).
C : curlew.core.CSet | curlew.fields.BaseAF | curlew.fields.BaseNF
Either a pre-constructed neural field, explicit field or a set of constraints to use to construct a new GeoEvent representing this stratigraphic series.
base : float | str
Scalar value or isosurface name representing the basement contact of this stratigraphic series. This is important for unconformities, as only values greater than base are considered part of this event. Default is -infinity. If a string is provided, an isosurface with the corresponding name must be added to the returned GeoEvent.
mode : str
The overprinting mode. Options are: - "above": replace all regions greater than the provided threshold). Useful for e.g., unconformities. - "below": replace all regions less than than the provided threshold). Useful for e.g., intrusions.
onlap : bool
True if the geometry of the basal unconformity for the created stratigraphic package should be defined by the parent (older) field rather than the created one. If True, bedding generated by this event will onlap onto the basal unconformity.
lithoSharpness : float, optional
Sigmoid sharpness for softLithoID at the basal overprint and for isosurface lithology assignment on this event (see :class:~curlew.geology.interactions.Overprint). Default is 1.0.
structureSharpness : float, optional
Sigmoid sharpness for softStructureID at the basal overprint. Default is 1.0.

Keywords

All keywords are passed to GeoEvent(…), many of which are then used to initialise the underlying analytical or neural field. See GeoEvent(…) for further details.

Returns

A GeoEvent instance for the created structure.