Geometry and Coordinate Transformations

Geometry and Coordinate Transformations#

Overview#

parq-blockmodel uses a precise metadata-driven geometry model to define block model layout. All spatial coordinates are derived from:

  1. Geometry metadata (stored in Parquet key-value metadata)

  2. Logical indices (i, j, k) using canonical C-order

  3. Axis orientation (rotation transformation to world space)

The diagram below illustrates the complete flow from geometry definition to world coordinates:

digraph G { rankdir=TD; newrank=true; splines=ortho; node [shape=box, style=rounded, fontname=Arial]; edge [fontname=Arial, fontsize=9]; // Input layer subgraph cluster_input { label="Local Geometry Definition"; style=filled; color=lightgrey; corner [label="corner (i₀, j₀, k₀)", style="rounded,filled", fillcolor=lightblue]; blocksize [label="block_size (dx, dy, dz)", style="rounded,filled", fillcolor=lightblue]; shape [label="shape (nᵢ, nⱼ, nₖ)", style="rounded,filled", fillcolor=lightblue]; } // Logical indices layer subgraph cluster_logical { label="Logical Indices (Canonical C-order)"; style=filled; color=lightyellow; ijk [label="(i, j, k) 0 ≤ i < nᵢ 0 ≤ j < nⱼ 0 ≤ k < nₖ", style="rounded,filled", fillcolor=lightyellow]; rowindex [label="Flat Row Index (block_id) r = i + nᵢ(j + nⱼ·k)", style="rounded,filled", fillcolor=lightyellow]; } // Local coordinate layer subgraph cluster_local { label="Local Coordinates"; style=filled; color=lightgreen; local [label="Local Grid Position u = u₀ + (i + 0.5)·dx v = v₀ + (j + 0.5)·dy w = w₀ + (k + 0.5)·dz", style="rounded,filled", fillcolor=lightgreen]; } // World frame definition layer subgraph cluster_worldframe { label="World Frame Definition"; style=filled; color=lightcyan; origin [label="origin (x₀, y₀, z₀)", style="rounded,filled", fillcolor=azure]; axes [label="Axis Vectors axis_u, axis_v, axis_w", style="rounded,filled", fillcolor=lightcyan]; } // Transformation layer subgraph cluster_transform { label="Axis Transformation"; style=filled; color=lightpink; rotation [label="Rotation Matrix R R = [axis_u | axis_v | axis_w]ᵀ R⁻¹ = Rᵀ (orthonormal)", style="rounded,filled", fillcolor=lightsalmon]; } // World coordinates layer subgraph cluster_world { label="World Coordinates"; style=filled; color=lavender; world [label="(x, y, z)\ncentroid positions\nin world space\n\nworld = origin + R @ local", style="rounded,filled", fillcolor=lavender]; world_id [label="world_id = encode_coordinates(x, y, z)\nstable 64-bit positional identifier", style="rounded,dashed,filled", fillcolor=aliceblue]; } // Edges - Input to Logical corner -> ijk [label="define grid"]; blocksize -> ijk [label="cell size"]; shape -> ijk [label="bounds"]; ijk -> rowindex [label="flatten C-order"]; // Edges - Logical to Local ijk -> local [label="apply offsets"]; blocksize -> local [label="scale"]; corner -> local [label="translate"]; // Edges - Local to World local -> rotation [label="rotate local frame"]; axes -> rotation [label="define orientation", constraint=false]; origin -> world [label="translate", constraint=false]; rotation -> world [label="R @ local"]; world -> world_id [style=dashed, label="optional", constraint=false]; // Inverse transformations (dashed) world -> rotation [style=dashed, label="R⁻¹ (inverse)"]; rotation -> local [style=dashed]; local -> ijk [style=dashed]; // Layout helpers – keep world frame left of transform {rank=same; corner; blocksize; shape} {rank=same; ijk; rowindex} {rank=same; local} {rank=same; origin; axes} {rank=same; world; world_id} }

Geometry metadata → Logical indices (i,j,k) → Local coordinates → Rotation → World coordinates (x,y,z)#

Geometry Definition#

The parq_blockmodel.geometry.RegularGeometry class encapsulates all geometry metadata by combining two internal components:

Together, they define all block model geometry. The public API exposes only RegularGeometry; the internal split is transparent for most use cases.

Key metadata fields:

corner (u₀, v₀, w₀)

The local-space position of the (i=0, j=0, k=0) block corner (managed by LocalGeometry).

origin (x₀, y₀, z₀)

The world-space position of local coordinate (0, 0, 0) (managed by WorldFrame).

block_size (dx, dy, dz)

The size of each block along the logical/local axes. The public tuple is written as (dx, dy, dz), but for rotated geometries these are spacings along local i/j/k (equivalently u/v/w), not fixed world X/Y/Z dimensions (managed by LocalGeometry).

shape (nᵢ, nⱼ, nₖ)

The number of blocks along each logical axis (managed by LocalGeometry).

axis_u, axis_v, axis_w

Orthonormal basis vectors that define the orientation of the logical axes in world space (managed by WorldFrame).

Default (unrotated):

  • axis_u = (1, 0, 0) → u-axis aligned with world X

  • axis_v = (0, 1, 0) → v-axis aligned with world Y

  • axis_w = (0, 0, 1) → w-axis aligned with world Z

srs (optional)

Spatial Reference System (CRS) identifier for geographic context (managed by WorldFrame).

Logical Indices (i, j, k)#

Block indices are stored and processed using C-order (row-major) conventions:

  • i varies fastest (column index)

  • j varies next (row index)

  • k varies slowest (depth index)

The flat row index is computed as:

\[r = i + n_i \cdot (j + n_j \cdot k)\]

This C-order layout:

✓ Matches NumPy’s default behavior (np.ravel, np.unravel_index) ✓ Aligns with Parquet row-oriented storage ✓ Enables efficient dense array operations

Note

Do not confuse C-order with lexicographic MultiIndex sorting (x, y, z). Lexicographic sorting is not C-order and should not define canonical storage.

Local Coordinates#

Given logical indices (i, j, k), the local grid position is computed as:

\[\begin{split}u &= u_0 + (i + 0.5) \cdot dx \\ v &= v_0 + (j + 0.5) \cdot dy \\ w &= w_0 + (k + 0.5) \cdot dz\end{split}\]

The offset 0.5 positions each coordinate at the block centroid (not corner).

Note

Keeping the familiar names dx, dy, and dz is intentional: they represent physical spacing values, but they belong to the local lattice. Rotation via axis_u, axis_v, and axis_w determines how those local spacings appear in world coordinates.

Axis Orientation and Rotation#

The logical (u, v, w) coordinates are transformed to world-space (x, y, z) via a rotation matrix:

\[\begin{split}R = \begin{bmatrix} \text{axis\_u} \\ \text{axis\_v} \\ \text{axis\_w} \end{bmatrix}^T\end{split}\]

The world coordinates are then:

\[\begin{split}\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} o_x \\ o_y \\ o_z \end{bmatrix} + R \cdot \begin{bmatrix} u \\ v \\ w \end{bmatrix}\end{split}\]

Key properties:

  • Since axis vectors are orthonormal, the matrix is orthogonal: \(R^{-1} = R^T\)

  • The inverse transformation (world → local) is: \(\begin{bmatrix} u \\ v \\ w \end{bmatrix} = R^T \cdot (\begin{bmatrix} x \\ y \\ z \end{bmatrix} - \begin{bmatrix} o_x \\ o_y \\ o_z \end{bmatrix})\)

  • For unrotated geometries, R is the identity matrix (no rotation applied)

Once centroid coordinates are available in world space, they can also be converted into world_id using parq_blockmodel.utils.spatial_encoding.encode_world_coordinates (or the generic encode_frame_coordinates alias).

world_id is a stable 64-bit positional identifier derived from world-space XYZ coordinates within a given CRS. It uniquely identifies a block by location and is intended for joins, tracking, and external references. It is not a spatial index and should not be used for spatial range queries.

This world_id derivation is shown in the diagram as an optional downstream step from centroid positions.

Cross-model world_id use case#

world_id is useful when you work with multiple block models in the same SRS and want a stable, compact positional key for joins.

Within a shared encoding contract, world_id is unique across those models when:

  • all models use the same SRS,

  • all models use the same world_id_encoding metadata contract (scale, bit layout, frame, and axis order), and

  • the models do not contain overlapping block centroids.

Important

world_id uniqueness is not a pure CRS property by itself. It depends on both centroid positions and a compatible encoding definition (notably offsets and scale).

Decode world_id from a DataFrame (no ParquetBlockModel)#

If a user only has a pandas DataFrame with a world_id column and metadata in df.attrs['parq-blockmodel'], they can decode centroids directly:

import pandas as pd
from parq_blockmodel.utils import get_id_encoding_params, decode_frame_coordinates

# Example shape:
# df.columns includes: world_id, grade, density, ...
# df.attrs["parq-blockmodel"] includes: world_id_encoding metadata payload

meta = df.attrs["parq-blockmodel"]
encoding = meta["world_id_encoding"]

offset, scale = get_id_encoding_params(encoding)
x, y, z = decode_frame_coordinates(
    df["world_id"].to_numpy(dtype="int64"),
    offset=offset,
    scale=scale,
)

df_decoded = df.copy()
df_decoded["x"] = x
df_decoded["y"] = y
df_decoded["z"] = z

# Optional: set xyz as index
df_decoded = df_decoded.set_index(["x", "y", "z"])

Conversion APIs#

The parq_blockmodel.geometry.RegularGeometry class provides several conversion methods:

Forward (logical → world):

Inverse (world → logical):

Index conversion:

Export:

Metadata Storage#

All geometry metadata is embedded in Parquet files under a reserved key (default: "parq-blockmodel").

This approach:

✓ Eliminates the need for external sidecar files ✓ Keeps data and metadata synchronized ✓ Enables seamless data interchange ✓ Supports efficient metadata-driven processing

Metadata is serialized as JSON and stored in the Parquet file’s key-value metadata:

{
  "schema_version": "1.0",
  "corner": [100.0, 200.0, 300.0],
  "origin": [0.0, 0.0, 0.0],
  "block_size": [10.0, 10.0, 5.0],
  "shape": [50, 50, 20],
  "axis_u": [1.0, 0.0, 0.0],
  "axis_v": [0.0, 1.0, 0.0],
  "axis_w": [0.0, 0.0, 1.0],
  "srs": null,
  "world_id_encoding": {
    "enabled": true,
    "column": "world_id",
    "frame": "world_xyz",
    "axis_order": ["x", "y", "z"],
    "quantization": {"scale": 10.0},
    "offset": {"x": 500000.0, "y": 7000000.0, "z": 0.0}
  }
}

Examples#

Create a regular (unrotated) geometry:

from parq_blockmodel.geometry import RegularGeometry, LocalGeometry, WorldFrame

geom = RegularGeometry(
    local=LocalGeometry(
        corner=(100.0, 200.0, 300.0),
        block_size=(10.0, 10.0, 5.0),
        shape=(50, 50, 20),
    ),
    world=WorldFrame(origin=(0.0, 0.0, 0.0))
)

# Or more simply, using keyword arguments (backward compatible):
geom = RegularGeometry(
    corner=(100.0, 200.0, 300.0),
    block_size=(10.0, 10.0, 5.0),
    shape=(50, 50, 20),
)

# Convert logical index to world coordinates
x, y, z = geom.xyz_from_ijk(i=10, j=5, k=2)

Create a rotated geometry:

import numpy as np
from parq_blockmodel.geometry import RegularGeometry

# Define 30-degree rotation around Z-axis
angle = np.radians(30)
cos_a, sin_a = np.cos(angle), np.sin(angle)

geom = RegularGeometry(
    corner=(0.0, 0.0, 0.0),
    block_size=(10.0, 10.0, 5.0),
    shape=(50, 50, 20),
    axis_u=(cos_a, sin_a, 0.0),    # Rotated u-axis
    axis_v=(-sin_a, cos_a, 0.0),   # Rotated v-axis
    axis_w=(0.0, 0.0, 1.0),        # Z-axis unchanged
)

See Also#