# geoML - machine learning models for geospatial data
# Copyright (C) 2019 Ítalo Gomes Gonçalves
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Block models: the `_blockdata` fan-out shared by `Blocks1D/2D/3D`, the
variable-size `BlockSet3D` on its integer lattice, `RotatedBlockSet3D`,
and the sub-block geometry behind the mesh assignments and `crossed_by`.
"""
import copy as _copy
import itertools as _iter
import warnings as _warnings
from collections.abc import Sequence
from typing import TYPE_CHECKING
if TYPE_CHECKING:
from geoml.data.geoh5 import Workspace as _GeoH5Workspace
import numpy as _np
import pandas as _pd
import pyvista as _pv
import vtk as _vtk
import geoml._types as _types
import geoml.math.geometry as _gmt
import geoml.storage as _storage
from geoml.data.base import *
from geoml.data.base import _Attribute, _closing_value
from geoml.data.variables import *
from geoml.data.variables import _Variable
from geoml.data.containers import *
from geoml.data.containers import _PointBased, _SpatialData
from geoml.data.grids import *
from geoml.data.grids import (
_TOL_COVER, _aggregate_onto, _cover_box, _fitted_rotation)
from geoml.data.meshes import *
from geoml.data.meshes import (
_below_sheet, _clip_to_box, _mesh_test, _separated, _uncovered_rule,
_within_body)
def _sub_block_shares(blocks, test, rows=None):
"""The share of each block's sub-blocks that `test` accepts.
Walked in chunks: a 5 x 5 x 5 discretization is 125 sub-blocks for
every location, and materializing all of them at once is what the
batching in `get_batched_coordinates` exists to avoid. `rows` narrows it
to the blocks worth asking about; the rest come back `nan`.
"""
per_block = blocks.rows_per_location
chunk = max(1, 1000000 // per_block)
if rows is None:
rows = _np.arange(blocks._n_data)
shares = _np.full(blocks._n_data, _np.nan)
for start in range(0, len(rows), chunk):
index = rows[start:start + chunk]
coordinates, _ = blocks.get_batched_coordinates(index)
# block-major, so every block's sub-blocks are one row
accepted = test(coordinates).reshape([len(index), per_block])
shares[index] = accepted.mean(axis=1)
return shares
def _blocks_from_surface(blocks, base, surface, name, labels, fraction,
uncovered):
"""`assign_from_surface` for anything made of blocks."""
base(blocks, surface, name, labels, uncovered=uncovered)
if fraction is not None:
shares = _sub_block_shares(blocks, _below_sheet(surface))
# a sheet that misses the centre describes the block no better
# than it describes a location, and left the flag empty for the
# same reason -- so the fraction says so rather than reading 0
missed = _np.asarray(blocks.metadata[name].values).ravel() < 0
shares[missed] = _uncovered_rule(uncovered)[1]
blocks.add_metadata(fraction, shares)
def _blocks_from_solid(blocks, base, solid, name, labels, fraction):
"""`assign_from_solid` for anything made of blocks."""
base(blocks, solid, name, labels)
if fraction is not None:
blocks.add_metadata(
fraction, _sub_block_shares(blocks, _within_body(solid)))
def _blockdata(cls):
# Decorator to extend functionality of some classes
# Used to avoid multiple inheritance
old_init = cls.__init__
def new_init(self, start, n, step=None, end=None, labels=None,
discretization=None, **kwargs):
# `end` only where it was given, and anything else straight through:
# a turned grid takes its angles here and has no `end` to be handed
if end is not None:
kwargs["end"] = end
old_init(self, start=start, n=n, step=step, labels=labels, **kwargs)
if discretization is None:
discretization = [1] * self.n_dim
self.discretization = discretization
half = _np.array(self.step_size) / 2
lo = _np.array([axis.min() for axis in self.grid]) - half
hi = _np.array([axis.max() for axis in self.grid]) + half
# a turned box reaches furthest at its corners, where it is turned
corners = _np.array(list(_iter.product(*zip(lo, hi))), dtype=float)
self._bounding_box = BoundingBox.from_array(self._to_world(corners))
sub_grid = _np.array(
list(_iter.product(
*[_np.arange(d) for d in self.discretization[::-1]]
)),
dtype=float)[:, ::-1]
sub_grid -= (_np.array(self.discretization)[None, :] - 1) / 2
sub_grid *= _np.array(self.step_size)[None, :]
sub_grid /= _np.array(self.discretization)[None, :]
self.sub_grid = sub_grid
# `sub_grid` is counted in the grid's own frame, as `grid` is, and turns
# with it on the way out
def discretized_coordinates(self, index):
center = _np.array([g[i] for g, i in zip(self.grid, index)])[None, :]
return self._to_world(self.sub_grid + center)
def inducing_grid(self, index):
center = _np.array([g[i] for g, i in zip(self.grid, index)])[None, :]
grid = self.sub_grid
discr = _np.array(self.discretization)[None, :]
# grid = grid * (discr + 1) / (discr - 1) + center
grid = grid * (discr + 2) / (discr + 1) + center
return PointData.from_array(self._to_world(grid))
def get_batched_coordinates(self, index):
if index is None:
index = _np.arange(self._n_data)
centers = _np.asarray(self.coordinates[index])
# the centres are placed already; an offset from one only turns
offsets = self.sub_grid if self._transform is None \
else _np.matmul(self.sub_grid, self._transform[1])
# block-major, one temporary instead of one per block: block b owns
# rows b*k to (b+1)*k, which is what `_aggregate` averages back
coords = (centers[:, None, :] + offsets[None, :, :])
coords = coords.reshape(-1, centers.shape[1])
splits = None if _np.prod(self.discretization) == 1 else len(index)
return coords, splits
@property
def rows_per_location(self):
return int(_np.prod(self.discretization))
def _batch_rows(self, index=None):
# one row per discretization point of every block in the batch
n, _ = _PointBased._batch_rows(self, index)
n_sub = int(_np.prod(self.discretization))
return n * n_sub, (None if n_sub == 1 else n)
# captured before being replaced, as `old_init` is: the block versions
# only add the partial blocks on top of what the grid already answers
base_assign_from_surface = cls.assign_from_surface
base_assign_from_solid = cls.assign_from_solid
def _sub_block_fraction(self, test):
"""The share of each block's sub-blocks that `test` accepts."""
return _sub_block_shares(self, test)
def assign_from_surface(self, surface, name, labels=("above", "below"),
fraction=None, uncovered=None):
"""
As `Grid3D.assign_from_surface`, measuring the partial blocks on
request.
The flag in `name` follows the block centre, as a whole-block code
does everywhere else. Name a `fraction` column as well and the share
of each block lying below the sheet is measured over the sub-blocks
`discretization` already defines — what a tonnage near surface needs,
where counting a half-buried block whole is the error.
Where the sheet reaches part of a block but not all of it, the
sub-blocks past its edge count as not below. Where it does not reach
the block at all — the centre included, which is what leaves the flag
empty — the fraction is `uncovered` instead of a measurement.
Parameters
----------
surface : Surface3D
The sheet to compare against.
name : str
Name of the metadata column holding the whole-block flag.
labels : tuple
What to call the two sides, in the order (above, below).
fraction : str, optional
Name of a second metadata column, to hold the share of each block
below the sheet. Costs `prod(discretization)` queries per block.
uncovered : float or "raise"
What the `fraction` column records for a block the sheet does not
reach: `numpy.nan` for `None`, the default, so it cannot pass for
a block genuinely above ground, or `0.0` to count it as nothing. Pass
`"raise"` to refuse a surface that does not cover every block.
"""
_blocks_from_surface(self, base_assign_from_surface, surface, name,
labels, fraction, uncovered)
def assign_from_solid(self, solid, name, labels=("outside", "inside"),
fraction=None):
"""
As `_SpatialData.assign_from_solid`, measuring the partial blocks on
request.
`fraction` behaves as it does in `assign_from_surface`: the flag in
`name` follows the block centre, while the column named here holds the
share of each block's sub-blocks falling inside the body — the share
of its volume, for the regular sub-blocks a discretization defines.
Parameters
----------
solid : Surface3D
The closed body to test against.
name : str
Name of the metadata column holding the whole-block flag.
labels : tuple
What to call the two sides, in the order (outside, inside).
fraction : str, optional
Name of a second metadata column, to hold the share of each block
inside the body. Costs `prod(discretization)` queries per block.
"""
_blocks_from_solid(self, base_assign_from_solid, solid, name, labels,
fraction)
cls.__init__ = new_init
cls.discretized_coordinates = discretized_coordinates
cls.inducing_grid = inducing_grid
cls.get_batched_coordinates = get_batched_coordinates
cls._batch_rows = _batch_rows
cls.rows_per_location = rows_per_location
cls._sub_block_fraction = _sub_block_fraction
cls.assign_from_surface = assign_from_surface
cls.assign_from_solid = assign_from_solid
return cls
[docs]
@_blockdata
class Blocks1D(Grid1D):
"""
Blocks along one axis, each averaged over its sub-blocks.
A block model: `Grid1D`'s lattice, each node the centre of a block of
the lattice's spacing.
A prediction is the mean over each block, taken from its `discretization`
sub-blocks, a regular lattice inside it; the variance within the block
is recorded as `dispersion`. What a resource is estimated on when every
block is the same size; `BlockSet3D` is the one whose blocks are not.
"""
[docs]
@_blockdata
class Blocks2D(Grid2D):
"""
Blocks in the plane, each averaged over its sub-blocks.
A block model: `Grid2D`'s lattice, each node the centre of a block of
the lattice's spacing.
A prediction is the mean over each block, taken from its `discretization`
sub-blocks, a regular lattice inside it; the variance within the block
is recorded as `dispersion`. What a resource is estimated on when every
block is the same size; `BlockSet3D` is the one whose blocks are not.
"""
[docs]
@_blockdata
class Blocks3D(Grid3D):
"""
Blocks in space, each averaged over its sub-blocks.
A block model: `Grid3D`'s lattice, each node the centre of a block of
the lattice's spacing.
A prediction is the mean over each block, taken from its `discretization`
sub-blocks, a regular lattice inside it; the variance within the block
is recorded as `dispersion`. What a resource is estimated on when every
block is the same size; `BlockSet3D` is the one whose blocks are not.
"""
[docs]
def as_pyvista(self, simulations=False, include="**"):
"""
Converts this object to a pyvista one, carrying its variables.
Parameters
----------
simulations
Which simulations to include: `False` for none (the default, since
each one is a full-length array in the exported object), `True` for
all of them, an `int` for the first n, or a sequence of indices.
"""
pv_blocks = _pv.ImageData(
dimensions=_np.array(self.grid_size) + 1,
spacing=self.step_size,
origin=_np.array([ax[0] for ax in self.grid]) - _np.array(self.step_size) / 2
)
return self._finish_pyvista(pv_blocks, "blocks", simulations,
include)
[docs]
def to_geoh5(self, workspace: "_types.PathLike | _GeoH5Workspace",
name: str = "Blocks", include: str = "**",
simulations: "bool | int | Sequence[int]" = False,
replace: bool = True,
folder: "str | None" = None) -> None:
"""
Writes this block model into a geoh5 workspace, as a BlockModel.
The geoh5 BlockModel is the regular-grid sibling of the Octree a
`BlockSet3D` writes; a uniform model maps onto it exactly.
Everything else is as `PointData.to_geoh5`: the workspace
accumulates, the columns ride with their path table, `replace`
makes a re-run export mean what it says. Needs the `geoh5py`
package: `pip install geoml[geoh5]`.
Parameters
----------
workspace
Path of the workspace to write into, or an open
`geoml.data.geoh5.Workspace`.
name
The name the object gets in the workspace.
include
A path pattern selecting the columns to carry.
simulations
Which simulations to include: `False` for none (the default),
`True` for all, an `int` for the first n, or a sequence of
indices.
replace
Whether an existing BlockModel of this name, in this folder,
makes way.
folder
Where the object sits in ANALYST's project tree, as a path —
each segment a group, created when it does not exist and
reused when it does. `None` is the root.
"""
import geoml.data.geoh5 as _geoh5io
_geoh5io.write_grid_blocks(self, workspace, name, include,
simulations, replace, folder)
[docs]
@classmethod
def from_geoh5(cls, workspace: "_types.PathLike | _GeoH5Workspace",
name: "str | None" = None
) -> "Blocks3D | RotatedBlockSet3D":
"""
Reads a BlockModel from a geoh5 workspace.
Only a **uniform** spacing has a geoML container: a true tartan
grid — uneven cells along an axis — is refused with the axis
named. An unrotated model comes back as a `Blocks3D`; a rotated
one as a `RotatedBlockSet3D` with `max_levels=0` — the blocks,
their rotation and block-support prediction all work, nothing can
be refined, and `as_blocks3d` makes it a `RotatedBlocks3D`.
Float cell data becomes continuous variables and referenced data
categorical ones, as in `PointData.from_geoh5`. Needs the
`geoh5py` package: `pip install geoml[geoh5]`.
Parameters
----------
workspace
Path of the workspace to read, or an open
`geoml.data.geoh5.Workspace`.
name
The BlockModel to read, when the file holds more than one.
Returns
-------
blocks : Blocks3D or RotatedBlockSet3D
Whichever the file's rotation calls for.
Raises
------
ValueError
If the workspace holds no such BlockModel, or its spacing is
uneven.
"""
import geoml.data.geoh5 as _geoh5io
raw = _geoh5io.read_grid_blocks(workspace, name)
azimuth = -raw["rotation"]
matrix = _gmt.rotation_matrix(azimuth, 0.0, 0.0)
start = (raw["offset"] + raw["step"] / 2.0) @ matrix \
+ raw["corner"]
if raw["rotation"] == 0.0:
blocks = Blocks3D(start, raw["shape"], raw["step"])
else:
blocks = RotatedBlockSet3D(start, raw["shape"], raw["step"],
azimuth=azimuth, max_levels=0)
for column, values in raw["floats"]:
blocks.add_continuous_variable(column, values)
for column, labels, values in raw["coded"]:
values = _np.asarray(values, dtype=object)
values[values == ""] = None
blocks.add_categorical_variable(column, labels=labels,
measurements=values)
return blocks
[docs]
@_blockdata
class RotatedBlocks3D(RotatedGrid3D):
"""
A regular block model turned about its first block.
`Blocks3D` with an azimuth, a dip and a rake, as `RotatedGrid3D` is
`Grid3D` with them. The blocks are counted in the model's own frame and
turned about the first block's centre where coordinates leave -- the
centres, the sub-blocks a prediction averages over, the exported cells
-- and turned back where they come in (`index_data`, and so
`aggregate`). `BlockSet3D.as_blocks3d` returns one for a
`RotatedBlockSet3D`, with the same angles and the same pivot.
Parameters
----------
start : array-like
Centre of the first block, which the model turns about.
n : array-like
Number of blocks along each of the model's own axes.
step : array-like
Size of a block along each of them.
labels : list, optional
Coordinate names.
discretization : array-like, optional
Sub-blocks per axis, averaged over to predict a block. One by
default, which predicts each block at its centre.
azimuth, dip, rake : float
The rotation in degrees, as `RotatedGrid3D` takes it. Given by
keyword.
"""
[docs]
def as_pyvista(self, simulations=False, include="**"):
"""
Converts this object to a pyvista one, carrying its variables.
The cells of a `Blocks3D` of the same shape, turned into place as
one piece, so that each cell's centre is its block's.
Parameters
----------
simulations
Which simulations to include: `False` for none (the default, since
each one is a full-length array in the exported object), `True` for
all of them, an `int` for the first n, or a sequence of indices.
"""
pv_blocks = _pv.ImageData(
dimensions=_np.array(self.grid_size) + 1,
spacing=self.step_size,
origin=_np.array([ax[0] for ax in self.grid])
- _np.array(self.step_size) / 2
)
return self._turned(self._finish_pyvista(
pv_blocks, "blocks", simulations, include))
def _ghost_shell(origin, size, shape):
"""Mirror images of the boundary cells, in all 26 directions.
A ghost is its partner's reflection across the face, edge or corner they
share, so their corners coincide exactly and the corner averaging cannot
tear along the box surface. What each holds is up to the caller: a copy
of its partner's value carries the field past the box as it stands at
the box, for the box to cut the surface off exactly, and `_ghost_values`
puts the cap on the box faces instead, where nothing will cut it.
**The diagonals are not optional** for the second, though they read as
though they should be: nothing crosses a box edge from inside, but the
cap itself runs along the boundary, and where two face ghosts meet
without the diagonal between them the cap reaches the shell's own edge
and stops there. That is a hole in the body, and it is the common case
rather than a corner one -- `close="below"` keeps the region under the
value, which on most models touches every face, so its cap ran the
length of all twelve box edges and tore along each.
"""
ghost_origin, ghost_size, ghost_parent = [], [], []
for direction in _iter.product((-1, 0, 1), repeat=3):
if not any(direction):
continue
touching = _np.ones(len(origin), dtype=bool)
for axis, step in enumerate(direction):
if step < 0:
touching &= origin[:, axis] == 0
elif step > 0:
touching &= origin[:, axis] + size[:, axis] == shape[axis]
if not touching.any():
continue
mirrored = origin[touching].copy()
for axis, step in enumerate(direction):
if step:
mirrored[:, axis] += step * size[touching][:, axis]
ghost_origin.append(mirrored)
ghost_size.append(size[touching])
ghost_parent.append(_np.flatnonzero(touching))
if not ghost_origin:
return (_np.zeros([0, 3], dtype=origin.dtype),
_np.zeros([0, 3], dtype=size.dtype),
_np.zeros([0], dtype=int))
return (_np.concatenate(ghost_origin), _np.concatenate(ghost_size),
_np.concatenate(ghost_parent))
def _ghost_values(parent_values, value, keep_above):
"""What each ghost holds, so the cap lands on the box face.
The way a closed contour is capped when the box cannot cut it: where
Manifold will not take the surface carried past the box as a body, and
in the welded mesh, which has no padding to close against. Every face
point of a kept block then sits exactly at the level, which rounds the
body's edges on the box by about half a boundary block.
A ghost is the mirror of the block it stands against, and the value it
carries decides where between the two centres the surface crosses. The
**reflection about the contour level**, `2 * value - v`, puts that
crossing exactly halfway -- on the shared face, which is the box
boundary. One shared constant far past the data's range put it hard
against the block's own centre instead, half a block short of the face,
so a body was cut back on the side it was kept and let out on the other.
Only the blocks the body occupies are reflected. A boundary block on the
far side of the level set has no cap to make and keeps its value, so
nothing crosses between it and its ghost. Either way every ghost ends up
outside the kept region, which is what stops the surface reaching the
shell's own outer face and leaving the body open.
Measured against Monte Carlo volumes on a ball meeting the box at a
face, an edge, a corner, and half outside it: `close="above"` was
reading -6.3%, -12.9%, -20.0% and -13.9%, and now reads within 1%;
`close="below"`, the complement, was +2.5% to +5.1% and is now within
0.2%. The worst of the eight went from 20.0% to 0.99%.
"""
kept = parent_values > value if keep_above else parent_values < value
return _np.where(kept, 2.0 * value - parent_values, parent_values)
def _corner_keys(origin, size, span):
"""Each block's eight corners as integer keys of the lattice points they
fall on, z-major, `span` points along each axis."""
keys = _np.empty((len(origin), 8), dtype=_np.int64)
for j, corner in enumerate(_gmt.HEX_CORNERS):
at = origin + corner * size
keys[:, j] = (at[:, 2] * span[1] + at[:, 1]) * span[0] + at[:, 0]
return keys
def _lattice_corners(origin, size, span):
"""The distinct lattice points the blocks' corners fall on, and each
block's eight corners as indices into them.
`origin` and `size` count lattice cells and `span` is the number of
points along each axis. The points come back as sorted z-major keys, so
one z-slab of them is one contiguous run. This is the welding of
`_hex_mesh`, done in integers, and on a big model it was most of a
contour's memory: the corners as an `n x 8 x 3` int64 array, their
keys, and `np.unique`'s sorted copy, argsort and int64 inverse, over
the 7 to 12 million cells of a Tom v6 contour, took 4 to 8 GB. Here the
keys are built a corner at a time, sorted once, and each key's rank
among the distinct points scattered back through the sort into 32-bit
integers wherever the keys number fewer than 2**31, every temporary
freed as the next one exists. Finding the ranks by `searchsorted` into
the distinct points instead is leaner still and was measured 17 times
slower than `np.unique`, the argsort 1.4 times faster.
"""
origin = _np.asarray(origin, dtype=_np.int64)
size = _np.asarray(size, dtype=_np.int64)
span = _np.asarray(span, dtype=_np.int64)
keys = _corner_keys(origin, size, span)
rank = _np.int32 if keys.size < 2 ** 31 else _np.int64
if keys.size == 0:
return _np.zeros(0, dtype=_np.int64), _np.zeros((0, 8), dtype=rank)
flat = keys.ravel()
order = _np.argsort(flat)
ordered = flat[order]
del keys, flat
fresh = _np.empty(len(ordered), dtype=bool)
fresh[0] = True
_np.not_equal(ordered[1:], ordered[:-1], out=fresh[1:])
distinct = ordered[fresh]
del ordered
group = _np.cumsum(fresh, dtype=rank)
group -= 1
del fresh
index = _np.empty(len(group), dtype=rank)
index[order] = group
return distinct, index.reshape(-1, 8)
# the largest slab _painted_contour will hold at once, in cells: two of
# these in float64 is ~130 MB, transient
_PAINT_BUDGET = 8_000_000
# how many interior points of coarse blocks it fills at once: 32-bit
# coordinates and float64 values, at most ~160 MB transient with the mask
# and the kept copies
_FILL_BUDGET = 4_000_000
def _painted_contour(origin, size, step, values, corner, value, label,
background, margin=0, fields=None):
"""The welded field painted onto a regular grid, contoured by flying
edges.
What `_cut_to_contour` hands over is a list of axis-aligned cells on
one lattice — origins and sizes in whole cells — and turning that list
into an unstructured mesh was most of a contour's cost: welding eight
corners per cell (`unique` over millions of rows) and VTK's
single-threaded unstructured contour. Measured at 912k blocks, 9.1 s
of an 11.9 s call against 3.7 painted.
What is painted is the *point* field the welded mesh reads, never the
cells: the cell-equal mean over the blocks meeting at each block
corner (`_at_corners`' rule — one vote per block, whatever its size),
and inside any block larger than one cell the trilinear reading of
its own eight corner means, which is exactly what VTK interpolates
across a hexahedron — and trilinear survives subdivision, so the fine
lattice reproduces the coarse cell's surface rather than resampling
it. Painting cell values and letting `cell_data_to_point_data`
average them was the first version of this route, and it disagrees
with the welded mesh at every fine/coarse interface (a
volume-weighted corner against a cell-equal one); on a near-binary
indicator, contoured within a whisker of its corners, that
disagreement pinched the surface into a few dozen same-way edges
nothing downstream could settle, so such fields had to be routed
around it. One field, one surface, no routing.
Painted in z-slabs of at most `_PAINT_BUDGET` cells. The point plane
two slabs share is painted by both from the same corner table,
identically, so the pieces weld back seamlessly and no ghost layers
are needed. Where two blocks' interior fills meet at a point that is
no block's corner — a coarse/coarse interface — the two trilinear
readings can differ and the later paint wins, but `_cut_to_contour`
keeps the surface a margin of finest cells away from any such point,
which is the same guarantee the welded mesh itself rests on.
`background` fills whatever the blocks do not cover. Without closing
ghosts the blocks tile their box exactly and it never shows; with them
the painted box is the hull of ghosts of different sizes, and the
hollows beyond the smaller ghosts must read as far outside for the cap
to stay shut. A block holding no value contributes to no corner
(`_at_corners`' rule again), and where a valued block has each of its
corners it is read across from them, as the welded mesh reads a
valueless hexahedron. What else it covers is absent ground, with no
level to draw, and flying edges put a crossing between a value and
none at NaN. A contour closing against `background` closes against it
too, painted far past the level on the same side, within a hundredth
of a cell of the last valued corner; an open one draws nothing in a
cell any of whose corners is absent, and stops there as at the box.
`margin` pads the painted box with that many cells of background on
every side, so a field carried out to the edge of the blocks still
closes before the image ends.
With `fields`, `values` holds several columns, all averaged onto the
corners, and `fields` turns those means into the fields contoured --
an array of one column per field, point by point -- painted and
contoured together, each slab `1 / n` of the budget.
Returns
-------
verts, faces : arrays
The contoured triangulation, in the lattice's own frame — a
rotated set turns the vertices afterwards. Empty when the surface
misses the value. With `fields`, a list of them, one per field.
"""
origin = _np.asarray(origin, dtype=_np.int64)
size = _np.asarray(size, dtype=_np.int64)
values = _np.asarray(values, dtype=float)
step = _np.broadcast_to(_np.asarray(step, dtype=float), (3,))
low = origin.min(axis=0) - int(margin)
dims = ((origin + size).max(axis=0) + int(margin) - low).astype(_np.int64)
origin = origin - low
world = _np.asarray(corner, dtype=float) + low * step
span = dims + 1
plane = int(span[0]) * int(span[1])
# the corner table: every finite block's eight corners as one integer
# key — z-major, so one slab's point planes are one contiguous key
# range — averaged with one vote per block
known = _np.isfinite(values) if values.ndim == 1 \
else _np.all(_np.isfinite(values), axis=1)
corner_key, inverse = _lattice_corners(origin[known], size[known], span)
votes = _np.bincount(inverse.ravel())
if fields is None:
corner_value = (_np.bincount(inverse.ravel(),
weights=_np.repeat(values[known], 8))
/ votes)[:, None]
else:
corner_value = _np.asarray(fields(_np.stack(
[_np.bincount(inverse.ravel(),
weights=_np.repeat(values[known, k], 8)) / votes
for k in range(values.shape[1])], axis=1)))
del votes
n_fields = corner_value.shape[1]
# a finite block bigger than one cell owns points no corner pass will
# visit; it reads its eight corner means back from the table (its own
# vote keeps them finite) to spread trilinearly. A one-cell block has
# no such points and needs no fill at all. A block holding no value
# fills the same way where the table has all eight of its corners,
# and paints its points absent, before any fill, where it does not.
big = _np.any(size > 1, axis=1)
big_rows = _np.flatnonzero(known & big)
eight = corner_value[inverse[big[known]]]
del inverse
absent_rows = _np.flatnonzero(~known)
if len(absent_rows) and len(corner_key):
keys = _corner_keys(origin[absent_rows], size[absent_rows], span)
at = _np.minimum(_np.searchsorted(corner_key, keys),
len(corner_key) - 1)
reached = _np.all(corner_key[at] == keys, axis=1)
bridged = reached & big[absent_rows]
big_rows = _np.concatenate([big_rows, absent_rows[bridged]])
eight = _np.concatenate([eight, corner_value[at[bridged]]])
absent_rows = absent_rows[~reached]
del keys, at
row_in_big = _np.full(len(values), -1, dtype=_np.int64)
row_in_big[big_rows] = _np.arange(len(big_rows))
def size_classes(rows):
# sizes as one integer key: `unique` over an axis sorts rows as void
# records, which measured ~1.5 s of this call at 912k blocks where
# the keyed form is milliseconds
size_key = (size[rows, 0] + (size[rows, 1] << 21)
+ (size[rows, 2] << 42))
_, first, member = _np.unique(size_key, return_index=True,
return_inverse=True)
return rows, size[rows[first]], member.ravel()
# the absent ground first, so a fill sharing a point with it wins
passes = [size_classes(absent_rows) + (False,),
size_classes(big_rows) + (True,)]
thick = max(1, int(_PAINT_BUDGET // max(1, plane * n_fields)))
# absent ground reads a hundred times further past the level than any
# corner, on the side a closing contour leaves out, which puts the
# body's face against it within a hundredth of a cell of the corners
closing = bool(_np.isfinite(background))
gap = abs(background - value) if closing else 1.0
reach = _np.maximum(
_np.max(corner_value, axis=0, initial=value) - value,
value - _np.min(corner_value, axis=0, initial=value))
far = value + (_np.sign(background - value) if closing else -1.0) * (
100.0 * reach + gap)
# the fills' coordinates in 32 bits wherever the painted box allows
lattice = origin.astype(_np.int32) if int(span.max()) < 2 ** 31 \
else origin
verts = [[] for _ in range(n_fields)]
faces = [[] for _ in range(n_fields)]
count = [0] * n_fields
for z_start in range(0, int(dims[2]), thick):
z_stop = min(z_start + thick, int(dims[2]))
paint = _np.full((n_fields, z_stop - z_start + 1, int(span[1]),
int(span[0])), background, dtype=float)
for fill_rows, classes, member, fills in passes:
for index, shape in enumerate(classes):
rows = fill_rows[
(member == index)
& (origin[fill_rows, 2] <= z_stop)
& (origin[fill_rows, 2] + size[fill_rows, 2] >= z_start)]
if len(rows) == 0:
continue
offsets = _np.stack(_np.meshgrid(
_np.arange(shape[0] + 1), _np.arange(shape[1] + 1),
_np.arange(shape[2] + 1), indexing="ij"),
axis=-1).reshape(-1, 3)
t = offsets / shape
mix = _np.prod(_np.where(_gmt.HEX_CORNERS[None, :, :] == 1,
t[:, None, :], 1.0 - t[:, None, :]),
axis=2)
offsets = offsets.astype(lattice.dtype)
# a few rows at a time: every point of every block of a
# size class at once -- int64 coordinates, their values,
# the mask and both kept copies -- was the paint's peak,
# well over its corner table (729 points a block for a
# coarse block eight cells a side)
batch = max(1, _FILL_BUDGET // (len(offsets) * n_fields))
for first_row in range(0, len(rows), batch):
part = rows[first_row:first_row + batch]
filled = _np.full((n_fields, len(part), len(offsets)),
_np.nan)
if fills:
reading = eight[row_in_big[part]]
for m in range(n_fields):
filled[m] = reading[:, :, m] @ mix.T
points = (lattice[part, None, :]
+ offsets[None, :, :]).reshape(-1, 3)
filled = filled.reshape(n_fields, -1)
keep = ((points[:, 2] >= z_start)
& (points[:, 2] <= z_stop))
points = points[keep]
paint[:, points[:, 2] - z_start, points[:, 1],
points[:, 0]] = filled[:, keep]
# the corner means last, over whatever the fills wrote: every
# point that is any block's corner reads the welded value
lo = _np.searchsorted(corner_key, z_start * plane)
hi = _np.searchsorted(corner_key, (z_stop + 1) * plane)
keys = corner_key[lo:hi]
rest = keys % plane
paint[:, keys // plane - z_start, rest // span[0],
rest % span[0]] = corner_value[lo:hi].T
# The image lives in the lattice's own units, and the world enters
# only after the contour, in float64: every VTK image contour
# writes float32 points with no say in the matter (flying edges
# has no output-precision setting at all), and at mine-grid
# coordinates that quantizes a northing of 2.8e6 to steps of 0.25
# and snaps any crossing within ~1e-4 of a lattice plane onto it,
# welding the surface into a pinch wherever it grazes one --
# measured twice on a real model, both at a bench boundary the
# surface hugged. At lattice magnitudes float32 resolves ~1e-5 of
# one cell, and a crossing collapsing onto a plane from that close
# is a sliver the degenerate-face drop already owns. The welded
# mesh never had the problem because an unstructured contour
# inherits its input's float64.
for m in range(n_fields):
absent = _np.isnan(paint[m])
holes = bool(absent.any())
if holes:
paint[m][absent] = far[m]
image = _pv.ImageData(
dimensions=(int(span[0]), int(span[1]),
z_stop - z_start + 1),
spacing=(1.0, 1.0, 1.0),
origin=(0.0, 0.0, float(z_start)))
image.point_data[label] = paint[m].ravel()
edges = _vtk.vtkFlyingEdges3D()
edges.SetInputData(image)
edges.SetValue(0, value)
edges.ComputeNormalsOff()
edges.ComputeGradientsOff()
edges.ComputeScalarsOff()
edges.Update()
piece = _pv.wrap(edges.GetOutput())
if piece.n_cells:
piece = piece.triangulate()
points = _np.asarray(piece.points, dtype=float)
triangles = piece.faces.reshape(-1, 4)[:, 1:]
if holes and not closing:
points, triangles = _clear_of(absent, points, triangles,
z_start)
if len(triangles):
verts[m].append(world + points * step)
faces[m].append(triangles + count[m])
count[m] += len(points)
found = [(_np.concatenate(verts[m]), _np.concatenate(faces[m]))
if verts[m] else (_np.zeros([0, 3]), _np.zeros([0, 3], dtype=int))
for m in range(n_fields)]
return found[0] if fields is None else found
def _clear_of(absent, points, triangles, z_start):
"""The triangles of a painted slab drawn in cells with no absent corner,
and the points they use.
Flying edges draws each triangle inside one cell, so the cell holding
its centroid is the one that drew it, and a cell with an absent corner
has no level to draw whatever its other corners read. `absent` marks
the slab's points, z first, and `points` are in the image's own frame,
whose slab starts at `z_start`.
"""
lost = _np.zeros([n - 1 for n in absent.shape], dtype=bool)
for dz, dy, dx in _iter.product((0, 1), repeat=3):
lost |= absent[dz:dz + lost.shape[0], dy:dy + lost.shape[1],
dx:dx + lost.shape[2]]
cell = _np.floor(points[triangles].mean(axis=1)).astype(_np.int64)
cell[:, 2] -= z_start
cell = _np.clip(cell, 0, _np.array(lost.shape[::-1]) - 1)
triangles = triangles[~lost[cell[:, 2], cell[:, 1], cell[:, 0]]]
used = _np.unique(triangles)
index = _np.zeros(len(points), dtype=triangles.dtype)
index[used] = _np.arange(len(used))
return points[used], index[triangles]
def _contour_column(blocks, path) -> "tuple[VariablePath, _Attribute]":
"""The column a contour path names, and that path in full.
A path naming a column (`"metals/zn/prediction"`) is taken as it stands;
one naming a variable or component defaults to its `prediction`. A bare
segment that resolves to nothing on its own is searched for anywhere in
the tree, which is what lets a component be named without spelling its
owner -- refused when the name belongs to more than one place.
"""
path = VariablePath(str(path))
try:
found = blocks.get(path)
except KeyError as err:
if len(path.parts) != 1:
raise ValueError(str(err))
matches = {p: thing
for p, thing in blocks.select("**/" + path.parts[0]).items()
if not isinstance(thing, _Attribute)}
if len(matches) > 1:
raise ValueError(
"%r sits in more than one place: %s; give the full path"
% (str(path), ", ".join(str(p) for p in matches)))
if len(matches) == 0:
# raises with the canonical listing of what there is to name
blocks._variable_or_component(path.parts[0])
raise ValueError(
"no variable or component named %r" % str(path))
(path, found), = matches.items()
if isinstance(found, _Attribute):
return path, found
try:
column = blocks.get(path / "prediction")
except KeyError:
parts = getattr(found, "components", None) or {}
if parts:
raise ValueError(
"%r is made of components and holds no prediction of its "
"own; name one of %s"
% (str(path), ", ".join(str(name) for name in parts)))
raise ValueError("%r holds no prediction to contour" % str(path))
return path / "prediction", column
def _add_rows(target, at, values):
"""`target[at] += values`, rows bound for the same place summed first."""
order = _np.argsort(at, kind="stable")
at = at[order]
starts = _np.flatnonzero(_np.r_[True, at[1:] != at[:-1]])
target[at[starts]] += _np.add.reduceat(values[order], starts, axis=0)
class _Grouping:
"""The blocks of a set gathered into the blocks of its coarsest level.
`cell` names the coarse block each block falls in and `volume` how many
base cells it holds, so its share of that block (`weight`) is a ratio of
integers. A coarse block that is one of the set's blocks whole is `whole`
and keeps it exactly -- `source` says which. The rest are gathered from
their parts, and what that means for a column is the column's to say
(`_Variable._coarsen_into`); these are the operations it says it with.
A value missing in any part is missing in the whole, which is what keeps
a partial answer from passing for one.
"""
def __init__(self, cell, volume, cell_volume, n_cells):
self.cell = cell
self.volume = volume
self.weight = volume / float(cell_volume)
self.n_cells = n_cells
self._alone = volume == cell_volume
self.source = _np.full(n_cells, -1, dtype=_np.int64)
self.source[cell[self._alone]] = _np.flatnonzero(self._alone)
self.whole = self.source >= 0
def mean(self, values, valid=None):
"""The volume-weighted mean of a column over each coarse block.
`valid` marks the parts holding a value, for a column that cannot
say so itself; a block kept whole is taken as it stands.
"""
values = _np.asarray(values, dtype=float)
if valid is not None:
values = _np.where(valid | self._alone, values, _np.nan)
return _np.bincount(self.cell, weights=self.weight * values,
minlength=self.n_cells)
def kept(self, values, missing):
"""A column where a coarse block is one block whole, `missing`
where it is gathered from parts."""
values = _np.asarray(values)
out = _np.full(self.n_cells, missing, dtype=values.dtype)
out[self.whole] = values[self.source[self.whole]]
return out
def realizations(self, store, target, valid=None, shift=None):
"""Each realization's volume-weighted mean, written into `target`.
Realization `i` of a coarse block is the mean of realization `i` of
its parts: pairing them any other way would invent a correlation the
model never produced. The source is read a band of rows at a time
and gathered into as many coarse blocks as a working array of the
store's spill size holds, so neither store is held whole; more coarse
blocks than that cost another pass over the source each.
With `shift` -- one number per coarse block, near its mean -- the
same pass returns how far the parts sit from their block, the
variance between them averaged over the realizations. The moments
are taken about the shift so that their difference keeps its digits,
and a block kept whole comes out at exactly zero.
"""
n_sim = int(store.shape[1])
# the sum, and with a shift the two moments about it
moments = 1 if shift is None else 3
rows = max(1, _storage.DEFAULT_THRESHOLD
// (8 * moments * max(n_sim, 1)))
keep = None if valid is None else valid | self._alone
spread = None if shift is None else _np.full(self.n_cells, _np.nan)
bands = store.row_bands()
for lo in range(0, self.n_cells, rows):
hi = min(lo + rows, self.n_cells)
sums = _np.zeros((moments, hi - lo, n_sim))
for band in bands:
cell = self.cell[band]
inside = (cell >= lo) & (cell < hi)
if not _np.any(inside):
continue
values = _np.asarray(store[band], dtype=float)[inside]
if keep is not None:
values[~keep[band][inside]] = _np.nan
weight = self.weight[band][inside][:, None]
at = cell[inside] - lo
_add_rows(sums[0], at, weight * values)
if shift is not None:
offset = values - shift[cell[inside]][:, None]
_add_rows(sums[1], at, weight * offset)
_add_rows(sums[2], at, weight * offset * offset)
target[lo:hi, :] = sums[0]
if spread is not None:
spread[lo:hi] = _np.mean(
_np.maximum(sums[2] - sums[1] * sums[1], 0.0), axis=1)
return spread
def gathered(self, column):
"""A metadata column over the coarse blocks, gathered as `aggregate`
gathers one but weighed by volume: a coded column keeps the label
holding most of the block, none on a tie; a boolean one holds where
it held throughout; a number averages."""
values = column.values.to_numpy()
if column.labels is not None:
return self._dominant(values)
if values.dtype == bool:
out = self.kept(values, False)
broken = _np.bincount(self.cell, weights=(~values).astype(float),
minlength=self.n_cells) > 0
out[~self.whole] = ~broken[~self.whole]
return out
return self.mean(values)
def _dominant(self, codes):
"""The code holding most of each gathered block's volume, counted in
base cells so that a tie is exact; missing codes do not vote."""
out = self.kept(codes, -1)
voting = ~self.whole[self.cell] & (codes >= 0)
if not _np.any(voting):
return out
votes = _pd.DataFrame({"cell": self.cell[voting],
"code": codes[voting],
"volume": self.volume[voting]})
votes = votes.groupby(["cell", "code"])["volume"].sum().reset_index()
most = votes.groupby("cell")["volume"].transform("max")
winners = votes[votes["volume"] == most].groupby("cell")["code"].agg(
["first", "size"])
settled = winners[winners["size"] == 1]
out[settled.index.to_numpy()] = settled["first"].to_numpy()
return out
[docs]
class BlockSet3D(PointData):
"""
Blocks of several sizes, on one integer lattice.
A block model where the interesting ground can be carried finely and the
rest coarsely. On a real deposit the ground worth resolving is a small
part of the volume, and a uniform model at the resolution that part needs
spends almost all of its cells saying nothing: refining 5 m only where it
is wanted takes a 29-million-cell model to under 700 000.
Every block's position and size are whole numbers of a **base cell**, the
finest the model may go, which is `step / discretization ** max_levels`.
Working in those integers rather than in metres is what makes it exact:
blocks meet without a tolerance, a block is a whole number of its own
children, and regrouping conserves mass to the last digit. It also means
the model can say which of its answers rest on coarse blocks, which is the
one thing a mixed-support model has to be able to prove -- see
`docs/variable-block-models.md`.
It is built **full**: the blocks tile their box exactly, and every
operation keeps them that way. Ground to leave out is filtered, not
removed, so that grouping is always safe -- a half-populated group would
average over blocks that are not there and quietly weigh the answer wrong.
`discretization` does two jobs, and they are the same job. It is how
finely a block is sampled to average it, and it is how a block splits:
each sub-block becomes a child. So the refinement ratio is the
discretization, per axis and not necessarily two -- `[2, 2, 1]` refines in
plan and leaves the bench height alone. Being the same at every level is
what lets a block of any size fan out into the same number of rows, so the
model sees one shape whatever it is looking at and nothing downstream has
to know that levels exist. It costs a coarse block some accuracy in its
own average, always by overstating how variable it is, which errs towards
splitting it -- and splitting is what removes the error.
Parameters
----------
start : array-like
Centre of the first (coarsest) block.
n : array-like
Number of blocks along each axis, at the coarsest level.
step : array-like
Size of a coarsest-level block.
discretization : array-like
Sub-blocks per axis, used at every level, and so also the ratio a
block splits by. An axis given 1 is never refined.
max_levels : int
How many times a block may be split. Fixes the base cell, and so the
lattice everything else is counted in.
labels : list
Coordinate names.
Attributes
----------
block_size : array
`(n_data, 3)`, each block's size in the coordinates' own units. Not
called `step_size`: that name means one size for the whole object, and
anything reading it would take the product of this array for a volume.
block_volume : array
`(n_data,)`, what each block is worth in a tonnage.
level : array
`(n_data,)`, 0 for a coarsest block up to `max_levels` for a base one.
"""
def __init__(self, start, n, step, discretization=(2, 2, 2), max_levels=3,
labels=("X", "Y", "Z")):
if int(max_levels) < 0:
raise ValueError("max_levels cannot be negative; got %r"
% max_levels)
self.max_levels = int(max_levels)
self.discretization = [int(d) for d in discretization]
if len(self.discretization) != 3 or any(
d < 1 for d in self.discretization):
raise ValueError(
"discretization needs three positive numbers; got %r"
% (discretization,))
if self.max_levels > 0 and _np.prod(self.discretization) == 1:
raise ValueError(
"a discretization of one sub-block cannot refine anything, so "
"max_levels=%d would have nothing to do; give a coarser "
"discretization or max_levels=0" % self.max_levels)
n = _np.asarray(n, dtype=_np.int64)
step = _np.asarray(step, dtype=float)
start = _np.asarray(start, dtype=float)
if not (n.shape == step.shape == start.shape == (3,)):
raise ValueError(
"start, n and step are three numbers each; got %s, %s and %s"
% (start.shape, n.shape, step.shape))
# Splitting a block hands each of its sub-blocks a block of its own, so
# the discretization is also the refinement ratio -- per axis, and not
# necessarily two. That is what keeps sub-block `j` and child `j` the
# same corner of the parent, which is what lets the criterion read off
# a coarse prediction mean anything about the children it would make.
# It also lets an axis be left alone: [2, 2, 1] refines in plan and
# keeps the bench height.
self.base_step = step / self._coarse_size
# the box corner, not the first centre: a block's origin is its lower
# corner, which is what makes the lattice arithmetic come out integer
self.box_corner = start - step / 2
self.lattice_shape = n * self._coarse_size
cells = _np.stack(_np.meshgrid(
*[_np.arange(k) for k in n], indexing="ij"), axis=-1
).reshape(-1, 3)
self._origin = cells * self._coarse_size
self._level = _np.zeros(len(cells), dtype=_np.int64)
# nothing has been predicted yet, so every block is new
self._fresh = _np.ones(len(cells), dtype=bool)
self._sub_grid = _gmt.unit_sub_grid(self.discretization)
super().__init__(
_pd.DataFrame(self._centres(), columns=list(labels)), list(labels))
self._bounding_box = self._lattice_bounding_box()
# ------------------------------------------------------------------ #
# geometry
# ------------------------------------------------------------------ #
def _centres(self):
return self.box_corner + (self._origin + self._size / 2) \
* self.base_step
@property
def _coarse_size(self):
"""A coarsest block's extent, in base cells, along each axis.
Not the number of sub-blocks, which is `prod(discretization)` and is
the same at every level -- this is how far a level-0 block reaches
across the lattice, and so how many base cells it would come to if it
were split all the way down.
"""
return _np.array(self.discretization,
dtype=_np.int64) ** self.max_levels
@property
def _size(self):
"""Each block's extent in base cells, from its level.
Held as the level rather than the size: a level is one small number
against three, and the two cannot then disagree.
"""
ratio = _np.array(self.discretization, dtype=_np.int64)
return self._coarse_size[None, :] \
// ratio[None, :] ** self._level[:, None]
def _ancestor(self, level):
"""Each block's ancestor at `level`, as one integer per block.
The flat lattice index of the ancestor's lower corner, which names it
among the blocks of that level. A block at `level` is its own
ancestor, and one coarser than `level` has none there and reads -1.
"""
ratio = _np.array(self.discretization, dtype=_np.int64)
size = self._coarse_size // ratio ** int(level)
corner = (self._origin // size[None, :]) * size[None, :]
shape = self.lattice_shape
key = (corner[:, 0] * shape[1] + corner[:, 1]) * shape[2] \
+ corner[:, 2]
return _np.where(self._level >= int(level), key, -1)
@property
def block_size(self):
return self._size * self.base_step
@property
def block_volume(self):
return _np.prod(self.block_size, axis=1)
@property
def level(self):
return self._level
@property
def rows_per_location(self):
return int(_np.prod(self.discretization))
[docs]
def is_full(self) -> bool:
"""Whether the blocks tile their box exactly -- no gap, no overlap.
Volume alone cannot answer: a gap and an overlap of the same size
cancel. So the base cells are counted, in an array the shape of the
lattice, which costs a byte per base cell (29 MB for a 30-million-cell
model). Construction is full and both `split` and `group` preserve it,
so this is for checking something that came from elsewhere rather than
for routine use.
"""
covered = int(_np.prod(self._size, axis=1).sum())
if covered != int(_np.prod(self.lattice_shape)):
return False
seen = _np.zeros(tuple(self.lattice_shape), dtype=_np.uint8)
for origin, size in zip(self._origin, self._size):
seen[origin[0]:origin[0] + size[0],
origin[1]:origin[1] + size[1],
origin[2]:origin[2] + size[2]] += 1
return bool(seen.min() == 1 and seen.max() == 1)
_NO_SUBSET = (
"a %s cannot be subsetted: it is structurally complete by design, and "
"removing blocks would leave a group averaging over children that are "
"no longer there -- the mass-conservation error the lattice exists to "
"prevent. Exclude ground by value instead, with a metadata column "
"(`add_metadata`); `predict(..., where=...)` visits part of a model "
"without making a smaller one. To hand the blocks to something else, "
"`as_data_frame` carries a size per row, `as_pyvista` writes them as "
"hexahedra, and `to_zarr` round-trips the lattice whole.")
def __getitem__(self, item):
# PointData's would hand back a plain PointData, silently: no origin,
# no level, no size, so the size columns vanish, `as_pyvista` draws
# points rather than blocks and `grade_tonnage` refuses. Better to say
# so here than to return something that looks usable.
raise TypeError(self._NO_SUBSET % type(self).__name__)
[docs]
def subset_region(self, min_val, max_val,
include_min=None, include_max=None):
raise TypeError(self._NO_SUBSET % type(self).__name__)
# ------------------------------------------------------------------ #
# refinement
# ------------------------------------------------------------------ #
[docs]
def split(self, mask: _types.ArrayLike, carry: bool = True,
labels: _types.Labels = ("X", "Y", "Z")
) -> "BlockSet3D":
"""A new block set with each marked block cut into its own sub-blocks.
A block becomes `prod(discretization)` children, one per sub-block and
in the same order, so the values a coarse prediction already holds for
those sub-blocks describe the blocks this makes.
A block that was **not** split keeps what was predicted for it: it is
the same block on the same support, so its value is still the right
answer and arriving at it again would be work for nothing. The
children start missing, and `unpredicted()` says which they are, so
`predict(..., where=...)` visits only them. A parent's value is never
handed down -- that would manufacture children agreeing exactly, which
is the one thing refining is meant to find out rather than assume.
Parameters
----------
mask : array-like
One boolean per block, or the indices of the blocks to split.
carry : bool
Whether to bring the variables and metadata across. `False` gives
bare geometry, for building a mesh to predict onto from scratch.
"""
mask = _np.asarray(mask)
if mask.dtype != bool:
index = _np.zeros(self.n_data, dtype=bool)
index[mask] = True
mask = index
if mask.shape != (self.n_data,):
raise ValueError(
"the mask needs one value per block: got %d for %d"
% (mask.size, self.n_data))
finest = self._level >= self.max_levels
if _np.any(mask & finest):
raise ValueError(
"%d block(s) are already at the finest level the lattice "
"allows (max_levels=%d); build the set with more levels to "
"go further"
% (int(_np.count_nonzero(mask & finest)), self.max_levels))
# one child per sub-block, in the sub-blocks' own order
corner = _gmt.sub_block_index(self.discretization)
child_size = self._size[mask] // _np.array(self.discretization)
children = (self._origin[mask][:, None, :]
+ corner[None, :, :] * child_size[:, None, :]
).reshape(-1, 3)
child_level = _np.repeat(self._level[mask] + 1, len(corner))
new = _copy.copy(self)
new._origin = _np.concatenate([self._origin[~mask], children])
new._level = _np.concatenate([self._level[~mask], child_level])
# the blocks this made, so that predicting can be told to visit only
# them without first being told what was predicted before
new._fresh = _np.concatenate([
_np.zeros(int(_np.count_nonzero(~mask)), dtype=bool),
_np.ones(len(children), dtype=bool)])
new.variables = {}
new.metadata = {}
new._init_coordinates(new._centres(), list(labels))
# `_init_coordinates` has just measured the box off the coordinates,
# which are block *centres*: that box is half a block short on every
# side, and short by a different amount after every split, since
# refining an edge block moves its centre nearer the boundary. Two
# sets over the same ground would then disagree about their extent.
# The real box is the lattice, which no refinement changes.
new._bounding_box = new._lattice_bounding_box()
if carry:
# A block that was not split is the same block on the same
# support, so what was predicted for it still stands; only the new
# children have nothing yet. Metadata goes to the children as
# well, and by inheritance rather than as missing: it describes the
# ground, and a child occupies its parent's ground.
parent = _np.concatenate(
[_np.flatnonzero(~mask),
_np.repeat(_np.flatnonzero(mask), len(corner))])
for name, variable in self.variables.items():
new.variables[name] = variable.carry_to(
new, ~mask, len(children))
for name, column in self.metadata.items():
values = _np.asarray(column.values)[parent]
fresh = _Attribute(new, values, dtype=values.dtype)
fresh.labels = column.labels
new.metadata[name] = fresh
return new
[docs]
def group(self, mask: _types.ArrayLike, carry: bool = True,
labels: _types.Labels = ("X", "Y", "Z")
) -> "BlockSet3D":
"""The inverse of `split`: whole families of children, back into the
parent they came from.
A block is grouped with its siblings, so the mask must name **every**
child of a parent or none of them. A partial family would average over
children that are not there and mis-weight the parent -- the
mass-conservation error the lattice exists to prevent -- so it is
refused rather than approximated. That check is what makes conversion
between supports two-directional: `group` undoes `split` exactly, and
the blocks tile their box afterwards as they did before.
A block that was **not** grouped keeps what it holds, exactly as in
`split` and for the same reason: it is the same block on the same
support, so its value is still the right answer for it. The parents
are the ones on new ground, and they come back missing --
`unpredicted()` names them and `predict(..., where=...)` fills them.
A parent's value is never averaged from its children. Coarsening is a
change of support and almost nothing survives it: a parent's spread is
not its children's, its within-block dispersion is larger by exactly
what the grouping absorbed, and a category has no mean. The one thing
that would come across exactly is a realization, and a variable is
more than its realizations. Metadata does come across, from the first
child -- it describes the ground rather than the model, and where the
children disagree about it there is no right answer to be had.
Parameters
----------
mask : array-like
One boolean per block, or the indices of the blocks to group.
carry : bool
Whether to bring across what the blocks that were not grouped
hold. `False` gives bare geometry.
labels : list
Coordinate names.
"""
mask = _np.asarray(mask)
if mask.dtype != bool:
index = _np.zeros(self.n_data, dtype=bool)
index[mask] = True
mask = index
if mask.shape != (self.n_data,):
raise ValueError(
"the mask needs one value per block: got %d for %d"
% (mask.size, self.n_data))
coarsest = self._level < 1
if _np.any(mask & coarsest):
raise ValueError(
"%d block(s) are already at the coarsest level the lattice "
"has and belong to no parent"
% int(_np.count_nonzero(mask & coarsest)))
ratio = _np.array(self.discretization, dtype=_np.int64)
family = self._size[mask] * ratio[None, :]
parent = (self._origin[mask] // family) * family
level = self._level[mask] - 1
# one key a parent, level included: two families of different sizes can
# share a lower corner, and they are not the same parent
shape = self.lattice_shape
key = (((parent[:, 0] * shape[1] + parent[:, 1]) * shape[2]
+ parent[:, 2]) * (self.max_levels + 1) + level)
_, first, counts = _np.unique(key, return_index=True,
return_counts=True)
whole = int(_np.prod(self.discretization))
if _np.any(counts != whole):
short = _np.flatnonzero(counts != whole)
raise ValueError(
"grouping takes whole families: %d parent(s) were named by "
"only some of their %d children, the first by %d. A partial "
"family averages over blocks that are not there, which is the "
"one thing the lattice exists to prevent"
% (len(short), whole, int(counts[short[0]])))
new = _copy.copy(self)
new._origin = _np.concatenate([self._origin[~mask], parent[first]])
new._level = _np.concatenate([self._level[~mask], level[first]])
# the blocks this made, on a support nothing has been predicted on
new._fresh = _np.concatenate([
_np.zeros(int(_np.count_nonzero(~mask)), dtype=bool),
_np.ones(len(first), dtype=bool)])
new.variables = {}
new.metadata = {}
new._init_coordinates(new._centres(), list(labels))
# the lattice is the box, not the spread of the centres -- see `split`
new._bounding_box = new._lattice_bounding_box()
if carry:
source = _np.concatenate(
[_np.flatnonzero(~mask), _np.flatnonzero(mask)[first]])
for name, variable in self.variables.items():
new.variables[name] = variable.carry_to(
new, ~mask, len(first))
for name, column in self.metadata.items():
values = _np.asarray(column.values)[source]
fresh = _Attribute(new, values, dtype=values.dtype)
fresh.labels = column.labels
new.metadata[name] = fresh
return new
[docs]
def as_blocks3d(self) -> "Blocks3D | RotatedBlocks3D":
"""This model at its coarsest level, as a regular block model.
One block per coarsest-level block, in `Blocks3D`'s own order and
with this model's discretization, so a block that was never split is
the same block on the same support and keeps every column exactly.
Where the model was refined, each coarse block is gathered from the
blocks inside it, weighted by volume:
- a column that is a mean over a block's sub-blocks -- `prediction`,
the latent moments, `noise_variance`, the shares below each
cut-off, a category's probability and entropy -- takes the mean of
its parts', which is that same mean over the whole block;
- every realization is averaged index by index, the parts'
realization `i` making the block's realization `i`;
- what is read off the realizations is read again: the quantiles,
the probabilities, the `dispersion` -- the parts' mean dispersion
plus the spread between them -- and the predicted category, the
winner of the averaged probabilities;
- anything else is missing, none of it following from the parts:
`divided`, the measurements, a binary variable's entropy.
A coarse block with any part that holds nothing -- a block `where=`
kept from being predicted, say -- holds nothing either, so a partial
answer never passes for a whole one. Metadata is gathered as
`aggregate` gathers it, by volume: a number averages, a coded column
keeps the label holding most of the block (none on a tie), and a
boolean one holds where it held throughout.
Returns
-------
Blocks3D or RotatedBlocks3D
A `RotatedBlocks3D` for a `RotatedBlockSet3D`, with the same
angles and turning about the same point.
See Also
--------
group : the coarsening that leaves a parent missing instead, since
the change of support is otherwise the model's to re-predict.
Notes
-----
A derived variable is averaged like any other, which suits a
quantity that adds up over a volume. One that does not needs `derive`
again on the result, where its function sees the coarse blocks.
"""
coarse = self._coarse_size
n = self.lattice_shape // coarse
index = self._origin // coarse[None, :]
# the regular model counts its blocks x fastest (`Grid3D._generate`)
cell = index[:, 0] + n[0] * (index[:, 1] + n[1] * index[:, 2])
grouping = _Grouping(cell, _np.prod(self._size, axis=1),
int(_np.prod(coarse)), int(_np.prod(n)))
blocks = self._coarsest(n, coarse * self.base_step)
for name, variable in self.variables.items():
fresh = variable.__class__.from_variable(blocks, variable)
variable._copy_attrs_into(fresh)
variable._coarsen_into(fresh, grouping)
blocks.variables[name] = fresh
for name, column in self.metadata.items():
values = grouping.gathered(column)
fresh = _Attribute(blocks, values, dtype=values.dtype)
fresh.labels = column.labels
blocks.metadata[name] = fresh
return blocks
def _coarsest(self, n, step) -> "Blocks3D | RotatedBlocks3D":
"""The empty regular model `as_blocks3d` fills: `n` blocks of `step`
over this model's box."""
return Blocks3D(start=self.box_corner + step / 2, n=n, step=step,
# added to the grid's arguments by `_blockdata`
discretization=list( # type: ignore[call-arg]
self.discretization),
labels=list(self.coordinate_labels))
# ------------------------------------------------------------------ #
# geoh5 interchange
# ------------------------------------------------------------------ #
[docs]
def to_geoh5(self, workspace: "_types.PathLike | _GeoH5Workspace",
name: str = "Blocks", include: str = "**",
simulations: "bool | int | Sequence[int]" = False,
replace: bool = True,
folder: "str | None" = None) -> None:
"""
Writes this block model into a geoh5 workspace, as an Octree.
The lattice maps one to one when `discretization` is `(2, 2, 2)`
— the default — since a geoh5 octree subdivides strictly two by
two by two; any other discretization is refused rather than
resampled. A workspace is what Geoscience ANALYST — a free
viewer for the format — opens as **one** project, so a model's
pieces belong in one file, as in `PointData.to_geoh5`; several
exports in a row go fastest through an open
`geoml.data.geoh5.Workspace`. Needs the `geoh5py` package:
`pip install geoml[geoh5]`.
Parameters
----------
workspace
Path of the workspace to write into, or an open
`geoml.data.geoh5.Workspace`.
name
The name the object gets in the workspace.
include
A path pattern selecting the columns to carry, as in
`as_pyvista`.
simulations
Which simulations to include: `False` for none (the default),
`True` for all, an `int` for the first n, or a sequence of
indices.
replace
Whether an existing Octree of this name, in this folder,
makes way — what a re-run export script means. `False` keeps
both, and reading that name back then requires saying which.
folder
Where the object sits in ANALYST's project tree, as a path —
each segment a group, created when it does not exist and
reused when it does. `None` is the root.
Raises
------
ValueError
If the discretization is not `(2, 2, 2)`, or the model
carries a dip or rake — geoh5 rotates about the vertical
axis alone.
"""
import geoml.data.geoh5 as _geoh5io
_geoh5io.write_blocks(self, workspace, name, include, simulations,
replace, folder)
[docs]
@classmethod
def from_geoh5(cls, workspace: "_types.PathLike | _GeoH5Workspace",
name: "str | None" = None) -> "BlockSet3D":
"""
Reads an Octree from a geoh5 workspace, as a block model.
The octree cells become blocks on a `(2, 2, 2)` lattice, in the
file's own row order, with float cell data as continuous
variables and referenced data as categorical ones. A rotated
octree comes back as a `RotatedBlockSet3D`; an octree written
with its vertical axis running downward is normalized on the way
in. `name` says which Octree to read when the workspace holds
more than one; `geoml.data.geoh5.contents` lists what there is
to name. Needs the `geoh5py` package: `pip install geoml[geoh5]`.
A foreign file is not ours to trust, and the always-full design
does not bend for it: cells must sit aligned to their own
power-of-two size and free of overlaps, and whatever the cells
do not cover inside their box is filled with unvalued blocks, a
metadata column named ``imported`` marking the difference — a
vendor octree usually carries only the domain of interest.
Parameters
----------
workspace
Path of the workspace to read, or an open `geoml.data.geoh5.Workspace`.
name
The Octree object to read, when the file holds more than one.
Returns
-------
blocks : BlockSet3D or RotatedBlockSet3D
Whichever the file's rotation calls for.
Raises
------
ValueError
If the workspace holds no such Octree, or its cells are not
an octree this lattice can hold (misaligned, overlapping, or
reflected under a rotation).
"""
import geoml.data.geoh5 as _geoh5io
raw = _geoh5io.read_blocks(workspace, name)
origin, size = raw["origin"], raw["size"]
safe = _np.where(size > 0, size, 1)
aligned = (size > 0) & ((size & (size - 1)) == 0) \
& (origin % safe[:, None] == 0).all(axis=1)
if not aligned.all():
raise ValueError(
"%d cell(s) of %s are not aligned power-of-two octree "
"cells, which is not a lattice this class can hold"
% (int(_np.count_nonzero(~aligned)), raw["filename"]))
offending = _gmt.dyadic_overlaps(origin, size)
if len(offending):
raise ValueError(
"%d cell(s) of %s sit inside larger cells of the same "
"file; an overlapping model has no one value per block"
% (len(offending), raw["filename"]))
# capacity from the writer's own record where there is one -- the
# cells alone cannot say it once every coarse block has been split
max_levels = int(raw["lattice"].get(
"max_levels", int(_np.log2(size.max()))))
max_levels = max(max_levels, int(_np.log2(size.max())))
root = 1 << max_levels
# the box on root multiples, in the file's own frame -- shifting
# by anything else would break the alignment just checked
low = (origin.min(axis=0) // root) * root
high = -(-(origin + size[:, None]).max(axis=0) // root) * root
origin = origin - low
shape = high - low
gap_origin, gap_size = _gmt.dyadic_complement(
origin, size, shape, root)
origin = _np.concatenate([origin, gap_origin])
size = _np.concatenate([size, gap_size])
level = max_levels - _np.log2(size).astype(_np.int64)
# the constructor takes the first coarse block's world *centre*;
# the corner offset turns with the rotation, the way every
# coordinate leaves the lattice frame
step = raw["step"]
azimuth = -raw["rotation"]
matrix = _gmt.rotation_matrix(azimuth, 0.0, 0.0)
centre_local = (low + root / 2.0) * step
start = centre_local @ matrix + raw["corner"]
if raw["rotation"] == 0.0:
blocks = BlockSet3D(start, shape // root, root * step,
max_levels=max_levels)
else:
blocks = RotatedBlockSet3D(start, shape // root, root * step,
azimuth=azimuth,
max_levels=max_levels)
# the mixed levels go in the way `split` writes its own
blocks._origin = origin
blocks._level = level
blocks._fresh = _np.ones(len(origin), dtype=bool)
blocks._init_coordinates(blocks._centres(), ["X", "Y", "Z"])
blocks._bounding_box = blocks._lattice_bounding_box()
for column, values in raw["floats"]:
blocks.add_continuous_variable(
column, _np.concatenate([values,
_np.full(len(gap_size), _np.nan)]))
for column, labels, values in raw["coded"]:
values = _np.asarray(values, dtype=object)
values[values == ""] = None
values = _np.concatenate(
[values, _np.full(len(gap_size), None, dtype=object)])
blocks.add_categorical_variable(column, labels=labels,
measurements=values)
if len(gap_size):
blocks.add_metadata(
"imported",
_np.concatenate([_np.ones(len(raw["origin"]), dtype=bool),
_np.zeros(len(gap_size), dtype=bool)]))
return blocks
[docs]
def block_shares(self, split_on: "str | Sequence[str] | None" = None
) -> "dict[str, _np.ndarray]":
"""How often each decision cuts a block in two.
A continuous variable contributes one column per cut-off it declares,
a categorical one per category, and both mean the same thing: the
share of realizations in which this block's sub-blocks fall on both
sides of something that matters. A grade is judged against the grades
someone declared; an indicator against zero, that being where one
category stops winning and its rival starts.
Returns a dict of name -> `(n_data,)` array, empty where nothing
declared a decision to make.
"""
chosen = list(self.variables) if split_on is None else \
[str(name) for name in _np.atleast_1d(split_on)]
shares = {}
for name in chosen:
if name not in self.variables:
raise ValueError(
"no variable named %r to split on; found %s"
% (name, ", ".join(sorted(self.variables)) or "none"))
# each variable reports its own -- see `_Variable.split_shares`
for key, values in self.variables[name].split_shares().items():
shares[("%s %s" % (name, key)).strip()] = values
return shares
[docs]
def needs_splitting(self, split_on: "str | Sequence[str] | None" = None,
tolerance: "float | None" = None) -> _np.ndarray:
"""Which blocks a decision's surface runs through, and so are worth
cutting.
A block is marked where the prediction at its sub-blocks falls on both
sides of a cut-off, or of a category's boundary: the surface a contour
of the prediction draws passes through it, and cutting the block lets
that surface bend there.
Note what this does *not* mark: a block the model is merely unsure
about. The prediction is the mean over the realizations, and where the
data do not reach, each realization crosses the cut-off somewhere of
its own while their mean crosses it nowhere -- cutting there would buy
blocks and no surface.
The test is over every decision at once and any one is enough, which
is the cautious way round on purpose: a block worth splitting for one
variable is worth splitting whatever the others say. Name `split_on`
to narrow it -- letting every element of a polymetallic deposit vote
marks most of the model and gives back the saving.
Blocks already at the finest level are never marked, there being
nothing to cut them into.
Parameters
----------
split_on : str or list, optional
Which variables get a say. All of them by default.
tolerance : float, optional
Deprecated, and without effect: a block is divided or it is not.
It was the share of realizations that had to find a block divided,
and will be removed.
"""
if tolerance is not None:
_warnings.warn(
"`tolerance` has no effect since 0.8.5 -- a block is divided "
"where the prediction's sub-blocks straddle a cut-off, which "
"is yes or no -- and will be removed", FutureWarning,
stacklevel=2)
shares = self.block_shares(split_on)
mask = _np.zeros(self.n_data, dtype=bool)
for values in shares.values():
# 0 or 1; a store written before 0.8.5 holds a share of the
# realizations, read here as a majority
mask |= _np.asarray(values, dtype=float) > 0.5
return mask & (self._level < self.max_levels)
def _by_level(self):
"""Where the blocks of each level sit, as a sorted table of origins.
Sorted rather than rasterized on purpose: naming the block covering a
point is then a `searchsorted`, where painting the base lattice would
cost a cell for every one the model exists to avoid carrying.
"""
shape = self.lattice_shape
tables = []
for k in range(self.max_levels + 1):
rows = _np.flatnonzero(self._level == k)
here = self._origin[rows]
key = ((here[:, 0] * shape[1] + here[:, 1]) * shape[2]
+ here[:, 2])
order = _np.argsort(key, kind="stable")
tables.append((key[order], rows[order]))
return tables
[docs]
def unbalanced(self, gap: int = 1) -> _np.ndarray:
"""Blocks with a neighbour more than `gap` levels finer than they are.
A block whose own sub-blocks agree is never marked by
`needs_splitting`, and rightly so -- cutting it would not change the
answer it gives. But the field can still turn sharply inside it, and
nothing in the block itself says so. That a *neighbour* was cut twice
while this block was not cut at all is the evidence, and it lives
outside the block.
It matters for what is drawn rather than for what is decided. A
contour reads a block through its eight corners, so a coarse block
beside much finer ones is a crude straight guess across a long span
laid right where the surface runs. Levelling the jump measured three
times closer to the true surface for 35% more blocks, where refining a
whole level deeper without it bought almost nothing for 2.6 times as
many -- deeper refinement widens the jumps as fast as it narrows the
blocks. `models.refine` therefore cuts these as it goes.
Blocks already at the finest level are never marked, there being
nothing to cut them into.
Parameters
----------
gap : int
How many levels of difference to tolerate. One is the usual 2:1
balance: a block may meet blocks one level finer, not two.
"""
ratio = _np.array(self.discretization, dtype=_np.int64)
shape = self.lattice_shape
tables = self._by_level()
size = self._size
marked = _np.zeros(self.n_data, dtype=bool)
# asked from the fine side: a fine block steps one cell past each of
# its faces and names the coarse block covering that point, which is
# exact where asking the coarse block what lies beyond its face is not
# -- a cell of its own size holds blocks that do not touch it
for fine in range(int(gap) + 1, self.max_levels + 1):
rows = _np.flatnonzero(self._level == fine)
if len(rows) == 0:
continue
origin, extent = self._origin[rows], size[rows]
for coarse in range(fine - int(gap)):
key_of, owner = tables[coarse]
if len(owner) == 0:
continue
cell = self._coarse_size // ratio ** coarse
for axis in range(3):
for side in (-1, 1):
beyond = origin.copy()
beyond[:, axis] += (extent[:, axis] if side > 0
else -1)
inside = ((beyond[:, axis] >= 0)
& (beyond[:, axis] < shape[axis]))
owning = (beyond // cell) * cell
key = ((owning[:, 0] * shape[1] + owning[:, 1])
* shape[2] + owning[:, 2])
at = _np.minimum(_np.searchsorted(key_of, key),
len(key_of) - 1)
hit = inside & (key_of[at] == key)
marked[owner[at[hit]]] = True
return marked & (self._level < self.max_levels)
def _lattice_bounding_box(self):
"""The lattice is the box, not the spread of the centres -- see
`split`. A rotated subclass answers with the rotated corners."""
return BoundingBox(
self.box_corner,
self.box_corner + self.lattice_shape * self.base_step)
[docs]
@classmethod
def from_data(cls, data, step, margin=0.1, decimals=0,
discretization=(2, 2, 2), max_levels=3):
"""
A block model covering another object's bounding box.
As `Grid3D.from_data`, counting blocks rather than nodes: the
margined box's lower *corner* is floored to `decimals`, and enough
blocks follow to cover the far side, so the corner is round and the
margin never shrinks.
Parameters
----------
data
Any spatial object, drillholes included.
step
The coarse block size, one number or one per direction.
margin : float or array
A fraction of the data's extent; see `Grid3D.from_data`.
decimals : int
Decimals to floor the box corner to.
discretization, max_levels
As the constructor takes them.
"""
corner, n, step, labels = _cover_box(
data, step, margin, decimals, n_dim=3, cells=True)
return cls(start=corner + step / 2, n=n, step=step,
discretization=discretization, max_levels=max_levels,
labels=labels if labels else ("X", "Y", "Z"))
# ------------------------------------------------------------------ #
# sample data
# ------------------------------------------------------------------ #
[docs]
def index_data(self, data: "_SpatialData") -> _np.ndarray:
"""Which block each of `data`'s locations falls in.
One row index per location, `-1` for anything outside the box. Note
this is **not** what a grid's `index_data` returns -- a cell index per
axis -- because blocks of several sizes have no per-axis index to
return. Which block *is* the answer here.
The lattice makes it cheap: a location's base cell is arithmetic, and
the block covering that cell is the one whose origin is the cell's
ancestor at that block's own level, so one `searchsorted` per level
finds it and every location is settled within `max_levels + 1` of them.
"""
if data.n_dim != self.n_dim:
raise DimensionMismatchError(
"Data dimension mismatch. Expected dimension %d and found %d."
% (self.n_dim, data.n_dim))
coordinates = _np.asarray(data.coordinates, dtype=float)
cell = _np.floor(
(coordinates - self.box_corner) / self.base_step).astype(_np.int64)
shape = self.lattice_shape
found = _np.full(len(coordinates), -1, dtype=_np.int64)
left = _np.flatnonzero(
_np.all((cell >= 0) & (cell < shape), axis=1))
ratio = _np.array(self.discretization, dtype=_np.int64)
tables = self._by_level()
for k in range(self.max_levels + 1):
if len(left) == 0:
break
key_of, owner = tables[k]
if len(owner) == 0:
continue
size = self._coarse_size // ratio ** k
ancestor = (cell[left] // size) * size
key = ((ancestor[:, 0] * shape[1] + ancestor[:, 1]) * shape[2]
+ ancestor[:, 2])
at = _np.minimum(_np.searchsorted(key_of, key), len(key_of) - 1)
hit = key_of[at] == key
found[left[hit]] = owner[at[hit]]
left = left[~hit]
return found
def _cell_of(self, data):
# which block, already this object's own row index
return self.index_data(data)
[docs]
def aggregate(self, data, variables=None, metadata=True):
"""Carries another object's measurements onto the blocks holding them.
As `Grid3D.aggregate` -- one method, the operation following each
variable's kind -- over blocks of several sizes.
"""
_aggregate_onto(self, data, variables, metadata)
return self
[docs]
def assign_from_surface(self, surface, name, labels=("above", "below"),
fraction=None, uncovered=None):
"""
As `Blocks3D.assign_from_surface`, over blocks of several sizes.
The flag in `name` follows the block centre; naming a `fraction`
column measures the share of each block below the sheet over the
sub-blocks `discretization` defines, scaled to each block's own size.
`crossed_by` asks the same question and answers with the blocks worth
cutting.
"""
_blocks_from_surface(self, PointData.assign_from_surface, surface,
name, labels, fraction, uncovered)
[docs]
def assign_from_solid(self, solid, name, labels=("outside", "inside"),
fraction=None):
"""
As `Blocks3D.assign_from_solid`, over blocks of several sizes.
`fraction` holds the share of each block's volume inside the body, and
`crossed_by` turns the same test into the blocks worth cutting.
"""
_blocks_from_solid(self, PointData.assign_from_solid, solid, name,
labels, fraction)
[docs]
def crossed_by(self, mesh: "Mesh3D") -> _np.ndarray:
"""Which blocks a mesh passes through, and so which are worth cutting.
A block is crossed when its sub-blocks fall on **both** sides of the
mesh -- the question `needs_splitting` asks of a cut-off, asked of
geometry instead. A topography, a vein wall, a lease boundary: a block
the surface runs through holds two answers whatever is predicted into
it, and no amount of prediction will separate them.
One entirely above or entirely below is left alone however close it
lies, which is what keeps this from refining a whole domain. Blocks
already at the finest level are never marked, as elsewhere.
Hand it to `split`, or give the mesh to `models.refine`, which unions
this with the other two criteria and repeats until nothing is left to
cut::
blocks = blocks.split(blocks.crossed_by(topography))
A sheet that covers only part of the model counts a sub-block past its
edge as not below, the way `fraction` does, so a block straddling the
sheet's own boundary reads as crossed. That is usually wanted -- the
edge is a real feature of the ground being described -- but it is why
a sheet should reach across the model when it is not.
Parameters
----------
mesh : Surface3D or Solid3D
The sheet or closed body to test against. A `Mesh3D` that is
neither has no sides and is refused.
Returns
-------
array
One boolean per block.
"""
splittable = _np.flatnonzero(self._level < self.max_levels)
if len(splittable) == 0:
return _np.zeros(self.n_data, dtype=bool)
shares = _sub_block_shares(self, _mesh_test(mesh), rows=splittable)
# nan everywhere else, and a comparison against it is False, so the
# blocks that were never asked about answer no on their own
with _np.errstate(invalid="ignore"):
return (shares > 0) & (shares < 1)
[docs]
def unpredicted(self, variable: str | None = None) -> _np.ndarray:
"""Which blocks have nothing in them yet.
What `split` leaves behind, and what a cancelled prediction did not
reach: hand it to `predict(..., where=...)` and only those blocks
are visited. Read off the missing values, as on every container,
once the set holds a variable; before anything has been predicted
into it, the blocks the last split made. Naming a variable reads
that one's missing values alone.
"""
if variable is not None:
return self.variables[variable].unpredicted()
# A split leaves its children missing in every variable the set
# holds, so the values say it all; the flags are what is left to say
# it with no variable at all. Their union, which this was in 0.7.0,
# kept every block of a new set unpredicted after predicting it in
# full, the flags never being cleared.
if any(v._PREDICTED_MARKER is not None
for v in self.variables.values()):
return _SpatialData.unpredicted(self)
return _np.array(self._fresh, dtype=bool)
# ------------------------------------------------------------------ #
# what the model asks for
# ------------------------------------------------------------------ #
[docs]
def get_batched_coordinates(self, index=None):
if index is None:
index = _np.arange(self._n_data)
centres = _np.asarray(self.coordinates[index])
size = self.block_size[index]
# block-major, one temporary rather than one per block: block b owns
# rows b*k to (b+1)*k, which is what `_aggregate` averages back. The
# only difference from a fixed-size block model is that the sub-grid
# is scaled by each block's own size.
coords = centres[:, None, :] + self._sub_grid[None, :, :] \
* size[:, None, :]
coords = coords.reshape(-1, centres.shape[1])
splits = None if self.rows_per_location == 1 else len(centres)
return coords, splits
def _batch_rows(self, index=None):
n, _ = _PointBased._batch_rows(self, index)
k = self.rows_per_location
return n * k, (None if k == 1 else n)
[docs]
def as_data_frame(self, metadata: bool = True,
**kwargs) -> _pd.DataFrame:
df = super().as_data_frame(metadata=metadata, **kwargs)
for i, label in enumerate(self.coordinate_labels):
df["_" + label] = self.block_size[:, i]
return df
# ------------------------------------------------------------------ #
# export
# ------------------------------------------------------------------ #
[docs]
def as_pyvista(self, simulations=False, include="**"):
"""
Converts this object to a pyvista one, carrying its variables.
One hexahedron per block, written out corner by corner, rather than
the `ImageData` a regular block model exports: implicit geometry can
only say one spacing, and the point here is that there is more than
one. The cells are welded, so blocks that touch share the corners they
meet at -- which is what lets anything be contoured across them.
Parameters
----------
simulations
Which simulations to include: `False` for none (the default, since
each one is a full-length array in the exported object), `True` for
all of them, an `int` for the first n, or a sequence of indices.
"""
mesh = self._hex_mesh(self._origin, self._size, self.base_step)
return self._finish_pyvista(mesh, "cells", simulations, include)
def _to_world(self, coordinates):
"""An unrotated set's lattice frame is the world frame; the rotated
subclass overrides this with its rotation, and the contour's
vertices go through it exactly as `_hex_mesh`'s points do."""
return coordinates
def _hex_mesh(self, origin, size, step):
"""One welded hexahedron per block, on any lattice of blocks -- the
object's own, or the finer one a contour is drawn on. `origin` and
`size` count that lattice's own cells, `step` says how big one is."""
low = self.box_corner + origin * step
size = size * step
points = (low[:, None, :]
+ _gmt.HEX_CORNERS[None, :, :] * size[:, None, :]).reshape(-1, 3)
connectivity = _np.arange(len(points)).reshape(-1, 8)
cells = _np.hstack(
[_np.full((len(connectivity), 1), 8), connectivity]).ravel()
kinds = _np.full(len(connectivity), _pv.CellType.HEXAHEDRON,
dtype=_np.uint8)
# Welded, and it matters: written corner by corner every cell owns its
# own eight, so neighbours share nothing, and anything reading values
# at a corner sees one block there instead of the several that meet.
# A contour over an unwelded mesh comes back empty for that reason.
return _pv.UnstructuredGrid(cells, kinds, points).clean(
tolerance=1e-9, produce_merge_map=False)
@staticmethod
def _shared_corners(origin, size, shape):
"""Each block's eight corners, as indices into the distinct lattice
points -- the welding of `_hex_mesh`, done in integers."""
return _lattice_corners(origin, size, _np.asarray(shape) + 1)[1]
def _base_corners(self, unit):
"""`_shared_corners` of the model's own blocks on a lattice `unit`
times finer, kept. It depends on the lattice alone, and every contour
starts from it: a mesh set's at every cut-off and realization, whose
forked workers read the parent's copy. The lattice is set when a set
is made and never changed in place -- `split` and `group` make new
sets -- so the arrays it was read from say whether it still holds."""
unit = _np.asarray(unit, dtype=_np.int64)
fingerprint = (id(self._origin), id(self._level),
tuple(int(u) for u in unit))
cached = getattr(self, "_corner_cache", None)
if cached is None or cached[0] != fingerprint:
cached = (fingerprint, self._shared_corners(
self._origin * unit, self._size * unit,
self.lattice_shape * unit))
self._corner_cache = cached
return cached[1]
@staticmethod
def _at_corners(corners, values):
"""The mean over the blocks meeting at each corner, read back per
block -- what `cell_data_to_point_data` puts there, and so what decides
where the surface runs. A block holding no value contributes nothing,
rather than carrying its absence into every corner it touches.
"""
return BlockSet3D._corner_means(corners, values)[corners]
@staticmethod
def _corner_means(corners, values):
"""`_at_corners` once per corner point rather than per block."""
known = _np.isfinite(values)
length = int(corners.max()) + 1
# bincount rather than `np.add.at`: the same sums, several times
# faster -- add.at goes through the buffered-ufunc machinery
total = _np.bincount(corners[known].ravel(),
weights=_np.repeat(values[known], 8),
minlength=length)
count = _np.bincount(corners[known].ravel(), minlength=length)
return _np.divide(total, count, out=_np.full(length, _np.nan),
where=count > 0)
def _cut_to_contour(self, values, value, margin=1, supersample=1,
fields=None):
"""The blocks a surface runs through, cut small -- in the mesh handed
to VTK, not in the model.
A hexahedron is contoured from its own eight corners, so where a
coarse block meets finer ones it cannot see the corners they place in
the middle of the shared face. The two sides then draw different
curves there and the surface tears along every interface it crosses.
Cutting the surface's neighbourhood to one size puts the interfaces
out of its way, and one ring of margin covers the surface shifting as
the values are read at the finer size.
Cutting `supersample` levels past the model's own finest also smooths
it. What VTK contours is a trilinear reading of the corner averages,
continuous but with a crease at every face it crosses, and on a field
with structure at a few blocks those creases are what looks blocky.
Averaging onto corners again at each finer level composes into a
rounder reconstruction, so the surface both turns less sharply and
sits closer to the level set it is meant to be: on a rough test field,
one level took the mean angle between neighbouring triangles from 9.1
to 5.7 degrees and halved the distance from the true surface, matching
a model of 3.4 times as many blocks that had to be predicted.
Nothing is predicted, at any level. A child takes its parent's value
plus what the corners say about the shape running across it, and that
correction averages to zero over the children, so a block's estimate
is exactly the mean of the children standing in for it.
With `fields`, `values` holds several columns, each cut the same
way, and the surfaces are those of the fields `fields` reads off the
columns' corner means (see `_contour_fields`); a block is cut where
any of them runs through it.
Returns
-------
origin, size, step, values
The finer lattice: origins and sizes counted in its own cells,
and how big one of those is. `(None, None, None, None)` if nothing
was cut and the blocks as they stand will do.
"""
ratio = _np.array(self.discretization, dtype=_np.int64)
offset = _gmt.sub_block_index(self.discretization)
weights = _gmt.trilinear_weights(self.discretization)
# a lattice `supersample` levels finer than the model's own, so that
# cutting below the base cell is still whole numbers of a cell
unit = ratio ** int(supersample)
origin = self._origin * unit
size = self._size * unit
step = self.base_step / unit
shape = self.lattice_shape * unit
values = _np.asarray(values, dtype=float)
columns = values.shape[1] if values.ndim == 2 else 1
cut = False
for level in range(self.max_levels + int(supersample)):
corners = self._base_corners(unit) if level == 0 \
else self._shared_corners(origin, size, shape)
means = _np.stack(
[self._corner_means(corners, values) if values.ndim == 1
else self._corner_means(corners, values[:, k])
for k in range(columns)], axis=1)
decided = means if fields is None else _np.asarray(fields(means))
marked = _np.zeros(len(origin), dtype=bool)
for field in decided.T:
at_corner = field[corners]
# a corner with no value must not decide anything either
# way -- `fmin` and `fmax` pass over a NaN, and an all-NaN
# row compares false, without an n x 8 copy each
marked |= (_np.fmin.reduce(at_corner, axis=1) < value) \
& (_np.fmax.reduce(at_corner, axis=1) > value)
del at_corner
near = _gmt.grow(corners, marked, margin)
near &= _np.all(size >= ratio, axis=1)
# a block holding no value has nothing to hand down: whole, the
# paint reads it across from its corners where valued blocks
# have them all, and cut, its middle would be corners no valued
# block reaches
near &= _np.isfinite(values) if values.ndim == 1 \
else _np.all(_np.isfinite(values), axis=1)
if not _np.any(near):
break
cut = True
child_size = size[near] // ratio
children = (origin[near][:, None, :]
+ offset[None, :, :] * child_size[:, None, :]
).reshape(-1, 3)
shaped = []
for k in range(columns):
column = values[near] if values.ndim == 1 else values[near, k]
at_corner = means[:, k][corners[near]]
parts = at_corner @ weights.T
parts += (column - at_corner.mean(axis=1))[:, None]
shaped.append(parts.ravel())
# this level's table goes before the next one is built, rather
# than the two of them standing side by side at the peak
del corners, means, decided, marked, at_corner, parts
origin = _np.concatenate([origin[~near], children])
size = _np.concatenate(
[size[~near], _np.repeat(child_size, len(offset), axis=0)])
values = _np.concatenate(
[values[~near], shaped[0] if values.ndim == 1
else _np.stack(shaped, axis=1)])
if not cut:
return None, None, None, None
return origin, size, step, values
[docs]
def get_contour(self, path, value, supersample=0, simplify=None,
close=False):
"""
Isosurface through blocks of more than one size.
The field is painted onto slabs of a regular grid and contoured by
flying edges (`_painted_contour`) -- the model itself never carries
every block at the finest size, but a transient slab of the field
can, and contouring it is what an unstructured mesh of welded
hexahedra used to be built for, at most of the call's cost. The
answer is the one a model carried at the finest size throughout
would have given.
Values live on the cells and an isosurface needs them on the
corners, so what is painted are the corners -- the blocks meeting
at each one decide where the surface passes, one vote per block
whatever its size.
The blocks the surface runs through are cut to the finest size the
lattice allows before any of that, in the mesh handed to VTK and not
in the model: a coarse block cannot see the corners its finer
neighbours place in the middle of the face they share, so the two
sides draw different curves there and the surface tears along every
such interface it crosses. Cutting its neighbourhood to one size puts
the interfaces out of the way. Nothing is predicted -- a child reads
its parent's value and the shape its corners carry -- and the surface
comes back the one a model carried at the finest size throughout would
have given. See `_cut_to_contour`.
Blocks holding no value -- ground a `where=` left out, or a split not
yet predicted -- vote at no corner. One whose corners valued blocks
all share is read across from them; beyond that the ground has no
level to draw, and the surface stops at it, or with `close` closes
against it.
Parameters
----------
path : str
The column to contour, named the way the tree names it:
`"grade/prediction"` is that column, `"grade"` alone defaults to
the variable's prediction, and `"Elements/Zn"` reaches a
component (only the components of a composition hold a grade). A
single bare segment that is no variable of its own is searched
for anywhere in the tree, so `"Zn"` still finds `"Elements/Zn"`
as long as only one variable holds a `Zn`.
value : float
The value to draw the surface at.
supersample : int
How many levels past the model's own finest block to cut the mesh
to. Costs `prod(discretization)` times the cells per level, around
the surface only. On a *smooth* field one level buys a rounder
and genuinely closer surface -- what VTK reads between block
corners is trilinear, creasing at every face, and averaging onto
corners again at each finer level composes into a reconstruction
worth roughly predicting a model several times the size. On a
near-binary field -- an indicator contoured close to zero -- the
surface is corner-locked either way, and the extra lattice buys
almost nothing: measured on a real 908k-block model, one level
cost 3.5-5.8x the time and moved the volumes by 0.03-0.17%. The
default is 0 for that reason; pass 1 when contouring a smooth
grade at a cut-off well inside its range.
simplify : float, optional
A geometric error budget, in coordinate units: the surface is
simplified until pushing further would move it more than this
(see `Mesh3D.simplify`). Pairs naturally with `supersample`,
which buys accuracy in triangles this then spends back where
the surface is flat. None returns the full triangulation.
close : bool or str
Whether to close the surface where it runs out of the model, so
that what comes back is a body rather than a sheet with a hole
in the side -- `"above"` (or `True`) keeps the region where the
values exceed `value`, a grade shell, and `"below"` the region
under it, as on a grid. Done with a shell of ghost cells
carrying the boundary blocks' values on past the box, each its
partner's own size, and the box then cutting the body off, so
its faces, edges and corners on the box are the box's own.
Returns
-------
surf : Solid3D, Surface3D or Mesh3D
Whichever the geometry calls for, as `get_contour` on a grid.
"""
# the column comes straight from the container by its resolved path
# -- `as_pyvista` used to be built here just to read one column out
# of it, which welded every block's corners and exported every
# variable's every attribute first; on a model of any size that was
# most of the call, and all of it thrown away
path, column = _contour_column(self, path)
if not column._has_content():
raise ValueError("nothing under %r to contour" % str(path))
values = _np.asarray(column.values, dtype=float).ravel()
surface = self._contour_values(values, value, str(path),
supersample=supersample, close=close)
if surface is None:
raise ValueError(
"no surface at %g; %r runs from %g to %g"
% (float(value), str(path), float(_np.nanmin(values)),
float(_np.nanmax(values))))
if simplify is not None:
surface = surface.simplify(simplify)
surface.provenance = {
"source": str(path), "value": float(value),
"close": None if not close
else ("above" if close is True else str(close)),
"supersample": int(supersample),
"simplify": None if simplify is None else float(simplify)}
return surface
def _contour_values(self, values, value, label, supersample=0,
close=False, fallback=True):
"""`get_contour` on an array of block values rather than a column.
What a mesh set contours every realization through, the column
already read. None where the field never reaches `value`: an empty
answer is a fact about a realization rather than an error, and the
caller decides which it is. `fallback=False` hands back a painted
surface that is no body as it is, rather than contouring the welded
mesh again: a mesh set moves the level instead, which the fallback
never beat, and on a big model the welded mesh is the costliest
thing a contour can do (14 million cells, 17 GB and two minutes an
attempt, measured on the Tom v6 model).
"""
value = float(value)
origin, size, step, cell_values = self._cut_to_contour(
values, value, supersample=supersample)
if origin is None:
origin, size = self._origin, self._size
step, cell_values = self.base_step, values
cell_values = _np.asarray(cell_values)
mirrored = None
if close:
# `_closing_value` raises on anything that is not a side, and
# says which one is kept
background = _closing_value(close, values)
keep_above = background < value
origin, size, cell_values, inside = self._with_ghosts(
origin, size, step, cell_values)
mirrored = cell_values[inside:].copy()
else:
background = _np.nan
# One field, one route: the painted grid carries the welded mesh's
# own point field -- cell-equal corner means, trilinear interiors
# -- so flying edges draws the surface the welded hexahedra would,
# at a fraction of the cost (see `_painted_contour`; an earlier
# version painted cell values instead, whose volume-weighted
# interface corners pinched near-binary indicators, and those
# fields had to be routed around it). The classification stays the
# arbiter regardless: a closed-but-inconsistent Mesh3D is never
# what a level set means, and sends the result back through the
# welded mesh, which remains the last word on what the field says.
verts, faces = _painted_contour(origin, size, step, cell_values,
self.box_corner, value, label,
background,
margin=1 if close else 0)
if close and len(faces):
# The ghosts carry the field on past the box as it stands at
# the box -- each a copy of the block it mirrors -- so the
# surface crosses the box's faces rather than lying in them,
# and closes a cell past the ghosts, against the padding. The
# box then cuts it, exactly, faces, edges and corners alike.
# The cap used to be painted onto the faces instead, by ghosts
# holding each block's reflection about the level: every face
# point of a kept block then sat exactly at the level, which
# rounded the box's edges by half a boundary block (5% of an
# 80 m box three slabs fill) and folded the cap flat onto
# itself wherever the kept ground thinned against it. That
# remains the way a body Manifold refuses is closed.
clipped = _clip_to_box(verts, faces, *self._lattice_box())
if clipped is not None:
verts, faces = clipped
else:
cell_values[inside:] = _ghost_values(mirrored, value,
keep_above)
verts, faces = _painted_contour(origin, size, step,
cell_values,
self.box_corner, value,
label, background)
mirrored = None
if len(faces) == 0:
return None
surface = self._surface_of(verts, faces)
if type(surface) is Mesh3D and surface.n_data > 0 and fallback:
if mirrored is not None:
# the welded mesh has no padding to close against, so its
# ghosts put the cap on the faces the old way
cell_values[inside:] = _ghost_values(mirrored, value,
keep_above)
# VTK averages cells onto points with a NaN in the sum, so the
# blocks holding no value stay out of the welded mesh, and the
# corners read the valued blocks alone, as the paint reads them
valued = _np.isfinite(cell_values)
mesh = self._hex_mesh(_np.asarray(origin)[valued],
_np.asarray(size)[valued], step)
mesh.cell_data[label] = cell_values[valued]
welded = mesh.cell_data_to_point_data().contour(
[value], scalars=label)
if welded.n_cells == 0:
return None
welded = welded.triangulate()
verts = _np.asarray(welded.points, dtype=float)
faces = _np.asarray(welded.faces).reshape(-1, 4)[:, 1:]
# the level set pinching at a corner writes its zero-area
# slivers out anyway, and they read as winding failures
verts, faces = _gmt.drop_degenerate_faces(verts, faces)
surface = mesh3d(verts, faces,
_gmt.vertex_normals(verts, faces))
return surface
def _contour_fields(self, values, fields, value, label, supersample=0):
"""Several bodies from one cut and one paint: where each of the
fields `fields` reads off the columns of `values` exceeds `value`,
closed against the box.
What a mesh set contours a categorical realization through, each
category's field its draw against the best of its rivals'. Read off
the draws' corner means, two categories meeting along a contact
take their surfaces through the same points, since there the one's
field is the other's negated. Each category contoured on a field
taken block by block, the best rival picked inside every block
before the corners average, drew each contact twice, and a
realization's rock bodies overlapped and left gaps.
`fields` takes an array of one column per column of `values` and
returns one per field, row by row, and must never read below the
lowest it reads off a block, as a mean of the columns against the
most of the means cannot. Returns one entry per field, as
`_contour_values(..., close="above", fallback=False)` returns it:
None where the field misses the level. A field whose surface
Manifold will not take is contoured on its own instead, from the
field `fields` reads off each block.
"""
value = float(value)
values = _np.asarray(values, dtype=float)
blockwise = _np.asarray(fields(values))
origin, size, step, cell_values = self._cut_to_contour(
values, value, supersample=supersample, fields=fields)
if origin is None:
origin, size = self._origin, self._size
step, cell_values = self.base_step, values
origin, size, cell_values, _ = self._with_ghosts(
origin, size, step, cell_values)
painted = _painted_contour(origin, size, step, cell_values,
self.box_corner, value, label,
_closing_value("above", blockwise),
margin=1, fields=fields)
found = []
for m, (verts, faces) in enumerate(painted):
if len(faces):
clipped = _clip_to_box(verts, faces, *self._lattice_box())
if clipped is None:
found.append(self._contour_values(
blockwise[:, m], value, label,
supersample=supersample, close=True,
fallback=False))
continue
verts, faces = clipped
found.append(self._surface_of(verts, faces) if len(faces)
else None)
return found
def _with_ghosts(self, origin, size, step, cell_values):
"""The blocks and the shell of ghosts mirroring them, each ghost
holding a copy of its partner's value, and where the ghosts start."""
# the box measured in whichever lattice the mesh is drawn on
shape = _np.asarray(self.lattice_shape) * _np.round(
_np.asarray(self.base_step) / _np.asarray(step)).astype(int)
ghost_origin, ghost_size, ghost_parent = _ghost_shell(
origin, size, shape)
cell_values = _np.asarray(cell_values)
return (_np.concatenate([_np.asarray(origin),
_np.asarray(ghost_origin)]),
_np.concatenate([_np.asarray(size), _np.asarray(ghost_size)]),
_np.concatenate([cell_values, cell_values[ghost_parent]]),
len(origin))
def _lattice_box(self):
"""The model's box in the lattice frame, lowest corner and highest."""
low = _np.asarray(self.box_corner, dtype=float)
return low, low + _np.asarray(self.lattice_shape) * _np.asarray(
self.base_step, dtype=float)
def _surface_of(self, verts, faces):
"""A painted contour as the mesh its geometry calls for."""
# the lattice frame is where the contour is exact; a rotated
# set turns the finished vertices, as its `_hex_mesh` does
verts = self._to_world(verts)
verts, faces = _gmt.drop_degenerate_faces(verts, faces)
surface = mesh3d(verts, faces,
_gmt.vertex_normals(verts, faces))
if type(surface) is Mesh3D and surface.n_data > 0:
# Where the kept ground thins to a layer against the box, the
# body can touch itself along an edge -- the layer's underside
# meeting a cap painted onto the faces, or two pieces of the
# cut meeting on a face -- which comes out as one edge four
# triangles share, and no winding repair settles that. Split,
# the two touch rather than share it -- the copies moved apart
# as the booleans move theirs.
split, parted = _gmt.split_touching_edges(verts, faces)
if split is not verts:
split = _separated(split, parted)
touching = mesh3d(split, parted,
_gmt.vertex_normals(split, parted))
if type(touching) is not Mesh3D:
surface = touching
return surface
[docs]
class RotatedBlockSet3D(BlockSet3D):
"""
A variable-size block model rotated about its starting block.
The lattice is `BlockSet3D`'s, untouched: splitting, grouping, the
refinement criteria and the integer arithmetic all happen in the
unrotated frame, which is what keeps them exact. The rotation is applied
where coordinates *leave* -- the block centres, the sub-block fan-out a
prediction reads, the exported hexahedra -- and removed where coordinates
*come in* (`index_data`, and so `aggregate`). Every mesh test and
assignment reads sub-block positions through `get_batched_coordinates`,
so geometry against surfaces and solids works in world coordinates with
nothing overridden.
"""
def __init__(self, start, n, step, azimuth=0.0, dip=0.0, rake=0.0,
discretization=(2, 2, 2), max_levels=3,
labels=("X", "Y", "Z")):
self.azimuth = float(azimuth)
self.dip = float(dip)
self.rake = float(rake)
# the pivot: the first coarse block's centre, as `RotatedGrid3D`
# turns about its own origin
self._pivot = _np.asarray(start, dtype=float)
super().__init__(start, n, step, discretization=discretization,
max_levels=max_levels, labels=labels)
[docs]
def rotation_matrix(self):
return _gmt.rotation_matrix(self.azimuth, self.dip, self.rake)
def _to_world(self, coordinates):
return _np.matmul(coordinates - self._pivot,
self.rotation_matrix()) + self._pivot
def _centres(self):
return self._to_world(super()._centres())
def _lattice_bounding_box(self):
corner = self.box_corner
far = corner + self.lattice_shape * self.base_step
corners = _np.array(list(_iter.product(
*zip(corner, far))), dtype=float)
return BoundingBox.from_array(self._to_world(corners))
[docs]
def index_data(self, data):
# into the lattice frame: the inverse of the map the centres left by
return super().index_data(
rotate(data, self._pivot, self.azimuth, self.dip, self.rake,
reverse=True))
[docs]
def get_batched_coordinates(self, index=None):
if index is None:
index = _np.arange(self._n_data)
# fan out in the lattice frame, where a sub-block is an axis-aligned
# offset, then rotate every row: the offsets turn with the blocks
centres = BlockSet3D._centres(self)[index]
size = self.block_size[index]
coords = centres[:, None, :] + self._sub_grid[None, :, :] \
* size[:, None, :]
coords = self._to_world(coords.reshape(-1, 3))
splits = None if self.rows_per_location == 1 else len(centres)
return coords, splits
def _hex_mesh(self, origin, size, step):
# built on the lattice, turned as one piece: the welding is by shared
# corner indices, which a rotation cannot tear
mesh = super()._hex_mesh(origin, size, step)
mesh.points = self._to_world(_np.asarray(mesh.points, dtype=float))
return mesh
def _coarsest(self, n, step):
# the same angles about the same point, so every block that was never
# split lands where it was
return RotatedBlocks3D(start=self._pivot, n=n, step=step,
azimuth=self.azimuth, dip=self.dip,
rake=self.rake,
# added to the grid's arguments by `_blockdata`
discretization=list( # type: ignore[call-arg]
self.discretization),
labels=list(self.coordinate_labels))
[docs]
@classmethod
def from_data(cls, data, step, margin=0.1, decimals=0,
discretization=(2, 2, 2), max_levels=3):
"""
A rotated block model fitted to another object's spread.
As `RotatedGrid3D.from_data` -- the rotation fitted to the points and
rounded to `decimals` (degrees) before the box is measured -- counting
blocks rather than nodes.
"""
points, azimuth, dip, rake = _fitted_rotation(data, decimals)
mat = _gmt.rotation_matrix(azimuth, dip, rake)
centre = _np.mean(points, axis=0, keepdims=True)
unrotated = _np.matmul(points - centre, mat.T)
margin = _np.asarray(margin, dtype=float)
if margin.shape == ():
margin = _np.full([2, 3], float(margin))
elif margin.shape == (2,):
margin = _np.stack([margin] * 3, axis=1)
low = unrotated.min(axis=0)
high = unrotated.max(axis=0)
extent = high - low
low = low - extent * margin[0]
high = high + extent * margin[1]
step = _np.broadcast_to(
_np.asarray(step, dtype=float).ravel(), (3,))
n = _np.maximum(
_np.ceil((high - low) / step - _TOL_COVER).astype(int), 1)
# `start` is the first block's centre, half a step past the corner
start = _np.squeeze(
_np.matmul((low + step / 2)[None, :], mat) + centre)
start = _np.round(start, decimals)
labels = getattr(data, "coordinate_labels", None)
return cls(start=start, n=n, step=_np.array(step), azimuth=azimuth,
dip=dip, rake=rake, discretization=discretization,
max_levels=max_levels,
labels=labels if labels else ("X", "Y", "Z"))