mortie.tools

Encoding, decoding, inspection, and buffering of packed morton words.

Functions for morton indexing.

`geo2mort(lats, lons, order=None, points=None)`

Compute morton indices from geographic coordinates.

The entire pipeline runs in Rust via the healpix crate — no Python HEALPix backend is needed.

lat/lon inputs are treated as points by default (indeterminate resolution, encoded at max precision), so a bare geo2mort(lats, lons) returns order-29 Kind::Point words. Passing an explicit order asks for an area cell at that resolution instead (points inferred False). The two flags resolve as:

order=None, points=None (bare call) -> order-29 point words;
an explicit order with points unset -> area cell at order;
points=True -> order-29 point words (order-29-only; an explicit order != 29 raises ValueError, matching :meth:MortonIndexArray.from_latlon);
points=False -> area cell at order (order=None -> 29).

Non-finite lat/lon encode to the reserved empty word 0 (base cell 0 is the null sentinel) on both the area and point routes.

Parameters:

Name	Type	Description	Default
`lats`	`array - like`	Latitude(s) in degrees.	required
`lons`	`array - like`	Longitude(s) in degrees.	required
`order`	`int`	HEALPix order (0-29). Defaults to 29. An explicit value implies an area cell unless `points=True` is also given.	`None`
`points`	`bool`	Encode `Kind::Point` (order-29) vs `Kind::Area` words. Defaults to `True` for a bare call and `False` when an `order` is given.	`None`

Returns:

Type	Description
`ndarray`	Packed `uint64` morton word(s), same shape family as the input (scalar in -> length-1 ndarray).

Raises:

Type	Description
`ValueError`	If `points=True` is combined with an explicit `order != 29`.

Source code in mortie/tools.py

def geo2mort(lats, lons, order=None, points=None):
    """Compute morton indices from geographic coordinates.

    The entire pipeline runs in Rust via the ``healpix`` crate — no
    Python HEALPix backend is needed.

    lat/lon inputs are treated as **points** by default (indeterminate
    resolution, encoded at max precision), so a bare ``geo2mort(lats, lons)``
    returns order-29 ``Kind::Point`` words. Passing an explicit ``order`` asks
    for an **area** cell at that resolution instead (``points`` inferred
    ``False``). The two flags resolve as:

    * ``order=None, points=None`` (bare call) -> order-29 **point** words;
    * an explicit ``order`` with ``points`` unset -> **area** cell at ``order``;
    * ``points=True`` -> order-29 point words (order-29-only; an explicit
      ``order != 29`` raises ``ValueError``, matching
      :meth:`MortonIndexArray.from_latlon`);
    * ``points=False`` -> area cell at ``order`` (``order=None`` -> 29).

    Non-finite ``lat``/``lon`` encode to the reserved empty word ``0`` (base
    cell 0 is the null sentinel) on both the area and point routes.

    Parameters
    ----------
    lats : array-like
        Latitude(s) in degrees.
    lons : array-like
        Longitude(s) in degrees.
    order : int, optional
        HEALPix order (0-29). Defaults to 29. An explicit value implies an area
        cell unless ``points=True`` is also given.
    points : bool, optional
        Encode ``Kind::Point`` (order-29) vs ``Kind::Area`` words. Defaults to
        ``True`` for a bare call and ``False`` when an ``order`` is given.

    Returns
    -------
    ndarray
        Packed ``uint64`` morton word(s), same shape family as the input
        (scalar in -> length-1 ndarray).

    Raises
    ------
    ValueError
        If ``points=True`` is combined with an explicit ``order != 29``.
    """
    # Resolve the point/area mode: a bare call encodes points; an explicit order
    # implies an area cell at that resolution unless the caller forces points.
    if points is None:
        points = order is None
    if order is None:
        order = MAX_ORDER
    if points and int(order) != MAX_ORDER:
        raise ValueError(
            "points=True encodes an order-29 point; pass order=29 "
            "(the default) or omit it"
        )
    # Ensure contiguous arrays for Rust FFI
    if not np.isscalar(lats):
        lats = np.ascontiguousarray(lats, dtype=np.float64)
        lons = np.ascontiguousarray(lons, dtype=np.float64)
    result = _rust_geo2mort(lats, lons, int(order), points)
    # Always return a contiguous uint64 ndarray. The scalar Rust path hands back
    # a Python int (which np would otherwise infer as int64), so coerce to keep
    # the dtype uint64 regardless of scalar-vs-array input or hemisphere.
    return np.ascontiguousarray(np.atleast_1d(result), dtype=np.uint64)

`mort2geo(morton)`

Convert morton index to lat/lon of pixel center.

This is the inverse of geo2mort, returning the center coordinates of the HEALPix cell identified by the morton index.

Mixed-order arrays are supported (issue #116): elements are grouped by order (:func:orders_of), each group runs the uniform kernel, and the results scatter back to input positions. Point words (spec §4) are order 29 by definition and group with order 29 — a point's location is exactly what mort2geo returns.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Morton index (mixed orders allowed).	required

Returns:

Name	Type	Description
`lat`	`float or array`	Latitude in degrees
`lon`	`float or array`	Longitude in degrees

Source code in mortie/tools.py

def mort2geo(morton):
    """Convert morton index to lat/lon of pixel center.

    This is the inverse of geo2mort, returning the center coordinates
    of the HEALPix cell identified by the morton index.

    Mixed-order arrays are supported (issue #116): elements are grouped by
    order (:func:`orders_of`), each group runs the uniform kernel, and the
    results scatter back to input positions. Point words (spec §4) are order
    29 by definition and group with order 29 — a point's location is exactly
    what mort2geo returns.

    Parameters
    ----------
    morton : int or array-like
        Morton index (mixed orders allowed).

    Returns
    -------
    lat : float or array
        Latitude in degrees
    lon : float or array
        Longitude in degrees
    """
    # Handle scalar vs array input to match geo2mort behavior
    input_is_scalar = np.isscalar(morton)

    # Group-by-order dispatch for mixed-order input (issue #116).
    if not input_is_scalar:
        words = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
        orders = orders_of(words)
        unique_orders = np.unique(orders)
        if unique_orders.size > 1:
            lat = np.empty(words.size, dtype=np.float64)
            lon = np.empty(words.size, dtype=np.float64)
            for order in unique_orders:
                mask = orders == order
                lat[mask], lon[mask] = mort2geo(words[mask])
            return lat, lon

    # Decode morton to normalized address and parent
    normed, parent, order = mort2norm(morton)

    # Convert to UNIQ
    uniq = norm2uniq(normed, parent, order)

    # Convert to lat/lon (uniq2geo decodes the order from the UNIQ value)
    lat, lon = uniq2geo(uniq)

    # Return array to match geo2mort behavior
    if input_is_scalar:
        return np.array([lat]), np.array([lon])
    return lat, lon

`mort2bbox(morton)`

Convert morton index to bounding box of the pixel.

For pixels touching the antimeridian, vertex longitudes at ±180° are normalized to use consistent representation based on hemisphere voting, preventing bbox misinterpretation as spanning the entire globe.

Mixed-order arrays are supported (issue #116): elements are grouped by order (:func:orders_of), each group runs the uniform kernel, and the results scatter back to input positions. Point words (spec §4) are order 29 by definition and group with order 29 — a point yields the bounding box of its containing order-29 cell (the cell that contains the point), which is exactly the bbox of the order-29 area word at the same location. A group of points therefore covers a well-defined area, element by element.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Morton index (mixed orders allowed).	required

Returns:

Name	Type	Description
`bbox`	`dict or list of dicts`	Bounding box in format suitable for STAC/CMR: {"west": min_lon, "south": min_lat, "east": max_lon, "north": max_lat}

Source code in mortie/tools.py

def mort2bbox(morton):
    """Convert morton index to bounding box of the pixel.

    For pixels touching the antimeridian, vertex longitudes at ±180° are
    normalized to use consistent representation based on hemisphere voting,
    preventing bbox misinterpretation as spanning the entire globe.

    Mixed-order arrays are supported (issue #116): elements are grouped by
    order (:func:`orders_of`), each group runs the uniform kernel, and the
    results scatter back to input positions. Point words (spec §4) are order
    29 by definition and group with order 29 — a point yields the bounding box
    of its containing order-29 cell (the cell that contains the point), which
    is exactly the bbox of the order-29 **area** word at the same location. A
    group of points therefore covers a well-defined area, element by element.

    Parameters
    ----------
    morton : int or array-like
        Morton index (mixed orders allowed).

    Returns
    -------
    bbox : dict or list of dicts
        Bounding box in format suitable for STAC/CMR:
        {"west": min_lon, "south": min_lat, "east": max_lon, "north": max_lat}
    """
    morton = np.atleast_1d(morton)
    is_scalar = len(morton) == 1

    words = np.asarray(morton, dtype=np.uint64)
    # Group-by-order dispatch for mixed-order input (issue #116).
    orders = orders_of(words)
    unique_orders = np.unique(orders)
    if unique_orders.size > 1:
        bboxes = [None] * words.size
        for order in unique_orders:
            (idx,) = np.nonzero(orders == order)
            group = mort2bbox(words[idx])
            if idx.size == 1:
                bboxes[idx[0]] = group  # length-1 call returns the bare dict
            else:
                for i, bbox in zip(idx, group):
                    bboxes[i] = bbox
        return bboxes

    # First get the pixel center
    normed, parent, order = mort2norm(morton)
    uniq = norm2uniq(normed, parent, order)

    nside = 2**order
    nest = uniq - 4 * (nside**2)

    # Get pixel boundaries: (N, 3, 4) — cell in axis 0, xyz in axis 1, the 4
    # corners in axis 2.  A single cell comes back 2-D (3, 4); promote it.
    boundaries = hp.boundaries(order, nest)
    if boundaries.ndim == 2:
        boundaries = boundaries[np.newaxis, ...]
    n = len(morton)

    # One batched vec2ang over every cell's corners (one Rust round-trip instead
    # of one per cell), then reshape to (N, 4).
    verts = np.transpose(boundaries, (0, 2, 1)).reshape(-1, 3)
    theta, phi = hp.vec2ang(verts)
    lats_all = (90 - np.degrees(theta)).reshape(n, 4)
    lons_all = np.degrees(phi)
    lons_all = np.where(lons_all > 180, lons_all - 360, lons_all).reshape(n, 4)

    bboxes = []
    for i in range(n):
        lats = lats_all[i]
        lons = lons_all[i]

        # Normalize antimeridian representation
        # Check if bbox touches antimeridian with mixed ±180°
        ANTIMERIDIAN_TOLERANCE = 1e-6
        on_antimeridian = np.abs(np.abs(lons) - 180.0) < ANTIMERIDIAN_TOLERANCE

        if np.any(on_antimeridian) and (np.max(lons) - np.min(lons)) > 180:
            # Count vertices in each hemisphere (excluding those on antimeridian)
            non_antimeridian = ~on_antimeridian
            if np.any(non_antimeridian):
                western_count = np.sum(lons[non_antimeridian] < -0.1)
                eastern_count = np.sum(lons[non_antimeridian] > 0.1)

                # Determine target longitude for antimeridian vertices
                if western_count > eastern_count:
                    target_lon = -180.0
                elif eastern_count > western_count:
                    target_lon = 180.0
                else:
                    # Use median of non-antimeridian lons
                    median_lon = np.median(lons[non_antimeridian])
                    target_lon = -180.0 if median_lon < 0 else 180.0

                # Normalize antimeridian vertices
                lons = lons.copy()
                lons[on_antimeridian] = target_lon

        # Create bounding box
        bbox = {
            "west": float(np.min(lons)),
            "south": float(np.min(lats)),
            "east": float(np.max(lons)),
            "north": float(np.max(lats))
        }
        bboxes.append(bbox)

    if is_scalar:
        return bboxes[0]
    return bboxes

`mort2polygon(morton, step=1)`

Convert morton index to polygon representation.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Morton index.	required
`step`	`int`	Points per side for the cell boundary (default 1 = 4 corners). Use step=32 for 128 boundary points that accurately trace curved cell edges, important for polar cells where 4-corner polygons poorly approximate the true HEALPix boundary.	`1`

Returns:

Name	Type	Description
`polygon`	`list or list of lists`	Polygon coordinates as [[lat, lon], ...] in standard geographic order. The polygon is closed (first point repeated at end). Note: Returns [lat, lon] pairs, NOT [lon, lat]. This is the standard geographic coordinate order used by most spatial analysis libraries.

Notes

Polygons that touch the antimeridian (±180° longitude) are automatically normalized to use consistent longitude representation (-180 or +180) based on which hemisphere contains the majority of vertices. This prevents spatial libraries from misinterpreting touching polygons as crossing polygons.

Mixed-order arrays are supported (issue #116): elements are grouped by order (:func:orders_of), each group runs the uniform kernel, and the results scatter back to input positions (rings are 4step+1 vertices at every order, so mixed orders do not change the output shape). Point words (spec §4) are order 29 by definition and group with order 29 — a point yields the polygon ring of its containing order-29 cell, exactly the ring of the order-29 area* word at the same location.

Source code in mortie/tools.py

def mort2polygon(morton, step=1):
    """Convert morton index to polygon representation.

    Parameters
    ----------
    morton : int or array-like
        Morton index.
    step : int, optional
        Points per side for the cell boundary (default 1 = 4 corners).
        Use step=32 for 128 boundary points that accurately trace
        curved cell edges, important for polar cells where 4-corner
        polygons poorly approximate the true HEALPix boundary.

    Returns
    -------
    polygon : list or list of lists
        Polygon coordinates as [[lat, lon], ...] in standard geographic order.
        The polygon is closed (first point repeated at end).

        **Note**: Returns [lat, lon] pairs, NOT [lon, lat]. This is the standard
        geographic coordinate order used by most spatial analysis libraries.

    Notes
    -----
    Polygons that touch the antimeridian (±180° longitude) are automatically
    normalized to use consistent longitude representation (-180 or +180) based
    on which hemisphere contains the majority of vertices. This prevents spatial
    libraries from misinterpreting touching polygons as crossing polygons.

    Mixed-order arrays are supported (issue #116): elements are grouped by
    order (:func:`orders_of`), each group runs the uniform kernel, and the
    results scatter back to input positions (rings are 4*step+1 vertices at
    every order, so mixed orders do not change the output shape). Point words
    (spec §4) are order 29 by definition and group with order 29 — a point
    yields the polygon ring of its containing order-29 cell, exactly the ring
    of the order-29 **area** word at the same location.
    """
    morton = np.atleast_1d(morton)
    is_scalar = len(morton) == 1

    words = np.asarray(morton, dtype=np.uint64)
    # Group-by-order dispatch for mixed-order input (issue #116).
    orders = orders_of(words)
    unique_orders = np.unique(orders)
    if unique_orders.size > 1:
        polygons = [None] * words.size
        for order in unique_orders:
            (idx,) = np.nonzero(orders == order)
            group = mort2polygon(words[idx], step=step)
            if idx.size == 1:
                polygons[idx[0]] = group  # length-1 call returns the bare ring
            else:
                for i, polygon in zip(idx, group):
                    polygons[i] = polygon
        return polygons

    # Get pixel information
    normed, parent, order = mort2norm(morton)
    uniq = norm2uniq(normed, parent, order)

    nside = 2**order
    nest = uniq - 4 * (nside**2)

    # Get pixel boundaries: (N, 3, 4*step) — cell in axis 0, xyz in axis 1, the
    # boundary points in axis 2.  A single cell comes back 2-D (3, ncols);
    # promote it.
    boundaries = hp.boundaries(order, nest, step=step)
    if boundaries.ndim == 2:
        boundaries = boundaries[np.newaxis, ...]
    n = len(morton)
    ncols = 4 * step

    # One batched vec2ang over every cell's boundary points (one Rust round-trip
    # instead of one per cell), then reshape to (N, ncols).
    verts = np.transpose(boundaries, (0, 2, 1)).reshape(-1, 3)
    theta, phi = hp.vec2ang(verts)
    lats_all = (90 - np.degrees(theta)).reshape(n, ncols)
    lons_all = np.degrees(phi)
    lons_all = np.where(lons_all > 180, lons_all - 360, lons_all).reshape(n, ncols)

    polygons = []
    for i in range(n):
        lats = lats_all[i]
        lons = lons_all[i]

        # Create polygon as list of [lat, lon] pairs (standard geographic order)
        # Close the polygon by repeating first point
        polygon = [[float(lats[j]), float(lons[j])] for j in range(len(lons))]
        polygon.append(polygon[0])  # Close the polygon

        # Normalize antimeridian representation to prevent misinterpretation
        polygon = _normalize_antimeridian_polygon(polygon)

        polygons.append(polygon)

    if is_scalar:
        return polygons[0]
    return polygons

`mort2healpix(morton)`

Convert morton index to HEALPix cell ID and order.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Morton index.	required

Returns:

Name	Type	Description
`cell_ids`	`int or ndarray`	HEALPix cell ID(s) in NESTED scheme
`order`	`int`	HEALPix order (resolution level)

Raises:

Type	Description
`ValueError`	If the words are at mixed orders (propagated from :func:`mort2norm`, which enforces the same-order precondition below).

Notes

The function converts morton indices to HEALPix NESTED scheme cell IDs. All input morton indices must be at the same order.

Examples:

>>> import mortie
>>> m = mortie.geo2mort(-80.0, 120.0, order=6)[0]
>>> cell_id, order = mortie.mort2healpix(m)
>>> print(f"HEALPix cell {cell_id} at order {order}")
HEALPix cell 37010 at order 6

Source code in mortie/tools.py

def mort2healpix(morton):
    """Convert morton index to HEALPix cell ID and order.

    Parameters
    ----------
    morton : int or array-like
        Morton index.

    Returns
    -------
    cell_ids : int or ndarray
        HEALPix cell ID(s) in NESTED scheme
    order : int
        HEALPix order (resolution level)

    Raises
    ------
    ValueError
        If the words are at mixed orders (propagated from :func:`mort2norm`,
        which enforces the same-order precondition below).

    Notes
    -----
    The function converts morton indices to HEALPix NESTED scheme cell IDs.
    All input morton indices must be at the same order.

    Examples
    --------
    >>> import mortie
    >>> m = mortie.geo2mort(-80.0, 120.0, order=6)[0]
    >>> cell_id, order = mortie.mort2healpix(m)
    >>> print(f"HEALPix cell {cell_id} at order {order}")
    HEALPix cell 37010 at order 6
    """
    # Check if input is scalar before converting to array
    is_scalar = np.isscalar(morton)
    morton = np.atleast_1d(morton)

    # Get normalized morton and order
    normed, parent, order = mort2norm(morton)

    # Convert to UNIQ indexing
    uniq = norm2uniq(normed, parent, order)

    # Convert UNIQ to HEALPix NESTED cell ID
    # UNIQ = 4 * nside^2 + nest_index
    nside = 2**order
    cell_ids = uniq - 4 * (nside**2)

    # Ensure arrays for consistent handling
    cell_ids = np.atleast_1d(cell_ids).astype(np.int64)
    order = np.atleast_1d(order)

    if is_scalar:
        return int(cell_ids[0]), int(order[0])

    # For array input, return single order if all are the same
    order_val = int(order[0]) if len(np.unique(order)) == 1 else order
    return cell_ids, order_val

`mort2norm(morton)`

Convert morton index back to normalized address and parent cell.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Packed morton word(s) (`uint64`; base cells 7-11 set bit 63).	required

Returns:

Name	Type	Description
`normed`	`int or array`	Normalized HEALPix address
`parent`	`int or array`	Parent base cell (0-11)
`order`	`int or array`	HEALPix order inferred from morton index

Raises:

Type	Description
`ValueError`	If the words are at mixed orders — the return contract carries a single scalar order, so use :func:`orders_of` for per-element orders.

Notes

Empty input returns two empty int64 arrays and order == 0.

Source code in mortie/tools.py

def mort2norm(morton):
    """Convert morton index back to normalized address and parent cell.

    Parameters
    ----------
    morton : int or array-like
        Packed morton word(s) (``uint64``; base cells 7-11 set bit 63).

    Returns
    -------
    normed : int or array
        Normalized HEALPix address
    parent : int or array
        Parent base cell (0-11)
    order : int or array
        HEALPix order inferred from morton index

    Raises
    ------
    ValueError
        If the words are at mixed orders — the return contract carries a single
        scalar order, so use :func:`orders_of` for per-element orders.

    Notes
    -----
    Empty input returns two empty ``int64`` arrays and ``order == 0``.
    """
    morton = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
    is_scalar = len(morton) == 1

    # Empty input: nothing to decode. Return empty int64 arrays (matching the
    # array-path dtype) and order 0.
    if morton.size == 0:
        empty = np.empty(0, dtype=np.int64)
        return empty, empty.copy(), 0

    # The packed-u64 kernel decodes each word to (nested, depth); the depth is
    # the HEALPix order (no decimal-digit scan). Reject mixed orders: the
    # return contract is a single scalar order (the geo kernels above this
    # dispatch group-by-order and never hit this — issue #116).
    nested, depths = _rust_mort2nested(np.ascontiguousarray(morton))
    if np.any(depths != depths[0]):
        raise ValueError(
            f"Mixed orders in morton array: {sorted(set(int(d) for d in depths))}; "
            "use orders_of for per-element orders"
        )

    order = int(depths[0])
    # nested ids are HEALPix cell ids (<< 2^58 for order <= 29), so int64 is safe
    # arithmetic here and keeps normed/parent signed for downstream callers.
    nested = nested.astype(np.int64)
    nside_sq = np.int64(1) << np.int64(2 * order)
    parent = nested // nside_sq
    normed = nested % nside_sq

    if is_scalar:
        return normed[0], parent[0], order
    return normed, parent, order

`norm2mort(normed, parent, order)`

Convert a normalized HEALPix address + base cell to a packed morton word.

The exact inverse of :func:mort2norm: mort2norm(norm2mort(n, p, o)) returns (n, p, o). Born order-29-native (issue #48) — there is no order cap beyond the kernel's MAX_ORDER of 29. The returned uint64 is the packed decimal_morton word (issue #58; the prefix is base+1, so bit 63 is set — a large unsigned value — for base cells 7-11), not the retired decimal encoding.

Parameters:

Name	Type	Description	Default
`normed`	`int or array`	Normalized HEALPix address (the in-base z-order, `0 <= normed < 4**order`).	required
`parent`	`int or array`	Parent base cell (0-11).	required
`order`	`int`	HEALPix order (0-29).	required

Returns:

Name	Type	Description
`morton`	`uint64 or ndarray`	Packed morton word(s).

Source code in mortie/tools.py

def norm2mort(normed, parent, order):
    """Convert a normalized HEALPix address + base cell to a packed morton word.

    The exact inverse of :func:`mort2norm`: ``mort2norm(norm2mort(n, p, o))``
    returns ``(n, p, o)``. Born order-29-native (issue #48) — there is no order
    cap beyond the kernel's ``MAX_ORDER`` of 29. The returned ``uint64`` is the
    packed ``decimal_morton`` word (issue #58; the prefix is ``base+1``, so bit 63
    is set — a large unsigned value — for base cells 7-11), not the retired
    decimal encoding.

    Parameters
    ----------
    normed : int or array
        Normalized HEALPix address (the in-base z-order, ``0 <= normed < 4**order``).
    parent : int or array
        Parent base cell (0-11).
    order : int
        HEALPix order (0-29).

    Returns
    -------
    morton : uint64 or ndarray
        Packed morton word(s).
    """
    normed = np.atleast_1d(np.asarray(normed, dtype=np.int64))
    parent = np.atleast_1d(np.asarray(parent, dtype=np.int64))
    is_scalar = normed.size == 1 and parent.size == 1
    # nested = parent * nside^2 + normed; pack via the kernel bridge.
    nested = (parent.astype(np.uint64) << np.uint64(2 * order)) | normed.astype(
        np.uint64
    )
    n = max(normed.size, parent.size)
    nested = np.ascontiguousarray(np.broadcast_to(nested, (n,)))
    depths = np.full(nested.size, order, dtype=np.uint8)
    morton = _rust_nested2mort(nested, depths)
    if is_scalar:
        return np.uint64(morton[0])
    return morton

`infer_order_from_morton(morton)`

Infer the single HEALPix order of packed morton word(s).

Decodes through the packed-u64 kernel (issue #48): the order is carried in the word's suffix, not in any decimal-digit count. The return is one scalar order, so array input must be uniform-order; mixed-order input raises, naming the distinct orders (issue #116 — previously the first element's order was returned silently). For per-element orders of a mixed array use :func:orders_of.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Packed morton word(s), all at one order.	required

Returns:

Type	Description
`int`	The HEALPix order.

Raises:

Type	Description
`ValueError`	If the words are at mixed orders.

Source code in mortie/tools.py

def infer_order_from_morton(morton):
    """Infer the single HEALPix order of packed morton word(s).

    Decodes through the packed-u64 kernel (issue #48): the order is carried in
    the word's suffix, not in any decimal-digit count. The return is one
    scalar order, so array input must be uniform-order; mixed-order input
    raises, naming the distinct orders (issue #116 — previously the first
    element's order was returned silently). For per-element orders of a mixed
    array use :func:`orders_of`.

    Parameters
    ----------
    morton : int or array-like
        Packed morton word(s), all at one order.

    Returns
    -------
    int
        The HEALPix order.

    Raises
    ------
    ValueError
        If the words are at mixed orders.
    """
    m = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
    _, depths = _rust_mort2nested(np.ascontiguousarray(m))
    distinct = np.unique(depths)
    if distinct.size > 1:
        raise ValueError(
            f"Mixed orders in morton array: {[int(d) for d in distinct]}; "
            "use orders_of for per-element orders"
        )
    return int(depths[0])

`orders_of(morton)`

Per-element HEALPix order of packed morton words.

Vectorized numpy decode of the 6-bit suffix (bits 5-0) per the spec page's suffix table (docs/specification.md §1):

suffix 0..=27 — variable-length area element; the order is the suffix value (0 = base-cell-only).
suffix 28..=47 — order-28/29 area cells in parent-first preorder suffix = 28 + t28*5 + (t29 present ? t29 + 1 : 0): each t28 owns a 5-block (the order-28 parent, then its four order-29 children), so (suffix - 28) % 5 == 0 is order 28 and everything else is order 29.
suffix 48..=63 — order-29 point (max-encoded, no area claim — spec §4); points are order 29 by definition.

Pure bit arithmetic — words are not validated (the empty sentinel 0 decodes as order 0; use :func:validate_morton to reject malformed words). This is the per-element, mixed-order-native counterpart of :func:infer_order_from_morton.

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Packed morton word(s) (`uint64`).	required

Returns:

Type	Description
`ndarray`	`uint8` order per element, 0-29 (scalar in -> length-1 ndarray, matching :func:`geo2mort`).

Source code in mortie/tools.py

def orders_of(morton):
    """Per-element HEALPix order of packed morton words.

    Vectorized numpy decode of the 6-bit suffix (bits 5-0) per the spec page's
    suffix table (``docs/specification.md`` §1):

    * suffix ``0..=27`` — variable-length area element; the order *is* the
      suffix value (``0`` = base-cell-only).
    * suffix ``28..=47`` — order-28/29 area cells in parent-first preorder
      ``suffix = 28 + t28*5 + (t29 present ? t29 + 1 : 0)``: each ``t28`` owns
      a 5-block (the order-28 parent, then its four order-29 children), so
      ``(suffix - 28) % 5 == 0`` is order 28 and everything else is order 29.
    * suffix ``48..=63`` — order-29 **point** (max-encoded, no area claim —
      spec §4); points are order 29 by definition.

    Pure bit arithmetic — words are not validated (the empty sentinel ``0``
    decodes as order 0; use :func:`validate_morton` to reject malformed
    words). This is the per-element, mixed-order-native counterpart of
    :func:`infer_order_from_morton`.

    Parameters
    ----------
    morton : int or array-like
        Packed morton word(s) (``uint64``).

    Returns
    -------
    ndarray
        ``uint8`` order per element, 0-29 (scalar in -> length-1 ndarray,
        matching :func:`geo2mort`).
    """
    m = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
    suffix = (m & np.uint64(0x3F)).astype(np.uint8)
    # 0..=27: order == suffix. 28..=47: order-28 on the 5-block parent slots,
    # order 29 otherwise. 48..=63: order-29 point.
    orders = suffix.copy()
    band = (suffix >= 28) & (suffix <= 47)
    orders[band] = np.where((suffix[band] - 28) % 5 == 0, 28, 29)
    orders[suffix >= 48] = 29
    return orders

`orders_of_uniq(uniq)`

Per-element HEALPix order decoded from UNIQ cell numbers.

UNIQ is self-describing: uniq = 4 * 4**order + nested with 0 <= nested < 12 * 4**order, so order-k values occupy exactly [4**(k+1), 4**(k+2)) and consecutive orders tile that line without gaps. The order is therefore a pure function of the value — no caller-supplied order is needed and mixed-resolution input decodes element by element (issue #136).

Implemented as an exact integer bucket search rather than the log2(uniq / 4) // 2 form this module used previously: the float64 round-trip is not exact above ~2**53, so e.g. 4**30 - 1 (the last order-28 value) rounds up to 4**30 and mis-decodes as order 29.

The UNIQ counterpart of :func:orders_of, and mirrors its contract: per-element, mixed-order-native, uint8 out, scalar in -> length-1 ndarray. One deliberate difference: :func:orders_of is pure bit arithmetic and never validates, because every 6-bit morton suffix decodes to some order. UNIQ has no such total decode -- a value outside [4, 4**31) names no cell at any order -- so this raises instead of inventing an answer.

Parameters:

Name	Type	Description	Default
`uniq`	`int or array - like`	UNIQ encoded cell number(s).	required

Returns:

Type	Description
`ndarray`	`uint8` order per element, 0-`MAX_ORDER` (scalar in -> length-1 ndarray, matching :func:`orders_of`).

Raises:

Type	Description
`ValueError`	If any value lies outside the UNIQ range for orders 0-`MAX_ORDER`.

Source code in mortie/tools.py

def orders_of_uniq(uniq):
    """Per-element HEALPix order decoded from UNIQ cell numbers.

    UNIQ is self-describing: ``uniq = 4 * 4**order + nested`` with
    ``0 <= nested < 12 * 4**order``, so order-``k`` values occupy exactly
    ``[4**(k+1), 4**(k+2))`` and consecutive orders tile that line without
    gaps. The order is therefore a pure function of the value — no
    caller-supplied order is needed and mixed-resolution input decodes element
    by element (issue #136).

    Implemented as an exact integer bucket search rather than the
    ``log2(uniq / 4) // 2`` form this module used previously: the float64
    round-trip is not exact above ~2**53, so e.g. ``4**30 - 1`` (the last
    order-28 value) rounds up to ``4**30`` and mis-decodes as order 29.

    The UNIQ counterpart of :func:`orders_of`, and mirrors its contract:
    per-element, mixed-order-native, ``uint8`` out, scalar in -> length-1
    ndarray. One deliberate difference: :func:`orders_of` is pure bit
    arithmetic and never validates, because every 6-bit morton suffix decodes
    to *some* order. UNIQ has no such total decode -- a value outside
    ``[4, 4**31)`` names no cell at any order -- so this raises instead of
    inventing an answer.

    Parameters
    ----------
    uniq : int or array-like
        UNIQ encoded cell number(s).

    Returns
    -------
    ndarray
        ``uint8`` order per element, 0-``MAX_ORDER`` (scalar in -> length-1
        ndarray, matching :func:`orders_of`).

    Raises
    ------
    ValueError
        If any value lies outside the UNIQ range for orders 0-``MAX_ORDER``.
    """
    # Cast defensively: `asarray(..., dtype=int64)` raises OverflowError for a
    # value above int64 and silently *truncates* a float, both of which would
    # bypass the ValueError this function documents. Normalize them here so the
    # contract holds for every input, not just int64-representable ones.
    arr = np.atleast_1d(np.asarray(uniq))
    if arr.dtype.kind == "f" and not np.all(np.equal(np.mod(arr, 1), 0)):
        raise ValueError(
            f"Not a valid UNIQ cell number for orders 0-{MAX_ORDER}: "
            f"{arr.ravel()[0]!r} is not an integer")
    if arr.dtype.kind == "u":
        # uint64 -> int64 *wraps* silently rather than raising, so an oversized
        # value would reach the range check as a meaningless negative and be
        # reported as such. Every wrap lands negative so nothing mis-decodes as
        # valid, but the message would name a number the caller never passed.
        over = arr > np.iinfo(np.int64).max
        if np.any(over):
            raise ValueError(
                f"Not a valid UNIQ cell number for orders 0-{MAX_ORDER}: "
                f"{int(arr[over].ravel()[0])} is out of the int64 range")
    try:
        u = np.atleast_1d(np.asarray(arr, dtype=np.int64))
    except (OverflowError, ValueError, TypeError) as exc:
        raise ValueError(
            f"Not a valid UNIQ cell number for orders 0-{MAX_ORDER}: "
            f"{uniq!r} is out of the int64 range") from exc
    # bounds[k] = 4**(k+1) is the first UNIQ value of order k; the trailing
    # entry closes order MAX_ORDER's range (4**31 still fits int64).
    bounds = np.int64(4) ** np.arange(1, MAX_ORDER + 3, dtype=np.int64)
    orders = (np.searchsorted(bounds, u, side='right') - 1).astype(np.int64)
    bad = (orders < 0) | (orders > MAX_ORDER)
    if np.any(bad):
        raise ValueError(
            f"Not a valid UNIQ cell number for orders 0-{MAX_ORDER}: "
            f"{int(u[bad][0])}")
    return orders.astype(np.uint8)

`is_point(morton)`

Per-element point-kind predicate for packed morton words.

Kind is carried by the encoding itself (spec §4): suffix 0..=47 decodes as an area word, suffix 48..=63 as an order-29 point (a location with no area claim — docs/specification.md §1 suffix table). Pure bit arithmetic; words are not validated (see :func:validate_morton).

Parameters:

Name	Type	Description	Default
`morton`	`int or array - like`	Packed morton word(s) (`uint64`).	required

Returns:

Type	Description
`ndarray`	`bool` per element, True for point words (scalar in -> length-1 ndarray, matching :func:`geo2mort`).

Source code in mortie/tools.py

def is_point(morton):
    """Per-element point-kind predicate for packed morton words.

    Kind is carried by the encoding itself (spec §4): suffix ``0..=47``
    decodes as an **area** word, suffix ``48..=63`` as an order-29 **point**
    (a location with no area claim — ``docs/specification.md`` §1 suffix
    table). Pure bit arithmetic; words are not validated (see
    :func:`validate_morton`).

    Parameters
    ----------
    morton : int or array-like
        Packed morton word(s) (``uint64``).

    Returns
    -------
    ndarray
        ``bool`` per element, True for point words (scalar in -> length-1
        ndarray, matching :func:`geo2mort`).
    """
    m = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
    return (m & np.uint64(0x3F)) >= np.uint64(48)

`validate_morton(morton, order=None)`

Validate that a packed morton word is well-formed.

The kernel decode rejects the empty sentinel (0) and any word with an invalid base-cell prefix; this also checks the decoded order matches order when one is supplied.

Parameters:

Name	Type	Description	Default
`morton`	`int`	Packed morton word to validate.	required
`order`	`int`	Expected HEALPix order. If None, no order check is made.	`None`

Returns:

Type	Description
`bool`	True if the word is a valid morton word.

Raises:

Type	Description
`ValueError`	If the word does not decode or its order disagrees with `order`.

Source code in mortie/tools.py

def validate_morton(morton, order=None):
    """Validate that a packed morton word is well-formed.

    The kernel decode rejects the empty sentinel (0) and any word with an
    invalid base-cell prefix; this also checks the decoded order matches
    ``order`` when one is supplied.

    Parameters
    ----------
    morton : int
        Packed morton word to validate.
    order : int, optional
        Expected HEALPix order. If None, no order check is made.

    Returns
    -------
    bool
        True if the word is a valid morton word.

    Raises
    ------
    ValueError
        If the word does not decode or its order disagrees with ``order``.
    """
    m = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
    # The kernel raises ValueError on the empty sentinel / an invalid prefix.
    _, depths = _rust_mort2nested(np.ascontiguousarray(m))
    decoded_order = int(depths[0])
    if order is not None and decoded_order != order:
        raise ValueError(
            f"Morton word decodes to order {decoded_order}, expected {order}"
        )
    return True

`clip2order(clip_order, midx)`

Coarsen packed morton words to a lower resolution.

Degrades each packed word to clip_order by coarsening it through the kernel (the inverse of refining): the base cell and the first clip_order tuples are kept, finer detail is dropped, and the suffix is rewritten. Words already at or below clip_order are returned unchanged.

The print_factor flag was removed for the 1.x freeze (issue #68). It returned 18 - clip_order, a level count anchored to the retired decimal encoding's order-18 ceiling, so it went negative for the order-19..29 words this package now encodes. There is no replacement: the levels a word actually drops is order - clip_order for its own decoded order, which :func:orders_of gives directly.

Parameters:

Name	Type	Description	Default
`clip_order`	`int`	HEALPix order to degrade to.	required
`midx`	`array-like of int`	Packed morton words (see :func:`res2display` for approximate resolutions).	required

Returns:

Type	Description
`ndarray`	Coarsened packed words, one per input word.

Source code in mortie/tools.py

def clip2order(clip_order, midx):
    """Coarsen packed morton words to a lower resolution.

    Degrades each packed word to ``clip_order`` by coarsening it through the
    kernel (the inverse of refining): the base cell and the first ``clip_order``
    tuples are kept, finer detail is dropped, and the suffix is rewritten. Words
    already at or below ``clip_order`` are returned unchanged.

    The ``print_factor`` flag was removed for the 1.x freeze (issue #68). It
    returned ``18 - clip_order``, a level count anchored to the retired
    decimal encoding's order-18 ceiling, so it went negative for the
    order-19..29 words this package now encodes. There is no replacement: the
    levels a word actually drops is ``order - clip_order`` for its own decoded
    order, which :func:`orders_of` gives directly.

    Parameters
    ----------
    clip_order : int
        HEALPix order to degrade to.
    midx : array-like of int
        Packed morton words (see :func:`res2display` for approximate resolutions).

    Returns
    -------
    ndarray
        Coarsened packed words, one per input word.
    """
    midx = np.ascontiguousarray(np.asarray(midx, dtype=np.uint64).ravel())
    return _rustie.rust_mi_coarsen(midx, int(clip_order))

`generate_morton_children(parent_morton, target_order)`

Generate all child morton indices at a target order.

Parameters:

Name	Type	Description	Default
`parent_morton`	`int`	Parent packed morton word.	required
`target_order`	`int`	Target order for children (must be >= parent order).	required

Returns:

Name	Type	Description
`children`	`ndarray`	Array of child packed morton words at target_order. If target_order equals parent_order, returns array with parent_morton.

Raises:

Type	Description
`ValueError`	If `target_order` is coarser than the parent word's own order.

Notes

Children are generated in HEALPix NESTED space — descending level_diff orders multiplies the cell count by 4**level_diff — then packed back to morton words via the kernel. If already at target_order, returns the parent itself.

Source code in mortie/tools.py

def generate_morton_children(parent_morton, target_order):
    """Generate all child morton indices at a target order.

    Parameters
    ----------
    parent_morton : int
        Parent packed morton word.
    target_order : int
        Target order for children (must be >= parent order).

    Returns
    -------
    children : ndarray
        Array of child packed morton words at target_order.
        If target_order equals parent_order, returns array with parent_morton.

    Raises
    ------
    ValueError
        If ``target_order`` is coarser than the parent word's own order.

    Notes
    -----
    Children are generated in HEALPix NESTED space — descending ``level_diff``
    orders multiplies the cell count by ``4**level_diff`` — then packed back to
    morton words via the kernel. If already at target_order, returns the parent
    itself.
    """
    # Decode the parent to its (nested, depth) via the packed kernel.
    parent_morton = np.uint64(parent_morton)
    nested, depths = _rust_mort2nested(
        np.ascontiguousarray(np.atleast_1d(parent_morton))
    )
    parent_order = int(depths[0])
    parent_nested = int(nested[0])

    if target_order < parent_order:
        raise ValueError(
            f"target_order ({target_order}) must be >= parent_order ({parent_order})"
        )

    if target_order == parent_order:
        return np.array([parent_morton], dtype=np.uint64)

    level_diff = target_order - parent_order
    # In NESTED space a cell's descendants at `target_order` are the contiguous
    # block `nested * 4**level_diff + [0 .. 4**level_diff)`.
    span = 4 ** level_diff
    child_nested = (parent_nested << (2 * level_diff)) + np.arange(
        span, dtype=np.uint64
    )
    depths = np.full(span, target_order, dtype=np.uint8)
    return _rust_nested2mort(np.ascontiguousarray(child_nested), depths)

`morton_buffer(morton_indices, k=1)`

Compute the k-cell border around a set of morton indices.

Returns only cells NOT in the input set (the expansion ring). User can union: np.union1d(morton_indices, border)

Parameters:

Name	Type	Description	Default
`morton_indices`	`array - like`	Morton indices, all at the same order.	required
`k`	`int`	Border width in cells (default 1, 8-connected neighbors). k=1 gives the immediate ring, k=2 gives a 2-cell border, etc.	`1`

Returns:

Name	Type	Description
`border`	`ndarray`	Sorted array of morton indices for the border cells.

Raises:

Type	Description
`ValueError`	If indices have mixed orders or k is out of range.

Source code in mortie/tools.py

def morton_buffer(morton_indices, k=1):
    """Compute the k-cell border around a set of morton indices.

    Returns only cells NOT in the input set (the expansion ring).
    User can union: ``np.union1d(morton_indices, border)``

    Parameters
    ----------
    morton_indices : array-like
        Morton indices, all at the same order.
    k : int, optional
        Border width in cells (default 1, 8-connected neighbors).
        k=1 gives the immediate ring, k=2 gives a 2-cell border, etc.

    Returns
    -------
    border : ndarray
        Sorted array of morton indices for the border cells.

    Raises
    ------
    ValueError
        If indices have mixed orders or k is out of range.
    """
    morton_indices = np.asarray(morton_indices, dtype=np.uint64)
    return _rustie.rust_morton_buffer(np.ascontiguousarray(morton_indices), k)

`morton_buffer_meters(morton_indices, width_m)`

Approximate meter-width buffer around a set of morton cells.

This is a convenience wrapper around :func:morton_buffer that picks k from the cells' HEALPix order so the resulting ring is roughly width_m meters wide. The input cells are assumed to all be at the same order.

.. warning:: This is an approximate buffer. The achieved width is rounded UP to the nearest whole HEALPix cell width — so the result always covers at least width_m meters, but may cover up to one cell width more. For order 18 cells (~30 m) the granularity is fine; at coarser orders it can be substantial. If you need a precise buffer, pick an order whose cell width is small relative to width_m and convert your input cells to that order first.

The cell width used for the calculation is the HEALPix angular resolution sqrt(pi/3) / nside converted to meters via the Earth's mean radius (6,371,008.77 m).

Parameters:

Name	Type	Description	Default
`morton_indices`	`array - like`	Morton indices, all at the same HEALPix order.	required
`width_m`	`float`	Desired buffer width in meters (must be > 0).	required

Returns:

Name	Type	Description
`border`	`ndarray`	Sorted array of morton indices for the border cells (NOT including the input cells). Union with the input if you want the filled ring: `np.union1d(morton_indices, border)`.

Raises:

Type	Description
`ValueError`	If `width_m` is non-positive, the input array is empty, or the cells are at mixed orders.

Examples:

>>> import mortie, numpy as np
>>> cells = mortie.linestring_coverage([10.0, 20.0], [30.0, 40.0], order=10)
>>> border = mortie.morton_buffer_meters(cells, width_m=5000.0)
>>> expanded = np.union1d(cells, border)

Source code in mortie/tools.py

def morton_buffer_meters(morton_indices, width_m):
    """Approximate meter-width buffer around a set of morton cells.

    This is a convenience wrapper around :func:`morton_buffer` that picks
    ``k`` from the cells' HEALPix order so the resulting ring is roughly
    *width_m* meters wide. The input cells are assumed to all be at the same
    order.

    .. warning::
       **This is an approximate buffer.** The achieved width is rounded
       UP to the nearest whole HEALPix cell width — so the result always
       covers *at least* ``width_m`` meters, but may cover up to one cell
       width more. For order 18 cells (~30 m) the granularity is fine; at
       coarser orders it can be substantial. If you need a precise buffer,
       pick an order whose cell width is small relative to ``width_m`` and
       convert your input cells to that order first.

    The cell width used for the calculation is the HEALPix angular
    resolution ``sqrt(pi/3) / nside`` converted to meters via the Earth's
    mean radius (6,371,008.77 m).

    Parameters
    ----------
    morton_indices : array-like
        Morton indices, all at the same HEALPix order.
    width_m : float
        Desired buffer width in meters (must be > 0).

    Returns
    -------
    border : ndarray
        Sorted array of morton indices for the border cells (NOT including
        the input cells). Union with the input if you want the filled ring:
        ``np.union1d(morton_indices, border)``.

    Raises
    ------
    ValueError
        If ``width_m`` is non-positive, the input array is empty, or the
        cells are at mixed orders.

    Examples
    --------
    >>> import mortie, numpy as np
    >>> cells = mortie.linestring_coverage([10.0, 20.0], [30.0, 40.0], order=10)
    >>> border = mortie.morton_buffer_meters(cells, width_m=5000.0)
    >>> expanded = np.union1d(cells, border)
    """
    morton_indices = np.asarray(morton_indices, dtype=np.uint64)
    if morton_indices.size == 0:
        raise ValueError("morton_indices must be non-empty")
    if not (width_m > 0):
        raise ValueError("width_m must be positive")

    # Infer order from the first cell. rust_morton_buffer will itself reject
    # mixed-order inputs downstream.
    order = infer_order_from_morton(int(morton_indices.flat[0]))
    if order < 1:
        raise ValueError("Could not infer a valid order from the input cells")

    nside = 1 << order
    cell_width_m = _EARTH_RADIUS_M * np.sqrt(np.pi / 3.0) / nside

    # Round UP so the buffer covers AT LEAST the requested width.
    k = int(np.ceil(width_m / cell_width_m))
    if k < 1:
        k = 1

    return morton_buffer(morton_indices, k=k)

`order2res(order)`

Approximate cell scale (km) at a HEALPix tessellation order.

The exact RMS cell spacing on the mean-radius HEALPix sphere: every order-k cell has identical area 4*pi*R**2 / (12 * 4**order) (HEALPix is equal-area), and the cell scale is the square root of that area. Derived from :data:EARTH_RADIUS_KM so code and the spec page (§3) share one Earth model (issue #119).

order may be a scalar (returns a float) or an array of orders such as :func:orders_of yields (returns an ndarray).

Parameters:

Name	Type	Description	Default
`order`	`int or array - like`	HEALPix tessellation order(s).	required

Returns:

Type	Description
`float or ndarray`	Approximate cell scale in kilometres (scalar in -> `float` out, array in -> `ndarray` out).

`res2display(max_order=MAX_ORDER)`

Resolution ladder for tessellation orders 0 through max_order.

Returns one record per order rather than printing (issue #68): each resolution is expressed in the largest sensible unit -- km at coarse orders, m once it drops below 1 km, cm once it drops below 1 m -- rounded to three decimals within that bracket, so fine orders read naturally (order 12 -> 1.592 km, order 13 -> 795.852 m) rather than as tiny km fractions.

Parameters:

Name	Type	Description	Default
`max_order`	`int`	Highest order to include, inclusive. Must lie in `0..MAX_ORDER` (default `MAX_ORDER` = 29, the finest order the packed-u64 kernel encodes).	`MAX_ORDER`

Returns:

Type	Description
`list of ResolutionLevel`	One named tuple `(order, value, unit, km)` per order, in ascending order. `value`/`unit` are the display pair; `km` is the unrounded resolution in kilometres for further arithmetic.

Raises:

Type	Description
`ValueError`	If `max_order` lies outside `0..MAX_ORDER`.