"""Lightweight terrain and bathymetry surfaces for generated farms."""
from __future__ import annotations
import json
import math
from abc import ABC, abstractmethod
from dataclasses import dataclass, field
from pathlib import Path
from typing import TYPE_CHECKING, ClassVar
import numpy as np
from ._random import as_rng, pin_seed, resolve_seed
from ._repr import _ArgumentRepr
from ._validation import (
_validate_finite,
_validate_nonnegative,
_validate_positive,
_validate_real_array,
)
from .boundaries.base import Boundary
if TYPE_CHECKING: # pragma: no cover
import matplotlib.axes
from numpy.typing import ArrayLike
__all__ = [
"Surface",
"Terrain",
"Bathymetry",
"SurfaceModel",
"PerlinTerrain",
"PerlinBathymetry",
"GaussianTerrain",
"GaussianBathymetry",
]
# Surface grids are limited to fifty million nodes and a 4 GiB working set,
# so a runaway domain fails fast rather than swapping, as layouts limit their
# fills. Perlin noise is computed in blocks, and its output, the copy a
# surface keeps and the blocks peak at about 32 bytes a node.
_SURFACE_NODE_BUDGET = 50_000_000
_SURFACE_BYTE_BUDGET = 4 * 1024**3
_PERLIN_OUTPUT_BYTES_PER_NODE = 4 * np.dtype(float).itemsize
#: Nodes of Perlin noise computed at once.
_PERLIN_BLOCK_NODES = 1 << 20
_PERLIN_GRADIENT_BYTES_PER_NODE = 2 * np.dtype(float).itemsize
_GAUSSIAN_OUTPUT_BYTES_PER_NODE = 3 * np.dtype(float).itemsize
_GAUSSIAN_INTERPOLATION_BYTES_PER_NODE = 2 * np.dtype(float).itemsize
def _grid_shape(
bounds: tuple[float, float, float, float],
spacing: float,
*,
extra: int = 1,
) -> tuple[int, int]:
if not np.isfinite(spacing) or spacing <= 0.0:
raise ValueError("surface grid dimensions exceed the allocation budget")
min_x, min_y, max_x, max_y = bounds
width_x, width_y = max_x - min_x, max_y - min_y
ratio_x, ratio_y = width_x / spacing, width_y / spacing
if not all(np.isfinite(ratio) and ratio >= 0.0 for ratio in (ratio_x, ratio_y)):
raise ValueError("surface grid dimensions exceed the allocation budget")
return (
max(2, math.ceil(ratio_x) + extra),
max(2, math.ceil(ratio_y) + extra),
)
def _check_grid_budget(
label: str,
shape: tuple[int, int],
*,
bytes_per_node: int,
parameter: tuple[str, float],
) -> None:
"""Raise ``ValueError`` when a grid of ``shape`` exceeds the budget.
``parameter`` is the name and value of the length that sets the grid's
spacing, which the message tells the caller to raise.
"""
nodes = math.prod(shape)
required_bytes = nodes * bytes_per_node
if nodes > _SURFACE_NODE_BUDGET or required_bytes > _SURFACE_BYTE_BUDGET:
name, value = parameter
raise ValueError(
f"{label} requires {nodes:,} nodes and approximately "
f"{required_bytes:,} bytes, exceeding the surface allocation budget "
f"of {_SURFACE_NODE_BUDGET:,} nodes or {_SURFACE_BYTE_BUDGET:,} "
f"bytes at {name} {value:g} m. Raise {name} or shrink the domain"
)
def _grid_axes(
bounds: tuple[float, float, float, float], shape: tuple[int, int]
) -> tuple[np.ndarray, np.ndarray]:
min_x, min_y, max_x, max_y = bounds
nx, ny = shape
return np.linspace(min_x, max_x, nx), np.linspace(min_y, max_y, ny)
def _validate_range(name: str, values) -> tuple[float, float]:
values = _validate_real_array(name, values)
if values.shape != (2,) or not np.all(np.isfinite(values)):
raise ValueError(f"{name} must contain two finite values")
if not values[0] < values[1]:
raise ValueError(f"{name} must be strictly increasing")
return float(values[0]), float(values[1])
def _validate_grid(x, y, values) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
x = _validate_real_array("surface x", x)
y = _validate_real_array("surface y", y)
values = _validate_real_array("surface values", values)
if x.ndim != 1 or y.ndim != 1 or len(x) < 2 or len(y) < 2:
raise ValueError(
"surface x and y axes must be one-dimensional with at least 2 points"
)
if not np.all(np.isfinite(x)) or not np.all(np.isfinite(y)):
raise ValueError("surface coordinates must be finite")
if not np.all(np.diff(x) > 0) or not np.all(np.diff(y) > 0):
raise ValueError("surface coordinates must be strictly increasing")
if values.shape != (len(y), len(x)):
raise ValueError("surface values must have shape (len(y), len(x))")
if not np.all(np.isfinite(values)):
raise ValueError("surface values must be finite")
return x, y, values
def _cell_index(axis: np.ndarray, values: np.ndarray) -> np.ndarray:
"""Lower node index of the grid cell holding each value in ``axis`` range.
Equals ``searchsorted(axis, values, side="right") - 1`` clipped to
``len(axis) - 2``. The index is guessed from even spacing, corrected by one
node and checked, and values a guess misses use ``searchsorted``.
"""
last = len(axis) - 1
scale = last / (axis[-1] - axis[0])
index = ((values - axis[0]) * scale).astype(np.intp)
np.clip(index, 0, last, out=index)
index -= axis[index] > values
upper = np.minimum(index + 1, last)
index += (axis[upper] <= values) & (upper > index)
upper = np.minimum(index + 1, last)
missed = (axis[index] > values) | ((upper > index) & (axis[upper] <= values))
if missed.any():
index[missed] = np.searchsorted(axis, values[missed], side="right") - 1
return np.minimum(index, last - 1, out=index)
[docs]
@dataclass(init=False)
class Surface:
"""A regular grid of physical values in meters.
``x`` and ``y`` are strictly increasing node coordinates and ``values`` has
shape ``(len(y), len(x))``. Invalid or nonfinite grids raise ``ValueError``.
"""
x: np.ndarray
y: np.ndarray
values: np.ndarray
metadata: dict = field(default_factory=dict)
kind: ClassVar[str] = "surface"
def __init__(
self,
x: ArrayLike,
y: ArrayLike,
values: ArrayLike,
*,
metadata: dict | None = None,
) -> None:
self.x, self.y, self.values = _validate_grid(x, y, values)
self.metadata = {} if metadata is None else metadata
if not isinstance(self.metadata, dict):
raise TypeError("surface metadata must be a dictionary")
def __repr__(self) -> str:
return (
f"{type(self).__name__}(shape={self.values.shape}, "
f"minimum={self.minimum:.6g}, maximum={self.maximum:.6g})"
)
def __eq__(self, other: object) -> bool:
"""Whether ``other`` is a surface of the same kind, grid and metadata."""
if other.__class__ is not self.__class__:
return NotImplemented
return (
np.array_equal(self.x, other.x)
and np.array_equal(self.y, other.y)
and np.array_equal(self.values, other.values)
and self.metadata == other.metadata
)
@property
def minimum(self) -> float:
"""Smallest surface value."""
return float(np.min(self.values))
@property
def maximum(self) -> float:
"""Largest surface value."""
return float(np.max(self.values))
[docs]
def at(self, points: ArrayLike) -> np.ndarray:
"""Bilinearly interpolate values at ``(N, 2)`` points, shape ``(N,)``.
Raises ``TypeError`` for non-numeric points and ``ValueError`` for
points outside the grid.
"""
points = _validate_real_array("points", points)
if points.ndim != 2 or points.shape[1] != 2:
raise ValueError("points must have shape (N, 2)")
if not np.all(np.isfinite(points)):
raise ValueError("points must contain only finite values")
if (
np.any(points[:, 0] < self.x[0])
or np.any(points[:, 0] > self.x[-1])
or np.any(points[:, 1] < self.y[0])
or np.any(points[:, 1] > self.y[-1])
):
raise ValueError("points lie outside the surface grid")
ix = _cell_index(self.x, points[:, 0])
iy = _cell_index(self.y, points[:, 1])
tx = (points[:, 0] - self.x[ix]) / (self.x[ix + 1] - self.x[ix])
ty = (points[:, 1] - self.y[iy]) / (self.y[iy + 1] - self.y[iy])
v00 = self.values[iy, ix]
v10 = self.values[iy, ix + 1]
v01 = self.values[iy + 1, ix]
v11 = self.values[iy + 1, ix + 1]
return (
(1.0 - tx) * (1.0 - ty) * v00
+ tx * (1.0 - ty) * v10
+ (1.0 - tx) * ty * v01
+ tx * ty * v11
)
[docs]
def plot(
self,
ax: "matplotlib.axes.Axes | None" = None,
*,
cmap: str | None = None,
title: str | None = None,
colorbar: bool = True,
) -> "matplotlib.axes.Axes":
"""Plot the surface and return its Matplotlib axes.
``cmap`` defaults to the style's ``image.cmap``.
"""
from .plotting import _new_axes, _style_axes
if ax is None:
ax = _new_axes()
# Gouraud shading treats values as nodes, so the mesh spans the grid exactly.
mesh = ax.pcolormesh(self.x, self.y, self.values, shading="gouraud", cmap=cmap)
if colorbar:
ax.figure.colorbar(
mesh,
ax=ax,
label="elevation [m]" if self.kind == "terrain" else "depth [m]",
)
_style_axes(ax, equal=True, grid=False)
ax.set(xlabel="x [m]", ylabel="y [m]")
if title:
ax.set_title(title)
return ax
[docs]
def save(self, path: "str | Path") -> "Path":
"""Write the surface to ``path``, a ``.nc`` file, and return the path.
The file is the NetCDF grid :func:`oowl.save` writes beside a site
with ``surface_format="netcdf"``, ``x``, ``y`` and ``depth`` or
``elevation`` in meters, which xarray and windIO read. The metadata
is JSON text in the global attribute ``oowl_metadata``, so
:meth:`load` rebuilds an equal surface. Raises ``ValueError`` for
another suffix or metadata that JSON cannot hold.
"""
from ._windio import _write_surface_netcdf
path = Path(path)
if path.suffix.lower() != ".nc":
raise ValueError(f"a surface is saved as a .nc file, not {path.name}")
try:
text = json.dumps(self.metadata, allow_nan=False)
except (TypeError, ValueError) as error:
raise ValueError(
"surface metadata must be finite and JSON-serializable to save"
) from error
_write_surface_netcdf(self, path, None, metadata=text)
return path
[docs]
@classmethod
def load(cls, path: "str | Path") -> "Surface":
"""Read the surface :meth:`save` wrote to ``path``.
A file holding ``depth`` is a :class:`Bathymetry` and one holding
``elevation`` a :class:`Terrain`, so ``Terrain.load`` of a depth grid
raises ``ValueError``. A grid saved by another tool loads with the
metadata ``{"kind": ..., "source": "netcdf"}``.
"""
from ._windio import _read_surface_netcdf_file, _surface_from_windio
path = Path(path)
entry, text = _read_surface_netcdf_file(path)
kind = "bathymetry" if "depth" in entry else "terrain"
surface = _surface_from_windio(entry, kind)
if not isinstance(surface, cls):
raise ValueError(
f"{path.name} holds a {type(surface).__name__}, not a {cls.__name__}"
)
surface.metadata = (
{"kind": kind, "source": "netcdf"} if text is None else json.loads(text)
)
return surface
[docs]
@dataclass(init=False, repr=False, eq=False)
class Terrain(Surface):
"""A surface of elevations relative to a common datum."""
kind: ClassVar[str] = "terrain"
[docs]
@dataclass(init=False, repr=False, eq=False)
class Bathymetry(Surface):
"""A surface of positive water depths, for offshore sites."""
kind: ClassVar[str] = "bathymetry"
def __init__(
self,
x: ArrayLike,
y: ArrayLike,
values: ArrayLike,
*,
metadata: dict | None = None,
) -> None:
super().__init__(x, y, values, metadata=metadata)
if np.any(self.values <= 0):
raise ValueError("bathymetry depths must be positive")
[docs]
class SurfaceModel(_ArgumentRepr, ABC):
"""Surface generator realized against a farm boundary.
A ``seed`` given as a ``numpy.random.Generator`` is pinned to one integer
at construction, so every farm drawn from the model replays alike. Grids
above 50,000,000 nodes or an estimated 4 GiB raise ``ValueError`` before
allocation, naming the length to raise.
"""
[docs]
@abstractmethod
def generate(
self,
boundary: "Boundary",
seed: int | np.random.Generator | None = None,
) -> Surface:
"""Generate a surface covering the boundary bounds.
``seed`` overrides the model's seed.
"""
class _PerlinSurface(SurfaceModel):
_surface_type: ClassVar[type[Surface]]
_range_name: ClassVar[str]
def __init__(
self,
value_range: tuple[float, float],
*,
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
self.value_range = _validate_range(self._range_name, value_range)
self.correlation_length = _validate_positive(
"correlation_length", correlation_length
)
self.resolution = _validate_positive("resolution", resolution)
self.seed = pin_seed(seed)
def _noise(
self, x: np.ndarray, y: np.ndarray, rng: np.random.Generator
) -> np.ndarray:
"""Perlin noise on the grid of ``x`` and ``y``, as ``(len(y), len(x))``.
Rows are computed in blocks of about :data:`_PERLIN_BLOCK_NODES`
nodes, so the working set beside the output stays bounded. The whole
grid at once held some 24 temporary arrays of its size.
"""
# Perlin noise vanishes at lattice nodes, so a random phase moves the
# output grid off them.
phase = rng.uniform(size=2)
ux = (x - x[0]) / self.correlation_length + phase[0]
uy = (y - y[0]) / self.correlation_length + phase[1]
ix = np.floor(ux).astype(int)
iy = np.floor(uy).astype(int)
tx = ux - ix
ty = uy - iy
gradients = rng.normal(size=(int(iy.max()) + 2, int(ix.max()) + 2, 2))
gradients /= np.linalg.norm(gradients, axis=2, keepdims=True)
# Quintic fade 6t^5 - 15t^4 + 10t^3 of Perlin (2002), "Improving noise",
# ACM Transactions on Graphics 21(3), https://doi.org/10.1145/566654.566636.
fade_x = tx * tx * tx * (tx * (tx * 6.0 - 15.0) + 10.0)
fade_y = ty * ty * ty * (ty * (ty * 6.0 - 15.0) + 10.0)
noise = np.empty((len(y), len(x)))
rows = max(1, _PERLIN_BLOCK_NODES // len(x))
for start in range(0, len(y), rows):
block = slice(start, start + rows)
row_index, row_t = iy[block, None], ty[block, None]
def dot(offset_x: int, offset_y: int) -> np.ndarray:
"""Dot the corner gradient with the offset to each sample point."""
gradient = gradients[row_index + offset_y, ix + offset_x]
return gradient[..., 0] * (tx - offset_x) + gradient[..., 1] * (
row_t - offset_y
)
dot_00, dot_01 = dot(0, 0), dot(0, 1)
lower = dot_00 + fade_x * (dot(1, 0) - dot_00)
upper = dot_01 + fade_x * (dot(1, 1) - dot_01)
noise[block] = lower + fade_y[block, None] * (upper - lower)
return noise
def generate(
self,
boundary: "Boundary",
seed: int | np.random.Generator | None = None,
) -> Surface:
"""Generate a Perlin surface covering ``boundary``'s bounds.
``seed`` overrides the model's seed and is recorded as the integer
``metadata["seed"]``. Values are in meters on a ``(len(y), len(x))``
grid spaced by ``resolution``. Amplitudes do not depend on
``resolution``. Raises ``TypeError`` for a non-``Boundary`` and
``ValueError`` above the allocation budget.
"""
if not isinstance(boundary, Boundary):
raise TypeError("surface generation requires a Boundary")
bounds = boundary.bounds
output_shape = _grid_shape(bounds, self.resolution)
gradient_shape = _grid_shape(bounds, self.correlation_length, extra=2)
_check_grid_budget(
"Perlin surface output grid",
output_shape,
bytes_per_node=_PERLIN_OUTPUT_BYTES_PER_NODE,
parameter=("resolution", self.resolution),
)
_check_grid_budget(
"Perlin gradient lattice",
gradient_shape,
bytes_per_node=_PERLIN_GRADIENT_BYTES_PER_NODE,
parameter=("correlation_length", self.correlation_length),
)
x, y = _grid_axes(bounds, output_shape)
resolved_seed = resolve_seed(self.seed if seed is None else seed)
values = self._noise(x, y, as_rng(resolved_seed))
low, high = self.value_range
# Corner dot products are bounded by sqrt(2) and the fade blend is
# convex, so |noise| <= sqrt(2) maps into the range.
values /= 2.0 * np.sqrt(2.0)
values += 0.5
values *= high - low
values += low
metadata = self.metadata()
metadata.update(
{
"seed": int(resolved_seed),
"bounds": [float(value) for value in bounds],
"grid_shape": [int(output_shape[1]), int(output_shape[0])],
"range_mapping": "bounded_unit_gradient",
}
)
return self._surface_type(x, y, values, metadata=metadata)
def metadata(self) -> dict:
"""Return the JSON-compatible model description."""
return {
"type": type(self).__name__,
self._range_name: list(self.value_range),
"correlation_length": float(self.correlation_length),
"resolution": float(self.resolution),
}
[docs]
class PerlinTerrain(_PerlinSurface):
"""Perlin-noise terrain model with elevations bounded by ``elevation_range``.
Elevations in meters stay within the range without necessarily reaching it.
``correlation_length`` is the noise-cell size in meters, not a statistical
correlation length, and ``resolution`` is the output grid spacing in meters.
Nodes shared by grids of different resolution have equal values for the
same boundary and seed.
"""
_surface_type = Terrain
_range_name = "elevation_range"
def __init__(
self,
elevation_range: tuple[float, float] = (0.0, 100.0),
*,
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
super().__init__(
elevation_range,
correlation_length=correlation_length,
resolution=resolution,
seed=seed,
)
[docs]
class PerlinBathymetry(_PerlinSurface):
"""Perlin-noise bathymetry model with depths bounded by ``depth_range``.
``depth_range`` must be strictly positive, else ``ValueError``. Other
parameters are as for :class:`PerlinTerrain`.
"""
_surface_type = Bathymetry
_range_name = "depth_range"
def __init__(
self,
depth_range: tuple[float, float] = (10.0, 50.0),
*,
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
super().__init__(
depth_range,
correlation_length=correlation_length,
resolution=resolution,
seed=seed,
)
if self.value_range[0] <= 0.0:
raise ValueError("depth_range must contain only positive depths")
# The canonical grid has eight nodes per correlation length. The field is drawn
# by circulant embedding, on a periodic grid padded by eight correlation lengths
# so wrapped covariance adds at most exp(-32), about 1e-14.
_GAUSSIAN_CANONICAL_FRACTION = 8.0
_GAUSSIAN_EMBEDDING_PAD = 64
# Normals, spectrum and transform are complex, 16 bytes a node each.
_GAUSSIAN_EMBEDDING_BYTES_PER_NODE = 4 * np.dtype(complex).itemsize
def _gaussian_embedding_shape(
shape: tuple[int, int], correlation_length: float
) -> tuple[int, int]:
"""The periodic grid embedding a canonical grid of ``shape``, checked to fit.
Raises ``ValueError`` naming what to change when it exceeds the budget.
"""
from scipy import fft
embedding = (
int(fft.next_fast_len(shape[0] + _GAUSSIAN_EMBEDDING_PAD)),
int(fft.next_fast_len(shape[1] + _GAUSSIAN_EMBEDDING_PAD)),
)
_check_grid_budget(
"Gaussian surface canonical grid",
embedding,
bytes_per_node=_GAUSSIAN_EMBEDDING_BYTES_PER_NODE,
parameter=("correlation_length", correlation_length),
)
return embedding
def _periodic_spectrum(nodes: int) -> np.ndarray:
"""Eigenvalues of the squared-exponential covariance wrapped onto ``nodes``.
Nodes lie ``1 / _GAUSSIAN_CANONICAL_FRACTION`` of a correlation length
apart. The covariance summed over its periodic images has, by Poisson
summation, a spectrum of positive Gaussian terms, so the circulant matrix
is a covariance. Its highest frequencies are zero to rounding and are set
to zero.
"""
from scipy import fft
index = np.arange(nodes, dtype=float)
scale = 2.0 * _GAUSSIAN_CANONICAL_FRACTION**2
covariance = sum(
np.exp(-((index - image * nodes) ** 2) / scale) for image in (-1, 0, 1, 2)
)
return np.maximum(fft.fft(covariance).real, 0.0)
def _bilinear_grid(
x_from: np.ndarray,
y_from: np.ndarray,
values: np.ndarray,
x_to: np.ndarray,
y_to: np.ndarray,
) -> np.ndarray:
along_x = np.stack([np.interp(x_to, x_from, row) for row in values], axis=0)
return np.stack(
[np.interp(y_to, y_from, along_x[:, column]) for column in range(len(x_to))],
axis=1,
)
class _GaussianSurface(SurfaceModel):
_surface_type: ClassVar[type[Surface]]
def __init__(
self,
mean: float,
standard_deviation: float,
*,
gradient: ArrayLike = (0.0, 0.0),
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
self.mean = _validate_finite("mean", mean)
self.standard_deviation = _validate_nonnegative(
"standard_deviation", standard_deviation
)
gradient = _validate_real_array("gradient", gradient)
if gradient.shape != (2,) or not np.all(np.isfinite(gradient)):
raise ValueError("gradient must contain two finite values")
self.gradient = tuple(float(value) for value in gradient)
self.correlation_length = _validate_positive(
"correlation_length", correlation_length
)
self.resolution = _validate_positive("resolution", resolution)
self.seed = pin_seed(seed)
def generate(
self,
boundary: "Boundary",
seed: int | np.random.Generator | None = None,
) -> Surface:
"""Generate a Gaussian random-field surface covering ``boundary``'s bounds.
``seed`` overrides the model's seed and is recorded as the integer
``metadata["seed"]``. Values are in meters on a ``(len(y), len(x))``
grid spaced by ``resolution``. The field is drawn on a canonical grid set
by ``correlation_length`` and interpolated, so ``resolution`` does not
change the realization. Raises ``TypeError`` for a non-``Boundary`` and
``ValueError`` above the allocation budget or for a nonpositive
canonical bathymetry depth.
"""
if not isinstance(boundary, Boundary):
raise TypeError("surface generation requires a Boundary")
bounds = boundary.bounds
canonical_spacing = self.correlation_length / _GAUSSIAN_CANONICAL_FRACTION
canonical_shape = _grid_shape(bounds, canonical_spacing)
output_shape = _grid_shape(bounds, self.resolution)
embedding = _gaussian_embedding_shape(canonical_shape, self.correlation_length)
_check_grid_budget(
"Gaussian surface output grid",
output_shape,
bytes_per_node=_GAUSSIAN_OUTPUT_BYTES_PER_NODE,
parameter=("resolution", self.resolution),
)
_check_grid_budget(
"Gaussian surface interpolation grid",
(output_shape[0], canonical_shape[1]),
bytes_per_node=_GAUSSIAN_INTERPOLATION_BYTES_PER_NODE,
parameter=("resolution", self.resolution),
)
# Nodes sit a fixed share of the correlation length apart from the
# bounds' corner. The canonical grid is the corner of a periodic grid
# whose circulant covariance the FFT diagonalizes, and the real part of
# the transform of complex normals scaled by the root of its spectrum
# has that covariance.
from scipy import fft
x = bounds[0] + canonical_spacing * np.arange(canonical_shape[0])
y = bounds[1] + canonical_spacing * np.arange(canonical_shape[1])
resolved_seed = resolve_seed(self.seed if seed is None else seed)
spectrum = np.sqrt(
np.outer(_periodic_spectrum(embedding[1]), _periodic_spectrum(embedding[0]))
/ math.prod(embedding)
)
normals = as_rng(resolved_seed).standard_normal((2, embedding[1], embedding[0]))
transformed = fft.fft2(spectrum * (normals[0] + 1j * normals[1]))
random_field = transformed.real[: len(y), : len(x)]
center_x = 0.5 * (bounds[0] + bounds[2])
center_y = 0.5 * (bounds[1] + bounds[3])
trend = self.mean + self.gradient[0] * (x[None, :] - center_x)
trend = trend + self.gradient[1] * (y[:, None] - center_y)
canonical = trend + self.standard_deviation * random_field
if self._surface_type is Bathymetry and np.any(canonical <= 0.0):
raise ValueError(
"Gaussian bathymetry depths must be positive on the canonical grid"
)
output_x, output_y = _grid_axes(bounds, output_shape)
values = _bilinear_grid(x, y, canonical, output_x, output_y)
metadata = self.metadata()
metadata.update(
{
"seed": int(resolved_seed),
"bounds": [float(value) for value in bounds],
"grid_shape": [int(output_shape[1]), int(output_shape[0])],
"canonical_spacing": float(canonical_spacing),
"canonical_shape": [
int(canonical_shape[1]),
int(canonical_shape[0]),
],
"covariance": "exp(-distance**2/(2*correlation_length**2))",
}
)
return self._surface_type(output_x, output_y, values, metadata=metadata)
def metadata(self) -> dict:
"""Return the JSON-compatible model description."""
return {
"type": type(self).__name__,
"mean": float(self.mean),
"standard_deviation": float(self.standard_deviation),
"gradient": list(self.gradient),
"correlation_length": float(self.correlation_length),
"resolution": float(self.resolution),
}
[docs]
class GaussianTerrain(_GaussianSurface):
"""Gaussian random-field terrain model around a sloping mean plane.
Elevation is ``mean`` plus ``gradient`` times the offset from the center of
the boundary bounds, plus a zero-mean field with squared-exponential
covariance of scale ``standard_deviation`` and length
``correlation_length``, output on a grid spaced by ``resolution``. Lengths
are in meters and ``gradient`` in meters per meter. Negative elevations are
permitted.
"""
_surface_type = Terrain
def __init__(
self,
mean: float = 0.0,
standard_deviation: float = 20.0,
*,
gradient: tuple[float, float] = (0.0, 0.0),
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
super().__init__(
mean,
standard_deviation,
gradient=gradient,
correlation_length=correlation_length,
resolution=resolution,
seed=seed,
)
[docs]
class GaussianBathymetry(_GaussianSurface):
"""Gaussian random-field bathymetry model around a sloping mean depth.
Parameters are as for :class:`GaussianTerrain`, with depth in place of
elevation. The field is never clipped, so a nonpositive depth on the
canonical grid makes ``generate`` raise ``ValueError``.
"""
_surface_type = Bathymetry
def __init__(
self,
mean: float = 50.0,
standard_deviation: float = 5.0,
*,
gradient: tuple[float, float] = (0.0, 0.0),
correlation_length: float = 1_000.0,
resolution: float = 100.0,
seed: int | np.random.Generator | None = None,
):
super().__init__(
mean,
standard_deviation,
gradient=gradient,
correlation_length=correlation_length,
resolution=resolution,
seed=seed,
)
