Skip to content

base.py

LevelSetSolverBase dataclass 

LevelSetSolverBase(
    Δ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.base.LevelSetSolverBase[LevelSetSolverBase]
              erosionfront.levelset.speed.LevelSetSpeedMethods[LevelSetSpeedMethods]
              erosionfront.levelset.grid.LevelSetGridMethods[LevelSetGridMethods]
              erosionfront.levelset.elementary.LevelSetElementaryMethods[LevelSetElementaryMethods]
              erosionfront.levelset.data.LevelSetData[LevelSetData]

                              erosionfront.levelset.speed.LevelSetSpeedMethods --> erosionfront.levelset.base.LevelSetSolverBase
                                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.base.LevelSetSolverBase href "" "erosionfront.levelset.base.LevelSetSolverBase"
              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"

Solver for propagating an erosion front using a level-set scheme.

Subclasses LevelSetSpeedMethods.

Attributes

ρ: NDArray Rasterized surface-line grid. ρ_padded: NDArray Rasterized surface-line grid with padding. ρ_offsets: NDArray Gridded surface-normal distances. T_slices: dict Dictionary recording time T slice data of erosion front evolution. T_arrival: MaskedArray | None = None Grid of interpolated arrival times T, with mask removing pixels in air or in uneroded substrate.

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).

  • evolve –

    Iteratively propagate the erosion front surface.

  • 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).

  • initialize –

    Initialize level-set solver.

  • 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_T_slices –

    Rasterize time T slices onto grids.

  • rasterize_line –

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

  • record_T_slice –

    Record selected sim data for a single time step.

  • redo_φ_everywhere –

    Reassign signed distance values to pixels beyond the narrow band.

  • update –

    Update method that must be implemented by the actual solver.

  • 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 
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
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 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
@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 
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
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 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
@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 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
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 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
@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 
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
@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,)

evolve 

evolve(
    Δt: float = 0,
    t_span: float | None = None,
    n_steps: int | None = None,
    Δt_reinitialize: float | None = None,
    κ_boost: float = 1,
    report_pc: int = 5,
    n_record_steps: int = 5,
    do_initial: bool = False,
    do_timing: bool = False,
) -> None

Iteratively propagate the erosion front surface.

Parameters

Δt: float=0 Simulation time step. t_span: float | None = None Duration of simulation. n_steps: int | None = None Number of time steps to take in the simulation. If neither t_span nor n_steps is specified, If neither are, a single step is assumed. Δt_reinitialize: float | None = None Period of reinitialization of φ signed-distance function grid. κ_boost: float = 1 Scale factor for boosting numerical viscosity. report_pc: int = 5 Reporting interval as percentage of runtime. n_record_steps: int = 5 Number of steps between T slicings. do_initial: bool=False Specify if this is the first update, and act accordingly. do_timing: bool=False Choose whether to perform timings of sim stages. Attributes

Most of the sim grids are modified.

Source code in erosionfront/levelset/base.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def evolve(
        self, 
        Δt: float=0, 
        t_span: float | None = None,
        n_steps: int | None = None,
        Δt_reinitialize: float | None = None, 
        κ_boost: float = 1,
        report_pc: int = 5,
        n_record_steps: int = 5,
        do_initial: bool=False,
        do_timing: bool=False,
    ) -> None:
    """
    Iteratively propagate the erosion front surface.

    Parameters
    ----------
    Δt: float=0
        Simulation time step.
    t_span: float | None = None
        Duration of simulation.
    n_steps: int | None = None
        Number of time steps to take in the simulation.
        If neither `t_span` nor `n_steps` is specified, 
        If neither are, a single step is assumed.
    Δt_reinitialize: float | None = None
        Period of reinitialization of φ signed-distance function grid.
    κ_boost: float = 1
        Scale factor for boosting numerical viscosity.
    report_pc: int = 5
        Reporting interval as percentage of runtime.
    n_record_steps: int = 5
        Number of steps between T slicings.
    do_initial: bool=False
        Specify if this is the first update, and act accordingly.
    do_timing: bool=False
        Choose whether to perform timings of sim stages.
    Attributes
    ----------
    Most of the sim grids are modified.
    """
    self.Δt = Δt
    if self.Δt>0 and self.model is not None:
        self.n_slabfailure = int(self.model.Δt_slabfailure/self.Δt)
    self.κ_boost = κ_boost
    if Δt_reinitialize is None:
        self.Δt_reinitialize = Δt
    else:
        self.Δt_reinitialize = Δt_reinitialize
    n_steps_: int
    t_span_: float
    if n_steps is None:
        if t_span is None:
            n_steps_ = 1
            t_span_ = Δt
        else:
            n_steps_ = int(t_span/Δt+0.5)
            t_span_ = t_span
    else:
        n_steps_ = n_steps
        t_span_ = n_steps*Δt
    print(
        rf"Run time T={t_span_} for n_steps={n_steps_} "
        + rf"with Δt={Δt} and Δx={self.Δx}"
    )
    do_report: Callable \
        = lambda step: (
            False if int(n_steps_*report_pc/100)==0
            else (step+1)%int(n_steps_*report_pc/100)==0
        )
    # do_record: Callable \
    #     = lambda step: (
    #         False if int(n_steps_*record_pc/100)==0 
    #         else (step+1)%int(n_steps_*record_pc/100)==0
    #     )
    rnd: int = int(np.round(-np.log10(Δt if Δt>0 else 1))+1)
    def record_report(step) -> None:
        progress_pc_: float
        if do_report(step):
            progress_pc_ = int(100*(step+1)/(n_steps_))
            print(f"   {progress_pc_}%:  step={self.i_step}"
                + f"  t={np.round(self.t_total, rnd)}")
        if ((step+1) % n_record_steps)==0:
            self.record_T_slice(
                self.i_step if self.i_step is not None else 1, 
                self.t_total, 
            )

    #######################################################################
    if self.t_total is None or do_initial:
        # Initial do-nothing step...
        self.record_T_slice(0, 0,)
        self.i_step = 0
        self.update(0, do_timing=do_timing,)
    else:
        # ...or a full run
        for step_ in range(n_steps_):
            self.i_step = self.i_step+1 if self.i_step is not None else 1
            self.update(Δt, Δt_reinitialize,)
            record_report(step_)
        self.update(Δt, do_finalize=True, do_timing=do_timing,)
    self.n_total = self.i_step

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 
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
@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 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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 
120
121
122
123
124
125
126
127
128
129
130
131
132
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 
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
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 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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 
134
135
136
137
138
139
140
141
142
143
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 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
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 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
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 
107
108
109
110
111
112
113
114
115
116
117
118
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
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 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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
    )

initialize 

initialize(do_timing: bool = False) -> None

Initialize level-set solver.

Parameters

do_timing: bool=False Perform timing estimates for each computational stage.

Attributes

Many class attributes are modified.

Source code in erosionfront/levelset/base.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def initialize(self, do_timing: bool=False,) -> None:
    """
    Initialize level-set solver.

    Parameters
    ----------
    do_timing: bool=False
        Perform timing estimates for each computational stage.

    Attributes
    ----------
    Many class attributes are modified.
    """
    Timer(do_timing, "Initializing",)

    Timer.start("create initial surface line; generate subsurface polygon",)
    self.domain = copy(self.raw_domain)
    if self.surface is None:
        raise ValueError("Surface not specified")
    self.points_xz = self.surface.get_points_xz()
    self.line = self.make_line(self.points_xz)
    if self.n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    self.in_rock_polygon \
        = self.make_polygon(
            self.points_xz, 
            self.Δx,
            self.n_pixels_xz, 
        )
    Timer.stop()

    Timer.start("rasterize line",)
    self.ρ = self.rasterize_line(self.points_xz, self.do_extend_line,)
    Timer.stop()
    Timer.start("pad the grid vertically",)
    if self.raw_domain is None:
        raise ValueError("Domain not specified")
    if self.n_pad_pixels is None:
        raise ValueError("Number of pad pixels not specified")
    self.domain = self.pad_domain(
        self.raw_domain, 
        self.Δx,
        self.n_pad_pixels,
    )
    self.n_pixels_xz \
        = self.compute_pixel_dimensions(self.domain, self.Δx,)
    self.n_pixels_xz0 \
        = self.compute_pixel_origin(self.domain, self.Δx,)
    if self.n_pad_pixels is None:
        raise ValueError("Number of padding pixels not specified")
    self.ρ_padded \
        = self.pad_grid(
            grid=self.ρ, 
            n_pad_pixels=self.n_pad_pixels,
        )
    Timer.stop()

    Timer.start("measure pixel offsets from line",)
    self.ρ_offsets \
        = self.compute_line_pixel_offsets(
            line=self.line,
            ρ=self.ρ_padded,
            do_φ_everywhere=self.do_φ_everywhere,
        )
    Timer.stop()

    Timer.start("fast-march signed distance from pixel offsets",)
    self.φ \
        = self.fm_distance(
            f=self.ρ_offsets, 
            band_width=self.band_width,
            resolution=self.Δx,
            order=self.fm_order,
        )
    Timer.stop()

    Timer.start("fill whole signed-distance grid with very-far values",)
    self.φ_everywhere \
        = self.redo_φ_everywhere(
            φ_everywhere=self.ρ_offsets, 
            φ=self.φ,
        )
    Timer.stop()

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 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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 
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
@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 
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
@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 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
@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 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
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 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
@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 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
@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_T_slices 

rasterize_T_slices(n_subset: int = 1) -> None

Rasterize time T slices onto grids.

Parameters

n_subset: int=1 Rate of subsetting of T slices (n=1 means all of them).

Source code in erosionfront/levelset/base.py 
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
def rasterize_T_slices(self, n_subset: int=1,) -> None:
    """
    Rasterize time T slices onto grids.

    Parameters
    ----------
    n_subset: int=1
        Rate of subsetting of T slices (n=1 means all of them).
    """
    print("Gridding T slices...", end="",)
    if self.domain is None:
        raise ValueError("Domain not defined")
    bounds: NDArray | None= self.domain.bounds()
    if bounds is None:
        raise ValueError("Domain bounds not defined")
    if self.n_pixels_xz is None:
        raise ValueError("Number of x,z pixels not defined")
    n_x: int = self.n_pixels_xz[0]
    n_z: int = self.n_pixels_xz[1]
    x_bounds: NDArray = (bounds.T)[0]
    z_bounds: NDArray = (bounds.T)[1]
    x_pts: NDArray = np.linspace(x_bounds[0], x_bounds[1], n_x,) + self.Δx/2
    z_pts: NDArray = np.linspace(z_bounds[0], z_bounds[1], n_z,) + self.Δx/2
    assert x_pts.shape[0]==n_x
    assert z_pts.shape[0]==n_z
    x_mg: NDArray
    z_mg: NDArray
    (z_mg, x_mg,) = np.meshgrid(z_pts, x_pts, indexing="ij",)
    x_points: NDArray = np.concatenate(
        [slice_["x"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    z_points: NDArray = np.concatenate(
        [slice_["z"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    T_values: NDArray = np.concatenate(
        [slice_["x"]*0+slice_["T"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    zx_points: NDArray = np.vstack((z_points, x_points,)).T
    assert type(zx_points) is np.ndarray
    assert type(T_values) is np.ndarray
    T_arrival: NDArray \
        = griddata(zx_points, T_values, (z_mg, x_mg), method="linear",)
    not_reached: NDArray = np.zeros_like(self.φ_everywhere).astype(np.bool)
    not_reached[self.φ_everywhere<0] = True
    self.T_arrival = np.ma.masked_array(T_arrival, mask=not_reached,)
    print(" done")

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 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
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

record_T_slice 

record_T_slice(step: int, t_total: float) -> None

Record selected sim data for a single time step.

Recorded data are placed in a 'T_slice' dictionary.

Parameters

step: int Simulation step number. t_total: float Total time in simulation so far.

Attributes

Data for this sim time step are appended to the time-slicing dict.

Source code in erosionfront/levelset/base.py 
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
def record_T_slice(
        self, 
        step: int, 
        t_total: float, 
    ) -> None:
    """
    Record selected sim data for a single time step.

    Recorded data are placed in a 'T_slice' dictionary.

    Parameters
    ----------
    step: int
        Simulation step number.
    t_total: float
        Total time in simulation so far.

    Attributes
    ----------
    Data for this sim time step are appended to the time-slicing dict.
    """
    if self.φ_next is None:
        d = self.φ/np.max(self.φ)
    else:
        d = self.φ_next/np.max(self.φ_next)
    # d: NDArray | MaskedArray = self.φ_everywhere/np.max(self.φ_everywhere)
    xz_profile: NDArray = self.get_front_points(d)
    x: NDArray = xz_profile[0]
    z: NDArray = xz_profile[1]
    β: NDArray = np.arctan2(np.gradient(z), np.gradient(x))
    self.T_slices[step] = {
        "x" : x,
        "z" : z,
        "φ": self.φ,
        "β" : β,
        "T" : float(t_total),
    }

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 
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
@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

update abstractmethod 

update(
    Δt: float,
    Δt_reinitialize: float | None = None,
    do_finalize: bool = False,
    do_timing: bool = False,
) -> None

Update method that must be implemented by the actual solver.

Parameters

Δt: float Time step. Δt_reinitialize: float | None=None Periodicity of reinitialization of φ signed-distance grid. do_finalize: bool=False Just finish up, don't do another loop of computation. do_timing: bool=False Measure and report compute times taken at steps during the loop.

Attributes

Most of the sim grids are modified.

Source code in erosionfront/levelset/base.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
@abstractmethod
def update(
        self, 
        Δt: float, 
        Δt_reinitialize: float | None=None,
        do_finalize: bool=False,
        do_timing: bool=False,
    ) -> None:
    """
    Update method that must be implemented by the actual solver.

    Parameters
    ----------
    Δt: float
        Time step.
    Δt_reinitialize: float | None=None
        Periodicity of reinitialization of φ signed-distance grid.
    do_finalize: bool=False
        Just finish up, don't do another loop of computation.
    do_timing: bool=False
        Measure and report compute times taken at steps during the loop.

    Attributes
    ----------
    Most of the sim grids are modified.
    """
    pass

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 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
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
    )