speed.py

LevelSetSpeedMethods dataclass

LevelSetSpeedMethods(
    Δx: float,
    Δt: float = 0,
    n_total: int | None = None,
    i_step: int | None = None,
    t_total: float = 0,
    Δt_reinitialize: float = 0.1,
    κ_boost: float = 1,
    n_slabfailure: int | None = None,
    n_pad_pixels: int | None = None,
    band_width: float | None = None,
    fm_order: int | None = None,
    n_pixels_xz0: tuple[int, int] | None = None,
    n_pixels_xz: tuple[int, int] | None = None,
    do_extend_line: bool = True,
    do_φ_everywhere: bool | None = None,
    domain: Domain | None = None,
    raw_domain: Domain | None = None,
    surface: (
        StraightLineSurface
        | ErrorFunctionSurface
        | SemiCircularArcSurface
        | QuarterCircularArcSurface
        | None
    ) = None,
    substrate: (
        OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate | None
    ) = None,
    model: Isotropic | Step | ExponentialActivation | None = None,
    gma: GeometricMechanicsAnalysis | None = None,
)

              flowchart TD
              erosionfront.levelset.speed.LevelSetSpeedMethods[LevelSetSpeedMethods]
              erosionfront.levelset.grid.LevelSetGridMethods[LevelSetGridMethods]
              erosionfront.levelset.elementary.LevelSetElementaryMethods[LevelSetElementaryMethods]
              erosionfront.levelset.data.LevelSetData[LevelSetData]

                              erosionfront.levelset.grid.LevelSetGridMethods --> erosionfront.levelset.speed.LevelSetSpeedMethods
                                erosionfront.levelset.elementary.LevelSetElementaryMethods --> erosionfront.levelset.grid.LevelSetGridMethods
                                erosionfront.levelset.data.LevelSetData --> erosionfront.levelset.elementary.LevelSetElementaryMethods
                




              click erosionfront.levelset.speed.LevelSetSpeedMethods href "" "erosionfront.levelset.speed.LevelSetSpeedMethods"
              click erosionfront.levelset.grid.LevelSetGridMethods href "" "erosionfront.levelset.grid.LevelSetGridMethods"
              click erosionfront.levelset.elementary.LevelSetElementaryMethods href "" "erosionfront.levelset.elementary.LevelSetElementaryMethods"
              click erosionfront.levelset.data.LevelSetData href "" "erosionfront.levelset.data.LevelSetData"
            

Speed-related methods for use in erosion front modeling.

Subclasses LevelSetGridMethods.

Attributes

φ: MaskedArray Signed-distance grid at this time step. φ_next: MaskedArray | None Signed-distance grid at next time step.

Methods:

• __post_init__ –

  Complete instantiation by computing domain, pixel size, offset.

• above_or_below_ji –

  Test whether pixel i,j lies above or below surface.

• above_or_below_xz –

  Test whether point x,z lies above or below surface.

• compute_contour –

  Generate a continuous (for now) contour of grid function f.

• compute_domain –

  (Re)build domain after resizing.

• compute_line_pixel_offsets –

  Measure signed offsets between rasterized line pixels and line curve.

• compute_pixel_dimensions –

  Calculate numbers of pixels spanning the grid.

• compute_pixel_origin –

  Calculate pixel offset of domain x,z position (0,0).

• fm_distance –

  Use fast marching to compute signed distance from a ref grid (zeroset).

• fm_distance_everywhere –

  FM compute a narrowband signed distance and then extend everywhere.

• get_extent –

  Get domain (x,z) extent.

• get_front_points –

  Fetch points along erosion front surface aka zeroset of f.

• get_line_points –

  Convert Shapely geometry Line into (x,z) points.

• get_points_xz –

  Get surface (x,z) coordinates.

• get_substrate_ξ0 –

  Generate relative η grid (effective speed ξ0).

• get_x_limits –

  Get domain x min, max.

• get_z_limits –

  Get domain z min, max.

• get_φ_next –

  Return next signed-distance grid φ_next or just φ if not available.

• ij_to_xz –

  Convert pixel indices [i,j] into pixel center coords (x,z).

• ji_to_xz –

  Convert reversed pixel indices [j,i] into pixel center coords (x,z).

• make_line –

  Convert surface (x,z) coordinates into a Shapely Line geometry.

• make_polygon –

  Create a 'bounding' polygon from surface x,z points.

• measure_normal_distance –

  Compute distance between line and point using Shapely.

• measure_φ_ji –

  Compute signed distance between line and i,j pixel using Shapely.

• measure_φ_xz –

  Compute signed distance between line and xz point using Shapely.

• pad_domain –

  Expand the domain vertically by pixel width amounts.

• pad_grid –

  Pad grid with specified number of pixels below and above.

• rasterize_line –

  Rasterize surface line as x,z points into grid pixels using Rasterio.

• redo_φ_everywhere –

  Reassign signed distance values to pixels beyond the narrow band.

• xz_to_ij –

  Convert pixel center coords (x,z) into pixel indices [i,j].

__post_init__

__post_init__()

Complete instantiation by computing domain, pixel size, offset.

Attributes

raw_domain: Domain Model domain instance. n_pixels_xz: tuple[int, int]: Grid pixel size. n_pixels_xz0: tuple[int, int]: Origin location in pixels.

Source code in erosionfront/levelset/elementary.py

def __post_init__(self):
    """
    Complete instantiation by computing domain, pixel size, offset.

    Attributes
    ----------
    raw_domain: Domain
        Model domain instance.
    n_pixels_xz: tuple[int, int]:
        Grid pixel size.
    n_pixels_xz0: tuple[int, int]:
        Origin location in pixels.
    """
    super().__post_init__()
    # Compute bounds, extent, pixel dimensions of grid.
    # These methods MUST be provided by any subclass.
    if self.surface is None:
        raise ValueError("Surface not specified")
    self.raw_domain \
        = self.compute_domain(
            self.surface.x, self.surface.z, self.Δx,
        )
    self.n_pixels_xz \
        = self.compute_pixel_dimensions(self.raw_domain, self.Δx,)
    self.n_pixels_xz0 \
        = self.compute_pixel_origin(self.raw_domain, self.Δx,)

above_or_below_ji

above_or_below_ji(ji: tuple[int, int] | NDArray, in_rock_polygon: Polygon) -> int

Test whether pixel i,j lies above or below surface.

Uses "in_rock_polygon" to test if pixel is within the bedrock or in the air. In other words, it tests which side of the surface zeroset the pixel lies.

The pixel center is used in this computation.

Parameters

ji: tuple[float,float] | NDArray Test pixel with reversed indexes. in_rock_polygon: Polygon Shapely polygon defining "below ground".

Atrributes

Uses ji_to_xz method.

Returns

-1 if "in the air" and +1 if "in the rock".

Source code in erosionfront/levelset/grid.py

def above_or_below_ji(
        self,
        ji: tuple[int,int] | NDArray,
        in_rock_polygon: Polygon,
    ) -> int:
    """
    Test whether pixel i,j lies above or below surface.

    Uses "in_rock_polygon" to test if pixel is within the bedrock or in 
    the air. In other words, it tests which side of the surface zeroset
    the pixel lies.

    The pixel center is used in this computation.

    Parameters
    ----------
    ji: tuple[float,float] | NDArray
        Test pixel with reversed indexes.
    in_rock_polygon: Polygon
        Shapely polygon defining "below ground".

    Atrributes
    ----------
    Uses ji_to_xz method.

    Returns
    -------
    -1 if "in the air" and +1 if "in the rock".
    """
    # integer pixel indices convert to pixel *center* x,z coords
    return (
        -1 if shapely.contains_xy(in_rock_polygon, *self.ji_to_xz(ji)) 
        else +1
    )

above_or_below_xz staticmethod

above_or_below_xz(xz: tuple[float, float] | NDArray, in_rock_polygon: Polygon) -> int

Test whether point x,z lies above or below surface.

Uses "in_rock_polygon" to test if point is within the bedrock or in the air. In other words, it tests which side of the surface zeroset the point lies.

Parameters

xz: tuple[float,float] | NDArray Test point. in_rock_polygon: Polygon Shapely polygon defining "below ground".

Returns

-1 if "in the air" and +1 if "in the rock".

Source code in erosionfront/levelset/grid.py

@staticmethod
def above_or_below_xz(
        xz: tuple[float,float] | NDArray,
        in_rock_polygon: Polygon,
    ) -> int:
    """
    Test whether point x,z lies above or below surface.

    Uses "in_rock_polygon" to test if point is within the bedrock or in 
    the air. In other words, it tests which side of the surface zeroset
    the point lies.

    Parameters
    ----------
    xz: tuple[float,float] | NDArray
        Test point.
    in_rock_polygon: Polygon
        Shapely polygon defining "below ground".

    Returns
    -------
    -1 if "in the air" and +1 if "in the rock".
    """
    return (
        -1 if shapely.contains_xy(in_rock_polygon, *xz,) else +1
    )

compute_contour

compute_contour(f: NDArray | MaskedArray, level: float = 0) -> NDArray

Generate a continuous (for now) contour of grid function f.

BUG: the contouring function actually returns a set of disjoint contour segments, often with more than one segment. Here it's incorrectly assumed these segments are contiguous and can simply be concatenated. This is not true! But to do better downstream use of this method needs to handle the segment set properly, which is going to be inconvenient.

Parameters

f: NDArray | MaskedArray Smooth grid (e.g., signed distance) on which to compute contour. level: float=0 Chosen contour level.

Attributes

Uses ij_to_xz method.

Returns

NDArray Single sequence of x,z contour coordinates as (2,N) array.

Source code in erosionfront/levelset/grid.py

def compute_contour(
        self, 
        f: NDArray | MaskedArray, 
        level: float=0,
    ) -> NDArray:
    """
    Generate a continuous (for now) contour of grid function f.

    BUG: the contouring function actually returns a set of
    disjoint contour segments, often with more than one segment.
    Here it's incorrectly assumed these segments are contiguous 
    and can simply be concatenated. This is not true! But to do
    better downstream use of this method needs to handle the segment set
    properly, which is going to be inconvenient.

    Parameters
    ----------
    f: NDArray | MaskedArray
        Smooth grid (e.g., signed distance) on which to compute contour.
    level: float=0
        Chosen contour level.

    Attributes
    ----------
    Uses ij_to_xz method.

    Returns
    -------
    NDArray
        Single sequence of x,z contour coordinates as (2,N) array.
    """
    contours: list = skimage.measure.find_contours(
        f, 
        level=level, 
        mask=~f.mask if type(f) is MaskedArray else None,
    )
    # print(len(contours))
    #HACK: drop second contour (presumably generated by front sliding
    #      off grid) if it's too short & thus likely erroneous
    # if len(contours)>1:
    #     # print(contours[0].shape, contours[1].shape)
    #     if contours[1].shape[0]<10:
    #         contours = contours[0]
    return self.ij_to_xz(np.fliplr(np.vstack(contours))).T

compute_domain staticmethod

compute_domain(x: tuple | NDArray, z: tuple | NDArray, resolution: float) -> Domain

(Re)build domain after resizing.

Parameters

x: tuple | NDArray Vectors of x coordinates for which a domain geometry needs to be measured. z: tuple | NDArray Vectors of z coordinates for which a domain geometry needs to be measured. resolution: float Grid pixel size.

Returns

Domain Domain geometry class instance.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def compute_domain(
        x: tuple | NDArray,
        z: tuple | NDArray,
        resolution: float,
    ) -> Domain:
    """
    (Re)build domain after resizing.

    Parameters
    ----------
    x: tuple | NDArray
        Vectors of x coordinates for which a domain geometry needs
        to be measured.
    z: tuple | NDArray
        Vectors of z coordinates for which a domain geometry needs
        to be measured.
    resolution: float
        Grid pixel size.

    Returns
    -------
    Domain
        Domain geometry class instance.
    """
    n_digits: int = int(round(-np.log10(resolution)+0.5))+1
    min: Callable = lambda xz: round(float(
        (np.min(np.array(xz)/resolution) - 0.5) 
    ) * resolution, n_digits)
    max: Callable = lambda xz: round(float(
        (np.max(np.array(xz)/resolution) + 0.5) 
    ) * resolution, n_digits)
    return Domain(min(x), max(x), min(z), max(z),)

compute_line_pixel_offsets

compute_line_pixel_offsets(
    line: LineString, ρ: NDArray, do_φ_everywhere: bool | None = False
) -> NDArray

Measure signed offsets between rasterized line pixels and line curve.

Spatial aliasing during rasterization means each line pixel is likely offset a bit from the line curve, so this function computes these offsets.

Parameters

line: LineString ρ: NDArray

bool

Compute distances for all grid pixels, not just those not masked.

Attributes

Uses measure_φ_ji and above_or_below_ji methods.

Returns

NDArray Grid with signed-distance normal offsets from the line curve at each rasterized line pixel.

Source code in erosionfront/levelset/grid.py

def compute_line_pixel_offsets(
        self, 
        line: LineString,
        ρ: NDArray,
        do_φ_everywhere: bool | None=False,
    ) -> NDArray:
    """
    Measure signed offsets between rasterized line pixels and line curve.

    Spatial aliasing during rasterization means each line pixel is likely
    offset a bit from the line curve, so this function computes these
    offsets.

    Parameters
    ----------
    line: LineString
    ρ: NDArray

    do_φ_everywhere: bool
        Compute distances for all grid pixels, not just those not masked.

    Attributes
    ----------
    Uses measure_φ_ji and above_or_below_ji methods.

    Returns
    -------
    NDArray
        Grid with signed-distance normal offsets from the line curve
        at each rasterized line pixel.
    """
    offsets_grid: NDArray \
        = np.zeros_like(ρ, dtype=np.float64,)
    for ji_ in np.argwhere(ρ!=0):
        # integer pixel indices convert to pixel *center* x,z coords
        offsets_grid[*ji_] = self.measure_φ_ji(line,ji_,)
    d_min: float = np.min(offsets_grid)
    d_max: float = np.max(offsets_grid)
    if do_φ_everywhere:
        for ji_ in np.argwhere(ρ==0):
            offsets_grid[*ji_] = (
                d_min 
                if self.above_or_below_ji(ji_, self.in_rock_polygon,)==-1 
                else d_max
            )
    return offsets_grid

compute_pixel_dimensions staticmethod

compute_pixel_dimensions(domain: Domain, resolution: float) -> tuple[int, int]

Calculate numbers of pixels spanning the grid.

Parameters

domain: Domain Grid domain geometry. resolution: float Grid pixel size.

Returns

tuple[int, int] Width and height of grid in pixel numbers.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def compute_pixel_dimensions(
        domain: Domain,
        resolution: float,
    ) -> tuple[int, int]:
    """
    Calculate numbers of pixels spanning the grid.

    Parameters
    ----------
    domain: Domain
        Grid domain geometry.
    resolution: float
        Grid pixel size.

    Returns
    -------
    tuple[int, int]
        Width and height of grid in pixel numbers.
    """
    n_i_pixels: int = int(
        (domain.x_max-domain.x_min)/resolution + 0.5
    )
    n_j_pixels: int = int(
        (domain.z_max-domain.z_min)/resolution + 0.5
    )
    return (n_i_pixels, n_j_pixels,)

compute_pixel_origin staticmethod

compute_pixel_origin(domain: Domain, resolution: float) -> tuple[int, int]

Calculate pixel offset of domain x,z position (0,0).

Parameters

domain: Domain Grid domain geometry. resolution: float Grid pixel size.

Returns

tuple[int, int] Pixel indexes of bottom-left origin of grid.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def compute_pixel_origin(
        domain: Domain,
        resolution: float,
    ) -> tuple[int, int]:
    """
    Calculate pixel offset of domain x,z position (0,0).

    Parameters
    ----------
    domain: Domain
        Grid domain geometry.
    resolution: float
        Grid pixel size.

    Returns
    -------
    tuple[int, int]
        Pixel indexes of bottom-left origin of grid.        
    """
    i_origin: int = int(domain.x_min/resolution+0.5)
    j_origin: int = int(domain.z_min/resolution+0.5)
    return (i_origin, j_origin,)

fm_distance staticmethod

fm_distance(
    f: NDArray, band_width: float | None, resolution: float, order: int | None
) -> MaskedArray

Use fast marching to compute signed distance from a ref grid (zeroset).

The reference grid (typically travel time, but may be rasterized line offsets) provides a zero contour from which the signed distance is computed by fast marching.

Parameters

f: NDArray Reference grid to provide a zero contour. band_width: float | None, Width of the narrow band. resolution: float Grid pixel width (and height). order: int Fast-marching order of computational stencil (1 or 2).

Returns

MaskedArray Signed distance grid.

Source code in erosionfront/levelset/grid.py

@staticmethod
def fm_distance(
        f: NDArray, 
        band_width: float | None,
        resolution: float,
        order: int | None,
    ) -> MaskedArray:
    """
    Use fast marching to compute signed distance from a ref grid (zeroset).

    The reference grid (typically travel time, but may be rasterized line
    offsets) provides a zero contour from which the signed distance is 
    computed by fast marching.

    Parameters
    ----------
    f: NDArray
        Reference grid to provide a zero contour.
    band_width: float | None,
        Width of the narrow band. 
    resolution: float
        Grid pixel width (and height).
    order: int
        Fast-marching order of computational stencil (1 or 2).

    Returns
    -------
    MaskedArray
        Signed distance grid.
    """
    φ: NDArray | MaskedArray \
          = skfmm.distance(
                f, 
                dx=resolution, 
                order=order, 
                narrow=band_width,
            )
    return (
        φ if type(φ)==MaskedArray
        else np.ma.masked_array(φ, np.zeros_like(φ, dtype=np.bool,))
    )

fm_distance_everywhere

fm_distance_everywhere(
    φ: MaskedArray, φ_everywhere: NDArray, resolution: float, order: int | None
) -> NDArray

FM compute a narrowband signed distance and then extend everywhere.

The previous 'everywhere' grid of signed distance is used to supply values beyond the narrowband.

Parameters

φ: MaskedArray Reference grid to provide a zero contour. φ_everywhere: NDArray Previous grid of signed distance everywhere. resolution: float Grid pixel width (and height). order: int Fast-marching order of computational stencil (1 or [2]).

Attributes

Uses the fm_distance method.

Returns

MaskedArray Grid with signed distance values everywhere, not just on the nb.

Source code in erosionfront/levelset/grid.py

def fm_distance_everywhere(
        self, 
        φ: MaskedArray, 
        φ_everywhere: NDArray,
        resolution: float,
        order: int | None,
    ) -> NDArray:
    """
    FM compute a narrowband signed distance and then extend everywhere.

    The previous 'everywhere' grid of signed distance is used to supply
    values beyond the narrowband.

    Parameters
    ----------
    φ: MaskedArray
        Reference grid to provide a zero contour.
    φ_everywhere: NDArray
        Previous grid of signed distance everywhere.
    resolution: float
        Grid pixel width (and height).
    order: int
        Fast-marching order of computational stencil (1 or [2]).

    Attributes
    ----------
    Uses the fm_distance method.

    Returns
    -------
    MaskedArray
        Grid with signed distance values everywhere, not just on the nb.
    """
    previous_φ_everywhere: NDArray \
        = φ_everywhere
    previous_distant_pixels: NDArray = φ.mask
    φ_everywhere_: NDArray \
        = self.fm_distance(φ, 0, resolution, order,).data
    φ_everywhere_[previous_distant_pixels] \
        = previous_φ_everywhere[previous_distant_pixels]
    return φ_everywhere_

get_extent

get_extent() -> NDArray

Get domain (x,z) extent.

Returns

NDArray: Domain bounds as [x_min, x_max, z_max, z_min].

Source code in erosionfront/levelset/elementary.py

def get_extent(self) -> NDArray:
    """
    Get domain (x,z) extent.

    Returns
    -------
    NDArray:
        Domain bounds as [x_min, x_max, z_max, z_min].

    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return self.domain.extent

get_front_points

get_front_points(f: NDArray | MaskedArray) -> NDArray

Fetch points along erosion front surface aka zeroset of f.

Parameters

f: NDArray | MaskedArray Presumably a signed distance grid (but not necessarily).

Attributes

Uses compute_contour method.

Returns

NDArray Coordinate sequence (x,z)*N along surface as (2,N) array.

Source code in erosionfront/levelset/grid.py

def get_front_points(
        self, 
        f: NDArray | MaskedArray,
    ) -> NDArray:
    """
    Fetch points along erosion front surface aka zeroset of f.

    Parameters
    ----------
    f: NDArray | MaskedArray
        Presumably a signed distance grid (but not necessarily).

    Attributes
    ----------
    Uses compute_contour method.

    Returns
    -------
    NDArray
        Coordinate sequence (x,z)*N along surface as (2,N) array.
    """
    return self.compute_contour(f, level=0,)

get_line_points

get_line_points() -> NDArray

Convert Shapely geometry Line into (x,z) points.

Attributes

Uses line attribute.

Returns

NDArray Line as N coordinates in (2,N) array format, i.e., as two rows of N-length sequences.

Source code in erosionfront/levelset/grid.py

def get_line_points(self) -> NDArray:
    """
    Convert Shapely geometry Line into (x,z) points.

    Attributes
    ----------
    Uses line attribute.

    Returns
    -------
    NDArray
        Line as N coordinates in (2,N) array format, 
        i.e., as two rows of N-length sequences.
    """
    return (shapely.get_coordinates(self.line)).T

get_points_xz

get_points_xz() -> NDArray

Get surface (x,z) coordinates.

Returns

NDArray Points along surface as (x,z) pairs.

Source code in erosionfront/levelset/elementary.py

def get_points_xz(self) -> NDArray:
    """
    Get surface (x,z) coordinates.

    Returns
    -------
    NDArray
        Points along surface as (x,z) pairs.
    """
    return self.points_xz.T

get_substrate_ξ0

get_substrate_ξ0(φ: MaskedArray) -> MaskedArray

Generate relative η grid (effective speed ξ0).

Parameters

φ: MaskedArray Any grid of the size being used in the sim.

Attributes

substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate Instance of class containing substrate info. surface.height: float Top height of surface above z=0 base. η_upperlayer: float Relative erodiblity of strong layer (weak layer assumed = 1). do_hardtop: bool Whether to fake a highly resistant, very thin, upper "caprock". do_hardbase: bool Whether to fake a highly resistant, very thin, lower "bedrock". η_hardtopbase: float Relative erodiblity of base/top layers (weak layer assumed = 1).

Returns

MaskedArray Modulated erodbility (ref speed) grid ξ0 across the narrow band.

Source code in erosionfront/levelset/speed.py

def get_substrate_ξ0(
        self, 
        φ: MaskedArray,
    ) -> MaskedArray:
    """
    Generate relative η grid (effective speed ξ0).

    Parameters
    ----------
    φ: MaskedArray
        Any grid of the size being used in the sim.

    Attributes
    ----------
    substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate
        Instance of class containing substrate info.
    surface.height: float
        Top height of surface above z=0 base.
    η_upperlayer: float
        Relative erodiblity of strong layer (weak layer assumed = 1).
    do_hardtop: bool
        Whether to fake a highly resistant, very thin, upper "caprock".
    do_hardbase: bool
        Whether to fake a highly resistant, very thin, lower "bedrock".
    η_hardtopbase: float
        Relative erodiblity of base/top layers (weak layer assumed = 1).

    Returns
    -------
    MaskedArray
        Modulated erodbility (ref speed) grid ξ0 across the narrow band.
    """
    substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate
    if self.substrate is not None:
        substrate = self.substrate 
    else:
        raise ValueError("Substrate not specified")
    surface: (
        StraightLineSurface | 
        ErrorFunctionSurface | 
        SemiCircularArcSurface | 
        QuarterCircularArcSurface
    )
    if self.surface is not None:
        surface = self.surface 
    else:
        raise ValueError("Surface not specified")
    z_max: float = surface.height
    z_change: Callable = lambda z: int(self.xz_to_ij((0, z,))[1])
    ξ0: MaskedArray = np.ones_like(φ)
    # Prevent erosion of top surface
    if type(substrate)==TwoLayerSubstrate:
        ξ0[z_change(0.5):,:] = substrate.η_upperlayer
    if substrate.do_hardtop:
        # Make the hard top either 5% of height, or 10 pixels high,
        #    - whichever is smaller
        ξ0[z_change(z_max-min(self.Δx*substrate.nz_hardtop, 
                              0.001*substrate.nz_hardtop)):,:] \
            = substrate.η_hardtopbase
    # Prevent erosion of bottom surface and thus motion to left
    if substrate.do_hardbase:
        #BUG: this doesn't work if the initial surface isn't z \in [0,1]
        ξ0[:(z_change(0)),:] = substrate.η_hardtopbase
    # Smear out a bit to stabilize
    ξ0_ = scipy.ndimage.uniform_filter(
        ξ0, 
        size=substrate.nz_hardtop,
    )
    ξ0 = np.ma.masked_array(ξ0_, mask=ξ0.mask,)
    return ξ0

get_x_limits

get_x_limits() -> tuple

Get domain x min, max.

Returns

tuple: x_min, x_max pair.

Source code in erosionfront/levelset/elementary.py

def get_x_limits(self) -> tuple:
    """
    Get domain x min, max.

    Returns
    -------
    tuple:
        x_min, x_max pair.
    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return tuple(self.domain.extent[:2])

get_z_limits

get_z_limits() -> tuple

Get domain z min, max.

Returns

tuple: z_min, z_max pair.

Source code in erosionfront/levelset/elementary.py

def get_z_limits(self) -> tuple:
    """
    Get domain z min, max.

    Returns
    -------
    tuple:
        z_min, z_max pair.
    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return tuple(self.domain.extent[2:][::-1])

get_φ_next

get_φ_next() -> MaskedArray

Return next signed-distance grid φ_next or just φ if not available.

Returns

MaskedArray: Next signed-distance φ grid or just current φ if next unavailable.

Source code in erosionfront/levelset/speed.py

def get_φ_next(self) -> MaskedArray:
    """
    Return next signed-distance grid φ_next or just φ if not available.

    Returns
    -------
    MaskedArray:
        Next signed-distance φ grid or just current φ if next unavailable.
    """
    return (
        self.φ if self.φ_next is None
        else self.φ_next
    )

ij_to_xz

ij_to_xz(ij: NDArray | tuple) -> NDArray

Convert pixel indices [i,j] into pixel center coords (x,z).

Parameters

ij: NDArray | tuple Pixel index pairs [i,j] as numpy array or tuple.

Returns

NDArray: Pixel (x,z) coordinate pairs as float64 numpy array.

Source code in erosionfront/levelset/elementary.py

def ij_to_xz(self, ij: NDArray | tuple,) -> NDArray:
    """
    Convert pixel indices [i,j] into pixel *center* coords (x,z).

    Parameters
    ----------
    ij: NDArray | tuple
        Pixel index pairs [i,j] as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel (x,z) coordinate pairs as float64 numpy array.
    """
    if self.domain is None or (origin:=self.domain.origin()) is None:
        raise ValueError("Domain/origin not specified")
    return (
        self.Δx * np.array(ij).astype(np.float64) 
        + origin 
        + np.array([1,1])*self.Δx/2
    )

ji_to_xz

ji_to_xz(ji: NDArray | tuple) -> NDArray

Convert reversed pixel indices [j,i] into pixel center coords (x,z).

Parameters

ji: NDArray | tuple Pixel (reversed) index pairs [j,i] as numpy array or tuple.

Returns

NDArray: Pixel (x,z) coordinate pairs as float64 numpy array.

Source code in erosionfront/levelset/elementary.py

def ji_to_xz(self, ji: NDArray | tuple,) -> NDArray:
    """
    Convert reversed pixel indices [j,i] into pixel *center* coords (x,z).

    Parameters
    ----------
    ji: NDArray | tuple
        Pixel (reversed) index pairs [j,i] as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel (x,z) coordinate pairs as float64 numpy array.

    """
    return self.ij_to_xz(np.array(ji)[::-1])

make_line staticmethod

make_line(points_xz: NDArray) -> LineString

Convert surface (x,z) coordinates into a Shapely Line geometry.

Parameters

points_xz: NDArray Stacked vectors of x and z surface point coordinates.

Returns

LineString: Shapely LineString geometry of surface line.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def make_line(points_xz: NDArray,) -> LineString:
    """
    Convert surface (x,z) coordinates into a Shapely Line geometry.

    Parameters
    ----------
    points_xz: NDArray
        Stacked vectors of x and z surface point coordinates.

    Returns
    -------
    LineString:
        Shapely LineString geometry of surface line.
    """
    return LineString(points_xz)

make_polygon staticmethod

make_polygon(
    points_xz: NDArray, resolution: float, n_pixels_xz: tuple | NDArray
) -> Polygon

Create a 'bounding' polygon from surface x,z points.

This polygon is used to decide whether a point lies within the bedrock or outside it.

BUG: Crudely done. Should be much tighter geometry close to domain.

Parameters

points_xz: NDArray Surface coordinates n_pixels_xz: tuple | NDArray Grid dimensions resolution: float Grid pixel size.

Returns

Polygon Shapely geometry of bounding polygon.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def make_polygon(
        points_xz: NDArray,
        resolution: float,
        n_pixels_xz: tuple | NDArray,
    ) -> Polygon:
    """
    Create a 'bounding' polygon from surface x,z points.

    This polygon is used to decide whether a point lies within the bedrock
    or outside it.

    BUG: Crudely done. Should be much tighter geometry close to domain.

    Parameters
    ----------
    points_xz: NDArray
        Surface coordinates
    n_pixels_xz: tuple | NDArray
        Grid dimensions
    resolution: float
        Grid pixel size.

    Returns
    -------
    Polygon
        Shapely geometry of bounding polygon.
    """
    if n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    x_offset: float = n_pixels_xz[0]*resolution*10
    z_offset: float = n_pixels_xz[1]*resolution*10
    left_anchor: NDArray \
        = np.array([points_xz[0][0]-x_offset, 
                    points_xz[0][1]-z_offset])
    right_anchor: NDArray \
        = np.array([points_xz[-1][0]+x_offset, 
                    points_xz[-1][1]-z_offset])
    poly_points: NDArray \
        = np.vstack((
            left_anchor, 
            points_xz[0]-np.array([resolution*10,0]),
            points_xz, 
            points_xz[-1]+np.array([resolution*10,0]),
            right_anchor,
        ))
    return Polygon(poly_points)

measure_normal_distance staticmethod

measure_normal_distance(
    line: tuple | NDArray | LineString, point: tuple | NDArray | Point
) -> float

Compute distance between line and point using Shapely.

Parameters

line: tuple | NDArray | LineString Line represented by various possible types. point: tuple | NDArray | Point Coordinate (x,z) represented by various possible types.

Returns

float Perpendicular distance from line curve to point.

Source code in erosionfront/levelset/grid.py

@staticmethod
def measure_normal_distance(
        line: tuple | NDArray | LineString, 
        point: tuple | NDArray | Point,
    ) -> float:
    """
    Compute distance between line and point using Shapely.

    Parameters
    ----------
    line: tuple | NDArray | LineString
        Line represented by various possible types.
    point: tuple | NDArray | Point
        Coordinate (x,z) represented by various possible types.

    Returns
    -------
    float
        Perpendicular distance from line curve to point.
    """
    return float(shapely.distance(
        line if type(line) is LineString else LineString(line), 
        point if type(point) is Point else Point(point)
    ))

measure_φ_ji

measure_φ_ji(line: LineString, ji: tuple[int, int] | NDArray) -> float

Compute signed distance between line and i,j pixel using Shapely.

The sign is -1 for "in the air" and +1 for "in the rock".

Parameters

line: LineString Surface line. ji: tuple[int,int] | NDArray Pixel indexes in reversed format.

Attributes

Uses measure_φ_xz method to test which side of the surface the pixel lies.

Returns

float Signed perpendicular distance from the line to the pixel.

Source code in erosionfront/levelset/grid.py

def measure_φ_ji(
        self,
        line: LineString,
        ji: tuple[int,int] | NDArray,
    ) -> float:
    """
    Compute signed distance between line and i,j pixel using Shapely.

    The sign is -1 for "in the air" and +1 for "in the rock".

    Parameters
    ----------
    line: LineString
        Surface line.
    ji: tuple[int,int] | NDArray
        Pixel indexes in reversed format.

    Attributes
    ----------
    Uses measure_φ_xz method to test which side of the 
    surface the pixel lies.

    Returns
    -------
    float
        Signed perpendicular distance from the line to the pixel.
    """
    return self.measure_φ_xz(line, self.ji_to_xz(ji),)

measure_φ_xz

measure_φ_xz(line: LineString, xz: tuple[float, float] | NDArray) -> float

Compute signed distance between line and xz point using Shapely.

The sign is -1 for "in the air" and +1 for "in the rock".

Parameters

line: LineString Surface line. xz: tuple[float,float] | NDArray Point coordinates.

Attributes

Uses above_or_below_xz method to test which side of the surface the point lies.

Returns

float Signed perpendicular distance from the line to the point.

Source code in erosionfront/levelset/grid.py

def measure_φ_xz(
        self,
        line: LineString,
        xz: tuple[float,float] | NDArray,
    ) -> float:
    """
    Compute signed distance between line and xz point using Shapely.

    The sign is -1 for "in the air" and +1 for "in the rock".

    Parameters
    ----------
    line: LineString
        Surface line.
    xz: tuple[float,float] | NDArray
        Point coordinates.

    Attributes
    ----------
    Uses above_or_below_xz method to test which side of the surface the 
    point lies.

    Returns
    -------
    float
        Signed perpendicular distance from the line to the point.
    """
    return (
        self.measure_normal_distance(line, xz)
        *self.above_or_below_xz(xz, self.in_rock_polygon,)
    )

pad_domain staticmethod

pad_domain(domain: Domain, resolution: float, n_pad_pixels: int) -> Domain

Expand the domain vertically by pixel width amounts.

Parameters

domain: Domain (Initial) domain to be padded. resolution: float Grid pixel size. n_pad_pixels: int Number of pixel widths to pad by below the current domain.

Returns

New domain instance padded vertically.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def pad_domain(
        domain: Domain, 
        resolution: float,
        n_pad_pixels: int,
    ) -> Domain:
    """
    Expand the domain vertically by pixel width amounts.

    Parameters
    ----------
    domain: Domain
        (Initial) domain to be padded.
    resolution: float
        Grid pixel size.
    n_pad_pixels: int
        Number of pixel widths to pad by *below* the current domain.

    Returns
    -------
    New domain instance padded vertically.
    """
    dz_bot: float = (n_pad_pixels)*resolution
    dz_top: float = (n_pad_pixels)*resolution
    # dz_top: float = resolution*5 if n_pad_pixels>0 else 0
    n_digits: int = int(round(-np.log10(resolution)+0.5))+1
    return Domain(
        round(domain.x_min, n_digits),
        round(domain.x_max, n_digits),
        round(domain.z_min - dz_bot, n_digits),
        round(domain.z_max + dz_top, n_digits),
    )

pad_grid staticmethod

pad_grid(grid: NDArray, n_pad_pixels: int) -> NDArray

Pad grid with specified number of pixels below and above.

Parameters

grid: NDArray Grid to be padded. n_pad_pixels: int Number of pixels to pad below.

Returns

NDArray Padded grid.

Source code in erosionfront/levelset/elementary.py

@staticmethod
def pad_grid(
        grid: NDArray, 
        n_pad_pixels: int, 
    ) -> NDArray:
    """
    Pad grid with specified number of pixels below and above.

    Parameters
    ----------
    grid: NDArray
        Grid to be padded.
    n_pad_pixels: int
        Number of pixels to pad below.

    Returns
    -------
    NDArray
        Padded grid.
    """
    pad_ztop: NDArray \
        = np.zeros((n_pad_pixels,grid.shape[1],), dtype=np.int64,)
    pad_zbot: NDArray \
        = np.zeros((n_pad_pixels,grid.shape[1],), dtype=np.int64,)
    expanded_grid: NDArray 
    if n_pad_pixels>0:
        expanded_grid =  np.vstack((pad_zbot, grid, pad_ztop,))
    else:
        expanded_grid = grid
    return expanded_grid

rasterize_line

rasterize_line(points_xz: NDArray, do_extend_line: bool = True) -> NDArray

Rasterize surface line as x,z points into grid pixels using Rasterio.

The x,z bounds of the raster grid are precisely set by the bounding box of the line.

Parameters

points_xz: NDArray Surface line as a series of x,z points. do_extend_line: bool=True Artificially extend the surface line so it reaches the boundary.

Attributes

Uses n_pixels_xz and xz_to_ij methods.

Returns

NDArray Grid of rasterized line.

Source code in erosionfront/levelset/grid.py

def rasterize_line(
        self, 
        points_xz: NDArray,
        do_extend_line: bool=True,
    ) -> NDArray:
    """
    Rasterize surface line as x,z points into grid pixels using Rasterio.

    The x,z bounds of the raster grid are precisely set by the
    bounding box of the line.

    Parameters
    ----------
    points_xz: NDArray
        Surface line as a series of x,z points.
    do_extend_line: bool=True
        Artificially extend the surface line so it reaches the boundary.

    Attributes
    ----------
    Uses n_pixels_xz and xz_to_ij methods.

    Returns
    -------
    NDArray
        Grid of rasterized line.
    """
    #HACK: artificially extend line horizontally to ensure
    #      we get correct bounds
    points_xz_extended: NDArray
    if do_extend_line:
        points_xz_extended  = np.vstack((
            [2*points_xz[0][0],0],
            points_xz,
            [2*points_xz[-1][0],1],
        ))
    else:
        points_xz_extended = points_xz.copy()
    pts_ij = self.xz_to_ij(points_xz_extended)
    line: LineString = LineString(pts_ij)
    if self.n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    # Specific caching actually slows rasterization down!
    #   with rasterio.Env(GDAL_CACHEMAX=4_096_000_000) as env:
    ρ: NDArray \
        = rasterio.features.rasterize(
            [line],
            out_shape=self.n_pixels_xz[::-1], 
            all_touched=True, 
            fill=0,
        )
    # assert ρ.shape==self.n_pixels_xz[::-1]
    return ρ #, line

redo_φ_everywhere staticmethod

redo_φ_everywhere(φ_everywhere: NDArray, φ: MaskedArray) -> NDArray

Reassign signed distance values to pixels beyond the narrow band.

Assumes that pixels beyond the narrow band have signed-distance signs that are currently correct: given these signs, each beyond-nb pixel is assigned either the min or max value of the nb signed-distance values.

Parameters

φ_everywhere: NDArray Beyond narrow-band signed-distance grid to be recomputed. φ: MaskedArray Narrow-band signed-distance grid.

Returns

NDArray Recomputed signed-distances beyond narrow band.

Source code in erosionfront/levelset/grid.py

@staticmethod
def redo_φ_everywhere(
        φ_everywhere: NDArray, 
        φ: MaskedArray,
    ) -> NDArray:
    """
    Reassign signed distance values to pixels beyond the narrow band.

    Assumes that pixels beyond the narrow band have signed-distance signs
    that are currently correct: given these signs, each beyond-nb pixel is 
    assigned either the min or max value of the nb signed-distance values.

    Parameters
    ----------
    φ_everywhere: NDArray
        Beyond narrow-band signed-distance grid to be recomputed.
    φ: MaskedArray
        Narrow-band signed-distance grid.

    Returns
    -------
    NDArray
        Recomputed signed-distances beyond narrow band.
    """
    mask: NDArray = np.array(φ.mask)
    far_recomputed: NDArray = φ_everywhere.copy()
    far_recomputed[φ_everywhere<0] = np.min(φ)
    far_recomputed[φ_everywhere>0] = np.max(φ)
    far_recomputed[~mask] = φ[~mask]
    return far_recomputed

xz_to_ij

xz_to_ij(xz: NDArray | tuple) -> NDArray

Convert pixel center coords (x,z) into pixel indices [i,j].

Parameters

xz: NDArray | tuple Pixel (x,z) coordinate pairs as numpy array or tuple.

Returns

NDArray: Pixel index pairs [i,j] as int64 numpy array

Source code in erosionfront/levelset/elementary.py

def xz_to_ij(self, xz: NDArray | tuple,) -> NDArray:
    """
    Convert pixel *center* coords (x,z) into pixel indices [i,j].

    Parameters
    ----------
    xz: NDArray | tuple
        Pixel (x,z) coordinate pairs as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel index pairs [i,j] as int64 numpy array
    """
    if self.domain is None or (origin:=self.domain.origin()) is None:
        raise ValueError("Domain/origin not specified")
    return np.int64(  #type: ignore
        np.array(xz)/self.Δx 
        - origin/self.Δx
    )