API reference

Every public function and class, generated from the source docstrings. The same names are importable directly from the top-level geoai3d package.

GEOAI_3D: geospatial-first AI workflows for 3D data.

This package is at a pre-alpha stage. No public API is stable yet; see the roadmap in the README for what is planned and in what order.

The first building blocks are the columnar :class:PointCloud representation and the :class:Provenance lineage record it carries.

AccuracyReport dataclass

A classification accuracy assessment.

Attributes:

Name	Type	Description
labels	tuple[int, ...]	The class labels, in the row/column order of confusion.
confusion	NDArray[Any]	The confusion matrix, rows = truth, columns = prediction.
overall_accuracy	float	Fraction of points classified correctly.
per_class	dict[int, dict[str, float]]	Per-label precision, recall, f1, iou, and support (the number of truth points of that class).
mean_iou	float	Mean of the per-class IoU values.
macro_f1	float	Mean of the per-class F1 values.

Source code in src/geoai3d/metrics.py

@dataclass(frozen=True)
class AccuracyReport:
    """A classification accuracy assessment.

    Attributes:
        labels: The class labels, in the row/column order of ``confusion``.
        confusion: The confusion matrix, rows = truth, columns = prediction.
        overall_accuracy: Fraction of points classified correctly.
        per_class: Per-label ``precision``, ``recall``, ``f1``, ``iou``, and
            ``support`` (the number of truth points of that class).
        mean_iou: Mean of the per-class IoU values.
        macro_f1: Mean of the per-class F1 values.
    """

    labels: tuple[int, ...]
    confusion: NDArray[Any]
    overall_accuracy: float
    per_class: dict[int, dict[str, float]]
    mean_iou: float
    macro_f1: float

Classifier dataclass

A trained classifier and the feature columns it expects.

Attributes:

Name	Type	Description
estimator	Any	The fitted scikit-learn estimator.
feature_names	tuple[str, ...]	The attribute names used as features, in order.

Source code in src/geoai3d/classify.py

@dataclass(frozen=True)
class Classifier:
    """A trained classifier and the feature columns it expects.

    Attributes:
        estimator: The fitted scikit-learn estimator.
        feature_names: The attribute names used as features, in order.
    """

    estimator: Any
    feature_names: tuple[str, ...]

Plane dataclass

A fitted plane a*x + b*y + c*z + d = 0 and its inliers.

Attributes:

Name	Type	Description
normal	tuple[float, float, float]	The unit normal (a, b, c).
offset	float	The plane offset d.
inliers	NDArray[Any]	Boolean mask over the fitted cloud, true for inlier points.

Source code in src/geoai3d/primitives.py

@dataclass(frozen=True)
class Plane:
    """A fitted plane ``a*x + b*y + c*z + d = 0`` and its inliers.

    Attributes:
        normal: The unit normal ``(a, b, c)``.
        offset: The plane offset ``d``.
        inliers: Boolean mask over the fitted cloud, true for inlier points.
    """

    normal: tuple[float, float, float]
    offset: float
    inliers: NDArray[Any]

    @property
    def num_inliers(self) -> int:
        """Return the number of inlier points."""
        return int(np.count_nonzero(self.inliers))

    def signed_distance(self, points: NDArray[Any]) -> NDArray[Any]:
        """Return the signed perpendicular distance of points to the plane.

        Args:
            points: An ``(n, 3)`` array of coordinates.

        Returns:
            An ``(n,)`` array of signed distances (the normal side is positive).
        """
        normal = np.asarray(self.normal, dtype=np.float64)
        return np.asarray(points, dtype=np.float64) @ normal + self.offset

num_inliers property

num_inliers: int

Return the number of inlier points.

signed_distance

signed_distance(points: NDArray[Any]) -> NDArray[Any]

Return the signed perpendicular distance of points to the plane.

Parameters:

Name	Type	Description	Default
points	NDArray[Any]	An (n, 3) array of coordinates.	required

Returns:

Type	Description
NDArray[Any]	An (n,) array of signed distances (the normal side is positive).

Source code in src/geoai3d/primitives.py

def signed_distance(self, points: NDArray[Any]) -> NDArray[Any]:
    """Return the signed perpendicular distance of points to the plane.

    Args:
        points: An ``(n, 3)`` array of coordinates.

    Returns:
        An ``(n,)`` array of signed distances (the normal side is positive).
    """
    normal = np.asarray(self.normal, dtype=np.float64)
    return np.asarray(points, dtype=np.float64) @ normal + self.offset

PointCloud

A georeferenced 3D point cloud stored column by column.

Parameters:

Name	Type	Description	Default
xyz	NDArray[Any]	Point coordinates as an array of shape (n_points, 3). Values are stored as float64; other numeric inputs are converted.	required
attributes	dict[str, NDArray[Any]] | None	Optional per-point attribute arrays keyed by name. Each array must be one-dimensional with length n_points (for example intensity or classification).	None
crs	object | None	Coordinate reference system the coordinates are expressed in. Optional at construction; the input/output and transform routines added in later stages require it and raise if it is missing. This slot will hold a pyproj CRS once the georeferencing layer lands.	None
provenance	Provenance | None	Lineage record describing how this cloud was produced.	None

Raises:

Type	Description
ValueError	If xyz is not a two-dimensional array with three columns, or if an attribute is not one-dimensional or does not match the number of points.

Example

import numpy as np from geoai3d import PointCloud xyz = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) cloud = PointCloud(xyz, attributes={"intensity": np.array([10, 20])}) len(cloud) 2 cloud.attribute_names ['intensity']

Source code in src/geoai3d/core/pointcloud.py

class PointCloud:
    """A georeferenced 3D point cloud stored column by column.

    Args:
        xyz: Point coordinates as an array of shape ``(n_points, 3)``. Values
            are stored as ``float64``; other numeric inputs are converted.
        attributes: Optional per-point attribute arrays keyed by name. Each
            array must be one-dimensional with length ``n_points`` (for
            example ``intensity`` or ``classification``).
        crs: Coordinate reference system the coordinates are expressed in.
            Optional at construction; the input/output and transform routines
            added in later stages require it and raise if it is missing. This
            slot will hold a ``pyproj`` CRS once the georeferencing layer lands.
        provenance: Lineage record describing how this cloud was produced.

    Raises:
        ValueError: If ``xyz`` is not a two-dimensional array with three
            columns, or if an attribute is not one-dimensional or does not
            match the number of points.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud
        >>> xyz = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
        >>> cloud = PointCloud(xyz, attributes={"intensity": np.array([10, 20])})
        >>> len(cloud)
        2
        >>> cloud.attribute_names
        ['intensity']
    """

    def __init__(
        self,
        xyz: NDArray[Any],
        attributes: dict[str, NDArray[Any]] | None = None,
        crs: object | None = None,
        provenance: Provenance | None = None,
    ) -> None:
        """Validate the inputs and store them in columnar form.

        The ``crs`` argument is coerced to a :class:`pyproj.CRS` so that
        :attr:`crs` always returns a resolved reference system (or ``None``),
        however the cloud was built.

        Raises:
            ValueError: If ``xyz`` is not ``(n_points, 3)``; if an attribute is
                not one-dimensional or its length does not match the point
                count; or if ``crs`` cannot be interpreted as a coordinate
                reference system.
        """
        coordinates: NDArray[np.float64] = np.ascontiguousarray(xyz, dtype=np.float64)
        if coordinates.ndim != 2 or coordinates.shape[1] != 3:
            msg = (
                "xyz must have shape (n_points, 3); "
                f"got an array with shape {coordinates.shape}."
            )
            raise ValueError(msg)
        self._xyz: NDArray[np.float64] = coordinates
        self._attributes: dict[str, NDArray[Any]] = {}
        if attributes is not None:
            for name, values in attributes.items():
                self._set_attribute_checked(name, values)
        self._crs = coerce_crs(crs) if crs is not None else None
        self._provenance = provenance

    def _set_attribute_checked(self, name: str, values: NDArray[Any]) -> None:
        """Validate an attribute array and store it under ``name``."""
        array: NDArray[Any] = np.asarray(values)
        if array.ndim != 1:
            msg = (
                f"Attribute {name!r} must be one-dimensional; "
                f"got an array with {array.ndim} dimensions."
            )
            raise ValueError(msg)
        if array.shape[0] != self._xyz.shape[0]:
            msg = (
                f"Attribute {name!r} has length {array.shape[0]}, but the cloud "
                f"has {self._xyz.shape[0]} points."
            )
            raise ValueError(msg)
        self._attributes[name] = array

    @property
    def xyz(self) -> NDArray[np.float64]:
        """Return the point coordinates as a ``(n_points, 3)`` ``float64`` array."""
        return self._xyz

    @property
    def crs(self) -> pyproj.CRS | None:
        """Return the coordinate reference system, or ``None`` if unset."""
        return self._crs

    @property
    def provenance(self) -> Provenance | None:
        """Return the lineage record, or ``None`` if unset."""
        return self._provenance

    @property
    def attribute_names(self) -> list[str]:
        """Return the attribute names, in insertion order."""
        return list(self._attributes)

    @property
    def bounds(self) -> tuple[float, float, float, float, float, float]:
        """Return axis-aligned bounds ``(minx, miny, minz, maxx, maxy, maxz)``.

        Raises:
            ValueError: If the cloud is empty.
        """
        if len(self) == 0:
            msg = "Cannot compute the bounds of an empty point cloud."
            raise ValueError(msg)
        mins = self._xyz.min(axis=0)
        maxs = self._xyz.max(axis=0)
        return (
            float(mins[0]),
            float(mins[1]),
            float(mins[2]),
            float(maxs[0]),
            float(maxs[1]),
            float(maxs[2]),
        )

    def attribute(self, name: str) -> NDArray[Any]:
        """Return the attribute array stored under ``name``.

        Args:
            name: Attribute name.

        Returns:
            The one-dimensional attribute array.

        Raises:
            KeyError: If no attribute with that name exists.
        """
        if name not in self._attributes:
            msg = (
                f"No attribute named {name!r}; "
                f"available attributes: {self.attribute_names}."
            )
            raise KeyError(msg)
        return self._attributes[name]

    def with_attribute(self, name: str, values: NDArray[Any]) -> PointCloud:
        """Return a new cloud with an attribute column added or replaced.

        The coordinate array is shared with this cloud rather than copied, so
        attaching a computed feature to a large cloud is inexpensive.

        Args:
            name: Attribute name to add or replace.
            values: One-dimensional array of length ``len(self)``.

        Returns:
            A new :class:`PointCloud` carrying the extra attribute. The original
            cloud is left unchanged.

        Raises:
            ValueError: If ``values`` is not one-dimensional or its length does
                not match the number of points.

        Example:
            >>> import numpy as np
            >>> from geoai3d import PointCloud
            >>> cloud = PointCloud(np.zeros((3, 3)))
            >>> planar = cloud.with_attribute("planarity", np.ones(3))
            >>> planar.attribute("planarity").tolist()
            [1.0, 1.0, 1.0]
            >>> cloud.attribute_names
            []
        """
        new = PointCloud(
            self._xyz,
            attributes=dict(self._attributes),
            crs=self._crs,
            provenance=self._copy_provenance(),
        )
        new._set_attribute_checked(name, values)
        return new

    def __len__(self) -> int:
        """Return the number of points in the cloud."""
        return int(self._xyz.shape[0])

    def __getitem__(self, selector: slice | NDArray[Any]) -> PointCloud:
        """Return a new cloud containing only the selected points.

        Args:
            selector: A slice, a boolean mask of length ``len(self)``, or an
                array of integer indices.

        Returns:
            A new :class:`PointCloud` with coordinates and every attribute
            filtered the same way, and the CRS and provenance carried through.

        Example:
            >>> import numpy as np
            >>> from geoai3d import PointCloud
            >>> cloud = PointCloud(np.arange(9.0).reshape(3, 3))
            >>> cloud[np.array([True, False, True])].xyz.shape
            (2, 3)
        """
        selected_attributes = {
            name: values[selector] for name, values in self._attributes.items()
        }
        return PointCloud(
            self._xyz[selector],
            attributes=selected_attributes,
            crs=self._crs,
            provenance=self._copy_provenance(),
        )

    def __repr__(self) -> str:
        """Return a concise developer-facing representation."""
        return (
            f"PointCloud(n_points={len(self)}, "
            f"attributes={self.attribute_names}, crs={self._crs!r})"
        )

    def _copy_provenance(self) -> Provenance | None:
        """Return an independent copy of the provenance, or ``None``."""
        return self._provenance.copy() if self._provenance is not None else None

attribute_names property

attribute_names: list[str]

Return the attribute names, in insertion order.

bounds property

bounds: tuple[float, float, float, float, float, float]

Return axis-aligned bounds (minx, miny, minz, maxx, maxy, maxz).

Raises:

Type	Description
ValueError	If the cloud is empty.

crs property

crs: CRS | None

Return the coordinate reference system, or None if unset.

provenance property

provenance: Provenance | None

Return the lineage record, or None if unset.

xyz property

xyz: NDArray[float64]

Return the point coordinates as a (n_points, 3) float64 array.

__getitem__

__getitem__(selector: slice | NDArray[Any]) -> PointCloud

Return a new cloud containing only the selected points.

Parameters:

Name	Type	Description	Default
selector	slice | NDArray[Any]	A slice, a boolean mask of length len(self), or an array of integer indices.	required

Returns:

Type	Description
PointCloud	A new :class:PointCloud with coordinates and every attribute
PointCloud	filtered the same way, and the CRS and provenance carried through.

Example

import numpy as np from geoai3d import PointCloud cloud = PointCloud(np.arange(9.0).reshape(3, 3)) cloud[np.array([True, False, True])].xyz.shape (2, 3)

Source code in src/geoai3d/core/pointcloud.py

def __getitem__(self, selector: slice | NDArray[Any]) -> PointCloud:
    """Return a new cloud containing only the selected points.

    Args:
        selector: A slice, a boolean mask of length ``len(self)``, or an
            array of integer indices.

    Returns:
        A new :class:`PointCloud` with coordinates and every attribute
        filtered the same way, and the CRS and provenance carried through.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud
        >>> cloud = PointCloud(np.arange(9.0).reshape(3, 3))
        >>> cloud[np.array([True, False, True])].xyz.shape
        (2, 3)
    """
    selected_attributes = {
        name: values[selector] for name, values in self._attributes.items()
    }
    return PointCloud(
        self._xyz[selector],
        attributes=selected_attributes,
        crs=self._crs,
        provenance=self._copy_provenance(),
    )

__init__

__init__(xyz: NDArray[Any], attributes: dict[str, NDArray[Any]] | None = None, crs: object | None = None, provenance: Provenance | None = None) -> None

Validate the inputs and store them in columnar form.

The crs argument is coerced to a :class:pyproj.CRS so that :attr:crs always returns a resolved reference system (or None), however the cloud was built.

Raises:

Type	Description
ValueError	If xyz is not (n_points, 3); if an attribute is not one-dimensional or its length does not match the point count; or if crs cannot be interpreted as a coordinate reference system.

Source code in src/geoai3d/core/pointcloud.py

def __init__(
    self,
    xyz: NDArray[Any],
    attributes: dict[str, NDArray[Any]] | None = None,
    crs: object | None = None,
    provenance: Provenance | None = None,
) -> None:
    """Validate the inputs and store them in columnar form.

    The ``crs`` argument is coerced to a :class:`pyproj.CRS` so that
    :attr:`crs` always returns a resolved reference system (or ``None``),
    however the cloud was built.

    Raises:
        ValueError: If ``xyz`` is not ``(n_points, 3)``; if an attribute is
            not one-dimensional or its length does not match the point
            count; or if ``crs`` cannot be interpreted as a coordinate
            reference system.
    """
    coordinates: NDArray[np.float64] = np.ascontiguousarray(xyz, dtype=np.float64)
    if coordinates.ndim != 2 or coordinates.shape[1] != 3:
        msg = (
            "xyz must have shape (n_points, 3); "
            f"got an array with shape {coordinates.shape}."
        )
        raise ValueError(msg)
    self._xyz: NDArray[np.float64] = coordinates
    self._attributes: dict[str, NDArray[Any]] = {}
    if attributes is not None:
        for name, values in attributes.items():
            self._set_attribute_checked(name, values)
    self._crs = coerce_crs(crs) if crs is not None else None
    self._provenance = provenance

__len__

__len__() -> int

Return the number of points in the cloud.

Source code in src/geoai3d/core/pointcloud.py

def __len__(self) -> int:
    """Return the number of points in the cloud."""
    return int(self._xyz.shape[0])

__repr__

__repr__() -> str

Return a concise developer-facing representation.

Source code in src/geoai3d/core/pointcloud.py

def __repr__(self) -> str:
    """Return a concise developer-facing representation."""
    return (
        f"PointCloud(n_points={len(self)}, "
        f"attributes={self.attribute_names}, crs={self._crs!r})"
    )

attribute

attribute(name: str) -> NDArray[Any]

Return the attribute array stored under name.

Parameters:

Name	Type	Description	Default
name	str	Attribute name.	required

Returns:

Type	Description
NDArray[Any]	The one-dimensional attribute array.

Raises:

Type	Description
KeyError	If no attribute with that name exists.

Source code in src/geoai3d/core/pointcloud.py

def attribute(self, name: str) -> NDArray[Any]:
    """Return the attribute array stored under ``name``.

    Args:
        name: Attribute name.

    Returns:
        The one-dimensional attribute array.

    Raises:
        KeyError: If no attribute with that name exists.
    """
    if name not in self._attributes:
        msg = (
            f"No attribute named {name!r}; "
            f"available attributes: {self.attribute_names}."
        )
        raise KeyError(msg)
    return self._attributes[name]

with_attribute

with_attribute(name: str, values: NDArray[Any]) -> PointCloud

Return a new cloud with an attribute column added or replaced.

The coordinate array is shared with this cloud rather than copied, so attaching a computed feature to a large cloud is inexpensive.

Parameters:

Name	Type	Description	Default
name	str	Attribute name to add or replace.	required
values	NDArray[Any]	One-dimensional array of length len(self).	required

Returns:

Type	Description
PointCloud	A new :class:PointCloud carrying the extra attribute. The original
PointCloud	cloud is left unchanged.

Raises:

Type	Description
ValueError	If values is not one-dimensional or its length does not match the number of points.

Example

import numpy as np from geoai3d import PointCloud cloud = PointCloud(np.zeros((3, 3))) planar = cloud.with_attribute("planarity", np.ones(3)) planar.attribute("planarity").tolist() [1.0, 1.0, 1.0] cloud.attribute_names []

Source code in src/geoai3d/core/pointcloud.py

def with_attribute(self, name: str, values: NDArray[Any]) -> PointCloud:
    """Return a new cloud with an attribute column added or replaced.

    The coordinate array is shared with this cloud rather than copied, so
    attaching a computed feature to a large cloud is inexpensive.

    Args:
        name: Attribute name to add or replace.
        values: One-dimensional array of length ``len(self)``.

    Returns:
        A new :class:`PointCloud` carrying the extra attribute. The original
        cloud is left unchanged.

    Raises:
        ValueError: If ``values`` is not one-dimensional or its length does
            not match the number of points.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud
        >>> cloud = PointCloud(np.zeros((3, 3)))
        >>> planar = cloud.with_attribute("planarity", np.ones(3))
        >>> planar.attribute("planarity").tolist()
        [1.0, 1.0, 1.0]
        >>> cloud.attribute_names
        []
    """
    new = PointCloud(
        self._xyz,
        attributes=dict(self._attributes),
        crs=self._crs,
        provenance=self._copy_provenance(),
    )
    new._set_attribute_checked(name, values)
    return new

ProcessStep dataclass

A single recorded step in a data product's processing history.

Attributes:

Name	Type	Description
description	str	Human-readable summary of what the step did, for example "voxel subsampling".
parameters	dict[str, Any]	The parameters that controlled the step, as a plain mapping of names to values, for example {"voxel_size": 0.5}.
software	str	Name and version of the software that ran the step, for example "geoai3d 0.1.0".
timestamp	datetime	When the step ran, as a timezone-aware UTC datetime.

Source code in src/geoai3d/core/provenance.py

@dataclass(frozen=True)
class ProcessStep:
    """A single recorded step in a data product's processing history.

    Attributes:
        description: Human-readable summary of what the step did, for example
            ``"voxel subsampling"``.
        parameters: The parameters that controlled the step, as a plain mapping
            of names to values, for example ``{"voxel_size": 0.5}``.
        software: Name and version of the software that ran the step, for
            example ``"geoai3d 0.1.0"``.
        timestamp: When the step ran, as a timezone-aware UTC datetime.
    """

    description: str
    parameters: dict[str, Any] = field(default_factory=dict)
    software: str = field(default_factory=_default_software)
    timestamp: datetime = field(default_factory=_utc_now)

    def to_dict(self) -> dict[str, Any]:
        """Return a JSON-serialisable representation of this step.

        Returns:
            A mapping with the description, parameters, software, and an ISO
            8601 timestamp string.
        """
        return {
            "description": self.description,
            "parameters": dict(self.parameters),
            "software": self.software,
            "timestamp": self.timestamp.isoformat(),
        }

    @classmethod
    def from_dict(cls, data: dict[str, Any]) -> ProcessStep:
        """Rebuild a :class:`ProcessStep` from :meth:`to_dict` output.

        Args:
            data: A mapping as produced by :meth:`to_dict`.

        Returns:
            The reconstructed :class:`ProcessStep`.
        """
        return cls(
            description=data["description"],
            parameters=dict(data.get("parameters", {})),
            software=data["software"],
            timestamp=datetime.fromisoformat(data["timestamp"]),
        )

from_dict classmethod

from_dict(data: dict[str, Any]) -> ProcessStep

Rebuild a :class:ProcessStep from :meth:to_dict output.

Parameters:

Name	Type	Description	Default
data	dict[str, Any]	A mapping as produced by :meth:to_dict.	required

Returns:

Type	Description
ProcessStep	The reconstructed :class:ProcessStep.

Source code in src/geoai3d/core/provenance.py

@classmethod
def from_dict(cls, data: dict[str, Any]) -> ProcessStep:
    """Rebuild a :class:`ProcessStep` from :meth:`to_dict` output.

    Args:
        data: A mapping as produced by :meth:`to_dict`.

    Returns:
        The reconstructed :class:`ProcessStep`.
    """
    return cls(
        description=data["description"],
        parameters=dict(data.get("parameters", {})),
        software=data["software"],
        timestamp=datetime.fromisoformat(data["timestamp"]),
    )

to_dict

to_dict() -> dict[str, Any]

Return a JSON-serialisable representation of this step.

Returns:

Type	Description
dict[str, Any]	A mapping with the description, parameters, software, and an ISO
dict[str, Any]	8601 timestamp string.

Source code in src/geoai3d/core/provenance.py

def to_dict(self) -> dict[str, Any]:
    """Return a JSON-serialisable representation of this step.

    Returns:
        A mapping with the description, parameters, software, and an ISO
        8601 timestamp string.
    """
    return {
        "description": self.description,
        "parameters": dict(self.parameters),
        "software": self.software,
        "timestamp": self.timestamp.isoformat(),
    }

Provenance dataclass

An ordered lineage record for a data product.

A :class:Provenance bundles an optional description of the original data source with an ordered list of :class:ProcessStep entries, oldest first.

Parameters:

Name	Type	Description	Default
source	str | None	Description of the original input, for example a file name or a dataset DOI. None if unknown.	None
steps	list[ProcessStep]	The processing steps applied so far, oldest first.	list()

Example

prov = Provenance(source="scan.laz") _ = prov.add_step("voxel subsampling", {"voxel_size": 0.5}) len(prov.steps) 1

Source code in src/geoai3d/core/provenance.py

@dataclass
class Provenance:
    """An ordered lineage record for a data product.

    A :class:`Provenance` bundles an optional description of the original data
    source with an ordered list of :class:`ProcessStep` entries, oldest first.

    Args:
        source: Description of the original input, for example a file name or a
            dataset DOI. ``None`` if unknown.
        steps: The processing steps applied so far, oldest first.

    Example:
        >>> prov = Provenance(source="scan.laz")
        >>> _ = prov.add_step("voxel subsampling", {"voxel_size": 0.5})
        >>> len(prov.steps)
        1
    """

    source: str | None = None
    steps: list[ProcessStep] = field(default_factory=list)

    def add_step(
        self,
        description: str,
        parameters: dict[str, Any] | None = None,
        software: str | None = None,
    ) -> ProcessStep:
        """Append a processing step to the lineage and return it.

        Args:
            description: Human-readable summary of the step.
            parameters: Parameters that controlled the step. A copy is stored,
                so later mutation of the passed mapping does not affect the
                record. Defaults to an empty mapping.
            software: Software name and version. Defaults to the running
                ``geoai3d`` version.

        Returns:
            The :class:`ProcessStep` that was appended.

        Example:
            >>> prov = Provenance()
            >>> step = prov.add_step("read LAS", {"path": "scan.laz"})
            >>> step.description
            'read LAS'
        """
        step = ProcessStep(
            description=description,
            parameters=dict(parameters) if parameters is not None else {},
            software=software if software is not None else _default_software(),
        )
        self.steps.append(step)
        return step

    def copy(self) -> Provenance:
        """Return an independent copy of this lineage record.

        Returns:
            A new :class:`Provenance` whose step list can be extended without
            affecting the original. Individual :class:`ProcessStep` entries are
            immutable and therefore safe to share.
        """
        return Provenance(source=self.source, steps=list(self.steps))

    def to_dict(self) -> dict[str, Any]:
        """Return a JSON-serialisable representation of this lineage record.

        Returns:
            A mapping with the source and a list of serialised steps, suitable
            for storing in file metadata or a STAC item.
        """
        return {
            "source": self.source,
            "steps": [step.to_dict() for step in self.steps],
        }

    @classmethod
    def from_dict(cls, data: dict[str, Any]) -> Provenance:
        """Rebuild a :class:`Provenance` from :meth:`to_dict` output.

        Args:
            data: A mapping as produced by :meth:`to_dict`.

        Returns:
            The reconstructed :class:`Provenance`.
        """
        steps = [ProcessStep.from_dict(entry) for entry in data.get("steps", [])]
        return cls(source=data.get("source"), steps=steps)

add_step

add_step(description: str, parameters: dict[str, Any] | None = None, software: str | None = None) -> ProcessStep

Append a processing step to the lineage and return it.

Parameters:

Name	Type	Description	Default
description	str	Human-readable summary of the step.	required
parameters	dict[str, Any] | None	Parameters that controlled the step. A copy is stored, so later mutation of the passed mapping does not affect the record. Defaults to an empty mapping.	None
software	str | None	Software name and version. Defaults to the running geoai3d version.	None

Returns:

Name	Type	Description
The	ProcessStep	class:ProcessStep that was appended.

Example

prov = Provenance() step = prov.add_step("read LAS", {"path": "scan.laz"}) step.description 'read LAS'

Source code in src/geoai3d/core/provenance.py

def add_step(
    self,
    description: str,
    parameters: dict[str, Any] | None = None,
    software: str | None = None,
) -> ProcessStep:
    """Append a processing step to the lineage and return it.

    Args:
        description: Human-readable summary of the step.
        parameters: Parameters that controlled the step. A copy is stored,
            so later mutation of the passed mapping does not affect the
            record. Defaults to an empty mapping.
        software: Software name and version. Defaults to the running
            ``geoai3d`` version.

    Returns:
        The :class:`ProcessStep` that was appended.

    Example:
        >>> prov = Provenance()
        >>> step = prov.add_step("read LAS", {"path": "scan.laz"})
        >>> step.description
        'read LAS'
    """
    step = ProcessStep(
        description=description,
        parameters=dict(parameters) if parameters is not None else {},
        software=software if software is not None else _default_software(),
    )
    self.steps.append(step)
    return step

copy

copy() -> Provenance

Return an independent copy of this lineage record.

Returns:

Type	Description
Provenance	A new :class:Provenance whose step list can be extended without
Provenance	affecting the original. Individual :class:ProcessStep entries are
Provenance	immutable and therefore safe to share.

Source code in src/geoai3d/core/provenance.py

def copy(self) -> Provenance:
    """Return an independent copy of this lineage record.

    Returns:
        A new :class:`Provenance` whose step list can be extended without
        affecting the original. Individual :class:`ProcessStep` entries are
        immutable and therefore safe to share.
    """
    return Provenance(source=self.source, steps=list(self.steps))

from_dict classmethod

from_dict(data: dict[str, Any]) -> Provenance

Rebuild a :class:Provenance from :meth:to_dict output.

Parameters:

Name	Type	Description	Default
data	dict[str, Any]	A mapping as produced by :meth:to_dict.	required

Returns:

Type	Description
Provenance	The reconstructed :class:Provenance.

Source code in src/geoai3d/core/provenance.py

@classmethod
def from_dict(cls, data: dict[str, Any]) -> Provenance:
    """Rebuild a :class:`Provenance` from :meth:`to_dict` output.

    Args:
        data: A mapping as produced by :meth:`to_dict`.

    Returns:
        The reconstructed :class:`Provenance`.
    """
    steps = [ProcessStep.from_dict(entry) for entry in data.get("steps", [])]
    return cls(source=data.get("source"), steps=steps)

to_dict

to_dict() -> dict[str, Any]

Return a JSON-serialisable representation of this lineage record.

Returns:

Type	Description
dict[str, Any]	A mapping with the source and a list of serialised steps, suitable
dict[str, Any]	for storing in file metadata or a STAC item.

Source code in src/geoai3d/core/provenance.py

def to_dict(self) -> dict[str, Any]:
    """Return a JSON-serialisable representation of this lineage record.

    Returns:
        A mapping with the source and a list of serialised steps, suitable
        for storing in file metadata or a STAC item.
    """
    return {
        "source": self.source,
        "steps": [step.to_dict() for step in self.steps],
    }

Raster

A 2D georeferenced raster: values with an affine transform and CRS.

Source code in src/geoai3d/core/raster.py

class Raster:
    """A 2D georeferenced raster: values with an affine transform and CRS."""

    def __init__(
        self,
        data: NDArray[Any],
        transform: tuple[float, ...],
        crs: object,
        nodata: float = float("nan"),
    ) -> None:
        """Validate and store the grid, transform, CRS, and nodata value.

        Args:
            data: A 2D array of cell values, shape ``(rows, cols)``.
            transform: The affine transform as a 6-tuple ``(a, b, c, d, e, f)``.
            crs: The coordinate reference system (EPSG code, WKT, or
                ``pyproj.CRS``); coerced to a :class:`pyproj.CRS`.
            nodata: Value marking empty cells. Defaults to ``NaN``.

        Raises:
            ValueError: If ``data`` is not 2D, ``transform`` is not length 6,
                or ``crs`` is ``None`` or uninterpretable.
        """
        array: NDArray[Any] = np.ascontiguousarray(data)
        if array.ndim != 2:
            msg = f"data must be 2D (rows, cols); got shape {array.shape}."
            raise ValueError(msg)
        if len(tuple(transform)) != 6:
            msg = "transform must be a 6-tuple (a, b, c, d, e, f)."
            raise ValueError(msg)
        if crs is None:
            msg = (
                "A raster requires a coordinate reference system; pass crs= "
                "(an EPSG code, WKT string, or pyproj.CRS)."
            )
            raise ValueError(msg)
        self._data = array
        self._transform: _Transform = (
            float(transform[0]),
            float(transform[1]),
            float(transform[2]),
            float(transform[3]),
            float(transform[4]),
            float(transform[5]),
        )
        self._crs = coerce_crs(crs)
        self._nodata = float(nodata)

    @property
    def data(self) -> NDArray[Any]:
        """Return the cell values as a 2D ``(rows, cols)`` array."""
        return self._data

    @property
    def transform(self) -> _Transform:
        """Return the affine transform as ``(a, b, c, d, e, f)``."""
        return self._transform

    @property
    def crs(self) -> pyproj.CRS:
        """Return the coordinate reference system."""
        return self._crs

    @property
    def nodata(self) -> float:
        """Return the value marking empty cells."""
        return self._nodata

    @property
    def shape(self) -> tuple[int, int]:
        """Return the grid shape as ``(rows, cols)``."""
        return (int(self._data.shape[0]), int(self._data.shape[1]))

    @property
    def height(self) -> int:
        """Return the number of rows."""
        return int(self._data.shape[0])

    @property
    def width(self) -> int:
        """Return the number of columns."""
        return int(self._data.shape[1])

    @property
    def resolution(self) -> tuple[float, float]:
        """Return the ``(x, y)`` cell size in CRS units (both positive)."""
        return (abs(self._transform[0]), abs(self._transform[4]))

    @property
    def bounds(self) -> tuple[float, float, float, float]:
        """Return axis-aligned bounds ``(min_x, min_y, max_x, max_y)``."""
        a, _, c, _, e, f = self._transform
        min_x = c
        max_y = f
        max_x = c + self.width * a
        min_y = f + self.height * e
        return (
            float(min(min_x, max_x)),
            float(min(min_y, max_y)),
            float(max(min_x, max_x)),
            float(max(min_y, max_y)),
        )

    def __repr__(self) -> str:
        """Return a concise, informative representation."""
        epsg = self._crs.to_epsg()
        crs_label = f"EPSG:{epsg}" if epsg is not None else "CRS"
        res_x, res_y = self.resolution
        return (
            f"Raster(shape=({self.height}, {self.width}), "
            f"resolution=({res_x:g}, {res_y:g}), {crs_label})"
        )

bounds property

bounds: tuple[float, float, float, float]

Return axis-aligned bounds (min_x, min_y, max_x, max_y).

crs property

crs: CRS

Return the coordinate reference system.

data property

data: NDArray[Any]

Return the cell values as a 2D (rows, cols) array.

height property

height: int

Return the number of rows.

nodata property

nodata: float

Return the value marking empty cells.

resolution property

resolution: tuple[float, float]

Return the (x, y) cell size in CRS units (both positive).

shape property

shape: tuple[int, int]

Return the grid shape as (rows, cols).

transform property

transform: _Transform

Return the affine transform as (a, b, c, d, e, f).

width property

width: int

Return the number of columns.

__init__

__init__(data: NDArray[Any], transform: tuple[float, ...], crs: object, nodata: float = float('nan')) -> None

Validate and store the grid, transform, CRS, and nodata value.

Parameters:

Name	Type	Description	Default
data	NDArray[Any]	A 2D array of cell values, shape (rows, cols).	required
transform	tuple[float, ...]	The affine transform as a 6-tuple (a, b, c, d, e, f).	required
crs	object	The coordinate reference system (EPSG code, WKT, or pyproj.CRS); coerced to a :class:pyproj.CRS.	required
nodata	float	Value marking empty cells. Defaults to NaN.	float('nan')

Raises:

Type	Description
ValueError	If data is not 2D, transform is not length 6, or crs is None or uninterpretable.

Source code in src/geoai3d/core/raster.py

def __init__(
    self,
    data: NDArray[Any],
    transform: tuple[float, ...],
    crs: object,
    nodata: float = float("nan"),
) -> None:
    """Validate and store the grid, transform, CRS, and nodata value.

    Args:
        data: A 2D array of cell values, shape ``(rows, cols)``.
        transform: The affine transform as a 6-tuple ``(a, b, c, d, e, f)``.
        crs: The coordinate reference system (EPSG code, WKT, or
            ``pyproj.CRS``); coerced to a :class:`pyproj.CRS`.
        nodata: Value marking empty cells. Defaults to ``NaN``.

    Raises:
        ValueError: If ``data`` is not 2D, ``transform`` is not length 6,
            or ``crs`` is ``None`` or uninterpretable.
    """
    array: NDArray[Any] = np.ascontiguousarray(data)
    if array.ndim != 2:
        msg = f"data must be 2D (rows, cols); got shape {array.shape}."
        raise ValueError(msg)
    if len(tuple(transform)) != 6:
        msg = "transform must be a 6-tuple (a, b, c, d, e, f)."
        raise ValueError(msg)
    if crs is None:
        msg = (
            "A raster requires a coordinate reference system; pass crs= "
            "(an EPSG code, WKT string, or pyproj.CRS)."
        )
        raise ValueError(msg)
    self._data = array
    self._transform: _Transform = (
        float(transform[0]),
        float(transform[1]),
        float(transform[2]),
        float(transform[3]),
        float(transform[4]),
        float(transform[5]),
    )
    self._crs = coerce_crs(crs)
    self._nodata = float(nodata)

__repr__

__repr__() -> str

Return a concise, informative representation.

Source code in src/geoai3d/core/raster.py

def __repr__(self) -> str:
    """Return a concise, informative representation."""
    epsg = self._crs.to_epsg()
    crs_label = f"EPSG:{epsg}" if epsg is not None else "CRS"
    res_x, res_y = self.resolution
    return (
        f"Raster(shape=({self.height}, {self.width}), "
        f"resolution=({res_x:g}, {res_y:g}), {crs_label})"
    )

build_kdtree

build_kdtree(cloud: PointCloud) -> cKDTree

Build a KD-tree over a cloud's coordinates.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to index.	required

Returns:

Type	Description
cKDTree	A scipy.spatial.cKDTree over the (n_points, 3) coordinates.

Example

import numpy as np from geoai3d import PointCloud, build_kdtree cloud = PointCloud(np.arange(300.0).reshape(100, 3), crs=28992) tree = build_kdtree(cloud) distances, indices = tree.query([0.0, 0.0, 0.0], k=3) indices.shape (3,)

Source code in src/geoai3d/index.py

def build_kdtree(cloud: PointCloud) -> cKDTree:
    """Build a KD-tree over a cloud's coordinates.

    Args:
        cloud: The cloud to index.

    Returns:
        A ``scipy.spatial.cKDTree`` over the ``(n_points, 3)`` coordinates.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, build_kdtree
        >>> cloud = PointCloud(np.arange(300.0).reshape(100, 3), crs=28992)
        >>> tree = build_kdtree(cloud)
        >>> distances, indices = tree.query([0.0, 0.0, 0.0], k=3)
        >>> indices.shape
        (3,)
    """
    return cKDTree(cloud.xyz)

confusion_matrix

confusion_matrix(truth: NDArray[Any], prediction: NDArray[Any], labels: NDArray[Any] | list[int] | tuple[int, ...] | None = None) -> tuple[NDArray[Any], tuple[int, ...]]

Build a confusion matrix (rows = truth, columns = prediction).

Parameters:

Name	Type	Description	Default
truth	NDArray[Any]	True integer labels, one per point.	required
prediction	NDArray[Any]	Predicted integer labels, one per point.	required
labels	NDArray[Any] | list[int] | tuple[int, ...] | None	The labels to include, in order. If omitted, the sorted union of the values present is used. Points whose truth or prediction falls outside labels are ignored.	None

Returns:

Type	Description
NDArray[Any]	The (k, k) integer confusion matrix and the tuple of labels giving
tuple[int, ...]	its row and column order.

Raises:

Type	Description
ValueError	If truth and prediction have different lengths, or are empty.

Source code in src/geoai3d/metrics.py

def confusion_matrix(
    truth: NDArray[Any],
    prediction: NDArray[Any],
    labels: NDArray[Any] | list[int] | tuple[int, ...] | None = None,
) -> tuple[NDArray[Any], tuple[int, ...]]:
    """Build a confusion matrix (rows = truth, columns = prediction).

    Args:
        truth: True integer labels, one per point.
        prediction: Predicted integer labels, one per point.
        labels: The labels to include, in order. If omitted, the sorted union of
            the values present is used. Points whose truth or prediction falls
            outside ``labels`` are ignored.

    Returns:
        The ``(k, k)`` integer confusion matrix and the tuple of labels giving
        its row and column order.

    Raises:
        ValueError: If ``truth`` and ``prediction`` have different lengths, or
            are empty.
    """
    true_labels = np.asarray(truth).ravel()
    predicted_labels = np.asarray(prediction).ravel()
    if true_labels.shape != predicted_labels.shape:
        msg = (
            f"truth and prediction must have the same length; got "
            f"{true_labels.shape[0]} and {predicted_labels.shape[0]}."
        )
        raise ValueError(msg)
    if true_labels.size == 0:
        msg = "Cannot assess accuracy on empty label arrays."
        raise ValueError(msg)

    if labels is None:
        label_values = np.unique(np.concatenate([true_labels, predicted_labels]))
    else:
        label_values = np.array(sorted({int(label) for label in labels}))
    n_labels = len(label_values)

    keep = np.isin(true_labels, label_values) & np.isin(predicted_labels, label_values)
    row = np.searchsorted(label_values, true_labels[keep])
    column = np.searchsorted(label_values, predicted_labels[keep])
    matrix = np.zeros((n_labels, n_labels), dtype=np.int64)
    np.add.at(matrix, (row, column), 1)
    return matrix, tuple(int(label) for label in label_values)

connected_components

connected_components(cloud: PointCloud, *, distance: float, min_size: int = 1, output_attribute: str = 'component') -> PointCloud

Label points into connected components by a distance threshold.

Two points are connected if they lie within distance of each other, and each connected group gets a distinct integer label. Components with fewer than min_size points are labelled -1 (noise); the rest are numbered 0, 1, 2, ... in order of first appearance.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to label.	required
distance	float	Maximum distance, in CRS units, between two points for them to be connected.	required
min_size	int	Minimum number of points for a component to be kept; smaller groups are labelled -1.	1
output_attribute	str	Name of the integer label column to add.	'component'

Returns:

Type	Description
PointCloud	A new :class:~geoai3d.PointCloud with the integer component label
PointCloud	column added, the CRS and provenance carried through.

Raises:

Type	Description
ValueError	If the cloud is empty, distance is not positive, or min_size is below 1.

Example

import numpy as np from geoai3d import PointCloud, connected_components a = np.random.default_rng(0).normal(0.0, 0.1, (50, 3)) b = np.random.default_rng(1).normal(10.0, 0.1, (50, 3)) cloud = connected_components( ... PointCloud(np.vstack([a, b]), crs=28992), distance=0.5 ... ) len(set(cloud.attribute("component").tolist())) 2

Source code in src/geoai3d/cluster.py

def connected_components(
    cloud: PointCloud,
    *,
    distance: float,
    min_size: int = 1,
    output_attribute: str = "component",
) -> PointCloud:
    """Label points into connected components by a distance threshold.

    Two points are connected if they lie within ``distance`` of each other, and
    each connected group gets a distinct integer label. Components with fewer
    than ``min_size`` points are labelled ``-1`` (noise); the rest are numbered
    ``0, 1, 2, ...`` in order of first appearance.

    Args:
        cloud: The cloud to label.
        distance: Maximum distance, in CRS units, between two points for them to
            be connected.
        min_size: Minimum number of points for a component to be kept; smaller
            groups are labelled ``-1``.
        output_attribute: Name of the integer label column to add.

    Returns:
        A new :class:`~geoai3d.PointCloud` with the integer component label
        column added, the CRS and provenance carried through.

    Raises:
        ValueError: If the cloud is empty, ``distance`` is not positive, or
            ``min_size`` is below 1.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, connected_components
        >>> a = np.random.default_rng(0).normal(0.0, 0.1, (50, 3))
        >>> b = np.random.default_rng(1).normal(10.0, 0.1, (50, 3))
        >>> cloud = connected_components(
        ...     PointCloud(np.vstack([a, b]), crs=28992), distance=0.5
        ... )
        >>> len(set(cloud.attribute("component").tolist()))
        2
    """
    n_points = len(cloud)
    if n_points == 0:
        msg = "Cannot cluster an empty cloud."
        raise ValueError(msg)
    if distance <= 0.0:
        msg = "distance must be a positive number."
        raise ValueError(msg)
    if min_size < 1:
        msg = "min_size must be at least 1."
        raise ValueError(msg)

    tree = cKDTree(cloud.xyz)
    pairs = tree.query_pairs(distance, output_type="ndarray")
    if pairs.size:
        edge_weights = np.ones(len(pairs))
        graph = coo_matrix(
            (edge_weights, (pairs[:, 0], pairs[:, 1])),
            shape=(n_points, n_points),
        )
        _, raw_labels = _sparse_components(graph, directed=False)
    else:
        raw_labels = np.arange(n_points, dtype=np.int64)

    counts = np.bincount(raw_labels)
    keep = counts >= min_size
    relabel = np.full(len(counts), -1, dtype=np.int64)
    relabel[np.flatnonzero(keep)] = np.arange(int(np.count_nonzero(keep)))
    component = relabel[raw_labels]

    return attach_attributes(
        cloud,
        {output_attribute: component},
        "connected_components",
        {"distance": distance, "min_size": min_size},
    )

dbscan

dbscan(cloud: PointCloud, *, eps: float, min_samples: int = 5, output_attribute: str = 'cluster') -> PointCloud

Cluster points by density with DBSCAN.

DBSCAN groups points that have at least min_samples neighbours within eps, growing clusters through those dense points and leaving sparse points unclustered. Unlike :func:connected_components, a thin bridge of isolated points does not merge two clusters.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to cluster.	required
eps	float	Neighbourhood radius, in CRS units.	required
min_samples	int	Minimum neighbours (including the point itself) for a point to be a dense core of a cluster.	5
output_attribute	str	Name of the integer label column to add. Cluster labels are 0, 1, 2, ... and unclustered points are -1.	'cluster'

Returns:

Type	Description
PointCloud	A new :class:~geoai3d.PointCloud with the integer cluster label column
PointCloud	added, the CRS and provenance carried through.

Raises:

Type	Description
ValueError	If the cloud is empty, eps is not positive, or min_samples is below 1.

Example

import numpy as np from geoai3d import PointCloud, dbscan a = np.random.default_rng(0).normal(0.0, 0.1, (50, 3)) b = np.random.default_rng(1).normal(10.0, 0.1, (50, 3)) out = dbscan(PointCloud(np.vstack([a, b]), crs=28992), eps=0.5) int(out.attribute("cluster").max()) 1

Source code in src/geoai3d/cluster.py

def dbscan(
    cloud: PointCloud,
    *,
    eps: float,
    min_samples: int = 5,
    output_attribute: str = "cluster",
) -> PointCloud:
    """Cluster points by density with DBSCAN.

    DBSCAN groups points that have at least ``min_samples`` neighbours within
    ``eps``, growing clusters through those dense points and leaving sparse
    points unclustered. Unlike :func:`connected_components`, a thin bridge of
    isolated points does not merge two clusters.

    Args:
        cloud: The cloud to cluster.
        eps: Neighbourhood radius, in CRS units.
        min_samples: Minimum neighbours (including the point itself) for a point
            to be a dense core of a cluster.
        output_attribute: Name of the integer label column to add. Cluster
            labels are ``0, 1, 2, ...`` and unclustered points are ``-1``.

    Returns:
        A new :class:`~geoai3d.PointCloud` with the integer cluster label column
        added, the CRS and provenance carried through.

    Raises:
        ValueError: If the cloud is empty, ``eps`` is not positive, or
            ``min_samples`` is below 1.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, dbscan
        >>> a = np.random.default_rng(0).normal(0.0, 0.1, (50, 3))
        >>> b = np.random.default_rng(1).normal(10.0, 0.1, (50, 3))
        >>> out = dbscan(PointCloud(np.vstack([a, b]), crs=28992), eps=0.5)
        >>> int(out.attribute("cluster").max())
        1
    """
    if len(cloud) == 0:
        msg = "Cannot cluster an empty cloud."
        raise ValueError(msg)
    if eps <= 0.0:
        msg = "eps must be a positive number."
        raise ValueError(msg)
    if min_samples < 1:
        msg = "min_samples must be at least 1."
        raise ValueError(msg)

    from sklearn.cluster import DBSCAN

    labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(cloud.xyz)
    return attach_attributes(
        cloud,
        {output_attribute: labels.astype(np.int64)},
        "dbscan",
        {"eps": eps, "min_samples": min_samples},
    )

difference

difference(minuend: Raster, subtrahend: Raster) -> Raster

Subtract one aligned raster from another, cell by cell.

difference(dsm, dtm) gives a normalised surface (object heights); difference(later, earlier) gives elevation change between epochs. A cell is nodata (NaN) where either input is nodata.

Parameters:

Name	Type	Description	Default
minuend	Raster	The raster to subtract from.	required
subtrahend	Raster	The raster to subtract.	required

Returns:

Name	Type	Description
A	Raster	class:~geoai3d.Raster of minuend - subtrahend.

Raises:

Type	Description
ValueError	If the rasters differ in shape, affine transform, or CRS (rasterise both on the same bounds and resolution first).

Example

import numpy as np from geoai3d import PointCloud, to_dtm, to_dsm, difference pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3)) cloud = PointCloud(pts, crs=28992) ndsm = difference(to_dsm(cloud, resolution=1.0), ... to_dtm(cloud, resolution=1.0)) ndsm.shape (10, 10)

Source code in src/geoai3d/dem.py

def difference(minuend: Raster, subtrahend: Raster) -> Raster:
    """Subtract one aligned raster from another, cell by cell.

    ``difference(dsm, dtm)`` gives a normalised surface (object heights);
    ``difference(later, earlier)`` gives elevation change between epochs. A cell
    is nodata (``NaN``) where either input is nodata.

    Args:
        minuend: The raster to subtract from.
        subtrahend: The raster to subtract.

    Returns:
        A :class:`~geoai3d.Raster` of ``minuend - subtrahend``.

    Raises:
        ValueError: If the rasters differ in shape, affine transform, or CRS
            (rasterise both on the same ``bounds`` and ``resolution`` first).

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_dtm, to_dsm, difference
        >>> pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3))
        >>> cloud = PointCloud(pts, crs=28992)
        >>> ndsm = difference(to_dsm(cloud, resolution=1.0),
        ...                   to_dtm(cloud, resolution=1.0))
        >>> ndsm.shape
        (10, 10)
    """
    if minuend.shape != subtrahend.shape:
        msg = f"Rasters are not aligned: shapes {minuend.shape} vs {subtrahend.shape}."
        raise ValueError(msg)
    if not np.allclose(minuend.transform, subtrahend.transform):
        msg = (
            "Rasters are not aligned: their affine transforms differ. "
            "Rasterise both on the same bounds and resolution."
        )
        raise ValueError(msg)
    if minuend.crs != subtrahend.crs:
        msg = "Rasters have different CRS; reproject to a common CRS first."
        raise ValueError(msg)
    result = _finite(minuend) - _finite(subtrahend)
    return Raster(result, minuend.transform, minuend.crs, nodata=float("nan"))

estimate_normals

estimate_normals(cloud: PointCloud, *, k: int = 16, orient_upward: bool = True) -> PointCloud

Estimate a surface normal at each point from its neighbours.

The normal is the direction of least variance of the k nearest neighbours (the smallest principal component of the local covariance).

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to add normals to.	required
k	int	Number of nearest neighbours used to fit each local plane. Must be at least 3.	16
orient_upward	bool	If true, flip every normal to point into the upper hemisphere (positive z), giving a consistent sign.	True

Returns:

Type	Description
PointCloud	A new :class:PointCloud with nx, ny, nz attributes added.

Raises:

Type	Description
ValueError	If k is less than 3.

Example

import numpy as np from geoai3d import PointCloud, estimate_normals flat = np.random.default_rng(0).random((200, 3)) flat[:, 2] = 0.0 cloud = estimate_normals(PointCloud(flat, crs=28992), k=8) bool(np.allclose(np.abs(cloud.attribute("nz")), 1.0)) True

Source code in src/geoai3d/geometry.py

def estimate_normals(
    cloud: PointCloud,
    *,
    k: int = 16,
    orient_upward: bool = True,
) -> PointCloud:
    """Estimate a surface normal at each point from its neighbours.

    The normal is the direction of least variance of the ``k`` nearest
    neighbours (the smallest principal component of the local covariance).

    Args:
        cloud: The cloud to add normals to.
        k: Number of nearest neighbours used to fit each local plane. Must be
            at least 3.
        orient_upward: If true, flip every normal to point into the upper
            hemisphere (positive z), giving a consistent sign.

    Returns:
        A new :class:`PointCloud` with ``nx``, ``ny``, ``nz`` attributes added.

    Raises:
        ValueError: If ``k`` is less than 3.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, estimate_normals
        >>> flat = np.random.default_rng(0).random((200, 3))
        >>> flat[:, 2] = 0.0
        >>> cloud = estimate_normals(PointCloud(flat, crs=28992), k=8)
        >>> bool(np.allclose(np.abs(cloud.attribute("nz")), 1.0))
        True
    """
    if k < 3:
        msg = "k must be at least 3 to fit a local plane."
        raise ValueError(msg)
    xyz = cloud.xyz
    k_effective = min(k, len(cloud))
    tree = cKDTree(xyz)
    _, neighbor_indices = tree.query(xyz, k=k_effective, workers=-1)
    neighbors = xyz[neighbor_indices]
    centered = neighbors - neighbors.mean(axis=1, keepdims=True)
    covariance = np.einsum("nki,nkj->nij", centered, centered) / k_effective
    _, vectors = np.linalg.eigh(covariance)
    normals = vectors[:, :, 0]
    if orient_upward:
        downward = normals[:, 2] < 0
        normals[downward] = -normals[downward]
    return attach_attributes(
        cloud,
        {
            "nx": np.ascontiguousarray(normals[:, 0]),
            "ny": np.ascontiguousarray(normals[:, 1]),
            "nz": np.ascontiguousarray(normals[:, 2]),
        },
        "estimate_normals",
        {"k": k, "orient_upward": orient_upward},
    )

estimate_radius

estimate_radius(cloud: PointCloud, *, target_neighbors: int = 20, sample_size: int = 10000, seed: int = 0) -> float

Estimate a feature radius that captures about target_neighbors points.

Measures, over a random sample of points, the distance to the target_neighbors-th nearest neighbour and returns the median. That is a radius which, on average, encloses roughly that many neighbours given the cloud's own density -- so the caller does not have to guess a number.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to measure.	required
target_neighbors	int	Desired average neighbour count within the radius.	20
sample_size	int	Number of points to sample for the estimate.	10000
seed	int	Seed for the random sample.	0

Returns:

Type	Description
float	The estimated radius in the cloud's CRS units.

Raises:

Type	Description
ValueError	If the cloud has fewer than two points, or target_neighbors is below 1.

Example

import numpy as np from geoai3d import PointCloud, estimate_radius cloud = PointCloud(np.random.default_rng(0).random((2000, 3)), crs=28992) radius = estimate_radius(cloud, target_neighbors=15) radius > 0 True

Source code in src/geoai3d/outofcore.py

def estimate_radius(
    cloud: PointCloud,
    *,
    target_neighbors: int = 20,
    sample_size: int = 10000,
    seed: int = 0,
) -> float:
    """Estimate a feature radius that captures about ``target_neighbors`` points.

    Measures, over a random sample of points, the distance to the
    ``target_neighbors``-th nearest neighbour and returns the median. That is a
    radius which, on average, encloses roughly that many neighbours given the
    cloud's own density -- so the caller does not have to guess a number.

    Args:
        cloud: The cloud to measure.
        target_neighbors: Desired average neighbour count within the radius.
        sample_size: Number of points to sample for the estimate.
        seed: Seed for the random sample.

    Returns:
        The estimated radius in the cloud's CRS units.

    Raises:
        ValueError: If the cloud has fewer than two points, or
            ``target_neighbors`` is below 1.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, estimate_radius
        >>> cloud = PointCloud(np.random.default_rng(0).random((2000, 3)), crs=28992)
        >>> radius = estimate_radius(cloud, target_neighbors=15)
        >>> radius > 0
        True
    """
    if target_neighbors < 1:
        msg = "target_neighbors must be at least 1."
        raise ValueError(msg)
    n_points = len(cloud)
    if n_points < 2:
        msg = "estimate_radius needs at least two points."
        raise ValueError(msg)
    xyz = cloud.xyz
    tree = cKDTree(xyz)
    k = min(target_neighbors + 1, n_points)
    rng = np.random.default_rng(seed)
    if n_points > sample_size:
        sample = xyz[rng.choice(n_points, size=sample_size, replace=False)]
    else:
        sample = xyz
    distances, _ = tree.query(sample, k=k, workers=-1)
    return float(np.median(distances[:, -1]))

evaluate

evaluate(truth: NDArray[Any], prediction: NDArray[Any], *, labels: NDArray[Any] | list[int] | tuple[int, ...] | None = None) -> AccuracyReport

Assess predicted labels against ground truth.

Parameters:

Name	Type	Description	Default
truth	NDArray[Any]	True integer labels, one per point.	required
prediction	NDArray[Any]	Predicted integer labels, one per point.	required
labels	NDArray[Any] | list[int] | tuple[int, ...] | None	The labels to include, in order. If omitted, the sorted union of the values present is used.	None

Returns:

Name	Type	Description
An	AccuracyReport	class:AccuracyReport with the confusion matrix, overall accuracy,
	AccuracyReport	per-class metrics, mean IoU, and macro F1.

Raises:

Type	Description
ValueError	If truth and prediction have different lengths, or are empty.

Example

import numpy as np from geoai3d import evaluate truth = np.array([0, 0, 1, 1, 1]) prediction = np.array([0, 1, 1, 1, 1]) report = evaluate(truth, prediction) round(report.overall_accuracy, 2) 0.8

Source code in src/geoai3d/metrics.py

def evaluate(
    truth: NDArray[Any],
    prediction: NDArray[Any],
    *,
    labels: NDArray[Any] | list[int] | tuple[int, ...] | None = None,
) -> AccuracyReport:
    """Assess predicted labels against ground truth.

    Args:
        truth: True integer labels, one per point.
        prediction: Predicted integer labels, one per point.
        labels: The labels to include, in order. If omitted, the sorted union of
            the values present is used.

    Returns:
        An :class:`AccuracyReport` with the confusion matrix, overall accuracy,
        per-class metrics, mean IoU, and macro F1.

    Raises:
        ValueError: If ``truth`` and ``prediction`` have different lengths, or
            are empty.

    Example:
        >>> import numpy as np
        >>> from geoai3d import evaluate
        >>> truth = np.array([0, 0, 1, 1, 1])
        >>> prediction = np.array([0, 1, 1, 1, 1])
        >>> report = evaluate(truth, prediction)
        >>> round(report.overall_accuracy, 2)
        0.8
    """
    matrix, label_values = confusion_matrix(truth, prediction, labels)
    total = int(matrix.sum())
    overall_accuracy = float(np.trace(matrix)) / total if total > 0 else 0.0

    per_class: dict[int, dict[str, float]] = {}
    iou_values: list[float] = []
    f1_values: list[float] = []
    for index, label in enumerate(label_values):
        true_positive = float(matrix[index, index])
        false_positive = float(matrix[:, index].sum()) - true_positive
        false_negative = float(matrix[index, :].sum()) - true_positive
        precision = (
            true_positive / (true_positive + false_positive)
            if true_positive + false_positive > 0
            else 0.0
        )
        recall = (
            true_positive / (true_positive + false_negative)
            if true_positive + false_negative > 0
            else 0.0
        )
        f1 = (
            2.0 * precision * recall / (precision + recall)
            if precision + recall > 0
            else 0.0
        )
        union = true_positive + false_positive + false_negative
        iou = true_positive / union if union > 0 else 0.0
        per_class[label] = {
            "precision": precision,
            "recall": recall,
            "f1": f1,
            "iou": iou,
            "support": float(matrix[index, :].sum()),
        }
        iou_values.append(iou)
        f1_values.append(f1)

    return AccuracyReport(
        labels=label_values,
        confusion=matrix,
        overall_accuracy=overall_accuracy,
        per_class=per_class,
        mean_iou=float(np.mean(iou_values)) if iou_values else 0.0,
        macro_f1=float(np.mean(f1_values)) if f1_values else 0.0,
    )

feature_matrix

feature_matrix(cloud: PointCloud, feature_names: tuple[str, ...] | list[str] | None = None) -> NDArray[Any]

Assemble a (n_points, n_features) matrix from named attributes.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to read features from.	required
feature_names	tuple[str, ...] | list[str] | None	Attribute names to stack, in order. If omitted, the geometric feature columns (and normalized_intensity) present on the cloud are used.	None

Returns:

Type	Description
NDArray[Any]	A float64 feature matrix.

Raises:

Type	Description
ValueError	If no features can be resolved, or a named attribute is absent.

Example

import numpy as np from geoai3d import PointCloud, feature_matrix cloud = PointCloud( ... np.zeros((4, 3)), ... attributes={"planarity": np.ones(4), "linearity": np.zeros(4)}, ... crs=28992, ... ) feature_matrix(cloud, ["planarity", "linearity"]).shape (4, 2)

Source code in src/geoai3d/classify.py

def feature_matrix(
    cloud: PointCloud,
    feature_names: tuple[str, ...] | list[str] | None = None,
) -> NDArray[Any]:
    """Assemble a ``(n_points, n_features)`` matrix from named attributes.

    Args:
        cloud: The cloud to read features from.
        feature_names: Attribute names to stack, in order. If omitted, the
            geometric feature columns (and ``normalized_intensity``) present on
            the cloud are used.

    Returns:
        A ``float64`` feature matrix.

    Raises:
        ValueError: If no features can be resolved, or a named attribute is
            absent.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, feature_matrix
        >>> cloud = PointCloud(
        ...     np.zeros((4, 3)),
        ...     attributes={"planarity": np.ones(4), "linearity": np.zeros(4)},
        ...     crs=28992,
        ... )
        >>> feature_matrix(cloud, ["planarity", "linearity"]).shape
        (4, 2)
    """
    names = _resolve_features(cloud, feature_names)
    columns = [cloud.attribute(name).astype(np.float64) for name in names]
    return np.asarray(np.column_stack(columns), dtype=np.float64)

fit_plane

fit_plane(cloud: PointCloud, *, distance_threshold: float = 0.1, max_iterations: int = 1000, seed: int = 0) -> Plane

Fit the dominant plane to a cloud with RANSAC.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to fit. Needs at least three points.	required
distance_threshold	float	Maximum perpendicular distance, in CRS units, for a point to count as an inlier.	0.1
max_iterations	int	Number of random three-point samples to try.	1000
seed	int	Seed for the random sampling, so the fit is reproducible.	0

Returns:

Type	Description
Plane	The fitted :class:Plane, refit to its inliers.

Raises:

Type	Description
ValueError	If the cloud has fewer than three points, or distance_threshold or max_iterations is not positive, or no plane could be fitted (all samples were collinear).

Example

import numpy as np from geoai3d import PointCloud, fit_plane rng = np.random.default_rng(0) xy = rng.uniform(0.0, 10.0, (500, 2)) z = 0.0 * xy[:, 0] # a horizontal plane at z = 0 pts = np.column_stack([xy, z]) plane = fit_plane(PointCloud(pts, crs=28992), distance_threshold=0.01) bool(abs(abs(plane.normal[2]) - 1.0) < 1e-6) True

Source code in src/geoai3d/primitives.py

def fit_plane(
    cloud: PointCloud,
    *,
    distance_threshold: float = 0.1,
    max_iterations: int = 1000,
    seed: int = 0,
) -> Plane:
    """Fit the dominant plane to a cloud with RANSAC.

    Args:
        cloud: The cloud to fit. Needs at least three points.
        distance_threshold: Maximum perpendicular distance, in CRS units, for a
            point to count as an inlier.
        max_iterations: Number of random three-point samples to try.
        seed: Seed for the random sampling, so the fit is reproducible.

    Returns:
        The fitted :class:`Plane`, refit to its inliers.

    Raises:
        ValueError: If the cloud has fewer than three points, or
            ``distance_threshold`` or ``max_iterations`` is not positive, or no
            plane could be fitted (all samples were collinear).

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, fit_plane
        >>> rng = np.random.default_rng(0)
        >>> xy = rng.uniform(0.0, 10.0, (500, 2))
        >>> z = 0.0 * xy[:, 0]  # a horizontal plane at z = 0
        >>> pts = np.column_stack([xy, z])
        >>> plane = fit_plane(PointCloud(pts, crs=28992), distance_threshold=0.01)
        >>> bool(abs(abs(plane.normal[2]) - 1.0) < 1e-6)
        True
    """
    n_points = len(cloud)
    if n_points < 3:
        msg = "fit_plane needs at least three points."
        raise ValueError(msg)
    if distance_threshold <= 0.0:
        msg = "distance_threshold must be a positive number."
        raise ValueError(msg)
    if max_iterations < 1:
        msg = "max_iterations must be at least 1."
        raise ValueError(msg)

    xyz = cloud.xyz
    rng = np.random.default_rng(seed)
    best_inliers: NDArray[Any] | None = None
    best_count = -1
    for _ in range(max_iterations):
        sample = xyz[rng.choice(n_points, size=3, replace=False)]
        normal = np.cross(sample[1] - sample[0], sample[2] - sample[0])
        norm = float(np.linalg.norm(normal))
        if norm == 0.0:  # collinear sample, no plane
            continue
        normal = normal / norm
        offset = -float(normal @ sample[0])
        inliers = np.abs(xyz @ normal + offset) < distance_threshold
        count = int(np.count_nonzero(inliers))
        if count > best_count:
            best_count = count
            best_inliers = inliers

    if best_inliers is None:
        msg = "Could not fit a plane: every sample was collinear."
        raise ValueError(msg)

    # Refit to all inliers by least squares: the normal is the direction of
    # least spread (the smallest right-singular vector of the centred inliers).
    inlier_points = xyz[best_inliers]
    centroid = inlier_points.mean(axis=0)
    _, _, right_vectors = np.linalg.svd(inlier_points - centroid)
    normal = right_vectors[-1]
    offset = -float(normal @ centroid)
    inliers = np.abs(xyz @ normal + offset) < distance_threshold
    return Plane(
        normal=(float(normal[0]), float(normal[1]), float(normal[2])),
        offset=offset,
        inliers=inliers,
    )

geometric_features

geometric_features(cloud: PointCloud, *, k: int | None = None, radius: float | None = None) -> PointCloud

Add eigenvalue-based geometric feature columns to a cloud.

Provide exactly one of k (a fixed number of neighbours) or radius (a fixed spatial scale). If neither is given, k=20 is used. Only the radius form carries the bit-identical out-of-core seam guarantee.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to describe.	required
k	int | None	Number of nearest neighbours per point. Must be at least 3.	None
radius	float | None	Neighbourhood radius in the cloud's CRS units. Points with fewer than three neighbours within it get NaN features.	None

Returns:

Type	Description
PointCloud	A new :class:PointCloud with the feature columns (linearity,
PointCloud	planarity, sphericity, anisotropy, omnivariance,
PointCloud	eigenentropy, surface_variation, verticality,
PointCloud	sum_eigenvalues) added.

Raises:

Type	Description
ValueError	If both k and radius are given, or k is below 3, or radius is not positive.

Example

import numpy as np from geoai3d import PointCloud, geometric_features grid = np.linspace(-1, 1, 21) gx, gy = np.meshgrid(grid, grid) flat = np.column_stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)]) cloud = geometric_features(PointCloud(flat, crs=28992), k=13) bool(np.nanmedian(cloud.attribute("planarity")) > 0.9) True

Source code in src/geoai3d/features.py

def geometric_features(
    cloud: PointCloud,
    *,
    k: int | None = None,
    radius: float | None = None,
) -> PointCloud:
    """Add eigenvalue-based geometric feature columns to a cloud.

    Provide exactly one of ``k`` (a fixed number of neighbours) or ``radius``
    (a fixed spatial scale). If neither is given, ``k=20`` is used. Only the
    ``radius`` form carries the bit-identical out-of-core seam guarantee.

    Args:
        cloud: The cloud to describe.
        k: Number of nearest neighbours per point. Must be at least 3.
        radius: Neighbourhood radius in the cloud's CRS units. Points with
            fewer than three neighbours within it get NaN features.

    Returns:
        A new :class:`PointCloud` with the feature columns (``linearity``,
        ``planarity``, ``sphericity``, ``anisotropy``, ``omnivariance``,
        ``eigenentropy``, ``surface_variation``, ``verticality``,
        ``sum_eigenvalues``) added.

    Raises:
        ValueError: If both ``k`` and ``radius`` are given, or ``k`` is below 3,
            or ``radius`` is not positive.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, geometric_features
        >>> grid = np.linspace(-1, 1, 21)
        >>> gx, gy = np.meshgrid(grid, grid)
        >>> flat = np.column_stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)])
        >>> cloud = geometric_features(PointCloud(flat, crs=28992), k=13)
        >>> bool(np.nanmedian(cloud.attribute("planarity")) > 0.9)
        True
    """
    if k is not None and radius is not None:
        msg = "Give either k= or radius=, not both."
        raise ValueError(msg)
    xyz = cloud.xyz
    tree = cKDTree(xyz)
    if radius is not None:
        if radius <= 0:
            msg = "radius must be a positive number."
            raise ValueError(msg)
        positions = np.arange(len(cloud), dtype=np.intp)
        eigenvalues, normals = radius_eigen(xyz, positions, tree, positions, radius)
    else:
        neighbors = _DEFAULT_K if k is None else k
        if neighbors < 3:
            msg = "k must be at least 3 to fit a local neighbourhood."
            raise ValueError(msg)
        _, neighbor_indices = tree.query(xyz, k=min(neighbors, len(cloud)), workers=-1)
        eigenvalues, normals = knn_eigen(xyz[neighbor_indices])

    features = features_from_eigen(eigenvalues, normals)
    return attach_attributes(
        cloud, features, "geometric_features", {"k": k, "radius": radius}
    )

geometric_features_stream

geometric_features_stream(source: str | PathLike[str], destination: str | PathLike[str], *, radius: float, tile_size: float, crs: object | None = None, attributes: Sequence[str] | None = None, chunk_size: int = 1000000, workers: int = 1, verbose: bool = False, max_points: int | None = None) -> None

Compute fixed-radius geometric features on a LAS/LAZ file out-of-core.

Streams the file in two passes (partition to per-tile temp files, then per-tile feature computation) so the whole cloud is never resident. The output matches :func:~geoai3d.geometric_features with the same radius, bit for bit, and carries every source point attribute through unchanged.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to the input .las or .laz file.	required
destination	str | PathLike[str]	Path to the output .parquet file.	required
radius	float	Neighbourhood radius and halo width, in the file's CRS units.	required
tile_size	float	Edge length of the square (x, y) tiles. Choose it so the number of tiles is modest (dozens to a few hundred); it must be at least radius.	required
crs	object | None	CRS to assign, overriding the file header. Required if the file has none.	None
attributes	Sequence[str] | None	Source point dimensions to carry through to the output. None (default) carries every dimension the file stores (intensity, returns, classification, GPS time, and so on). Pass a sequence of dimension names to carry only those and keep the output smaller.	None
chunk_size	int	Number of points read per streaming chunk in pass one.	1000000
workers	int	Number of processes for the per-tile feature pass. 1 runs serially.	1
verbose	bool	If true, print the time spent in the partition pass and the feature pass, and the tile count.	False
max_points	int | None	If set, process only about this many points (rounded up to a chunk). Mainly for benchmarking a slice of a very large file.	None

Raises:

Type	Description
ValueError	If radius or tile_size is not positive, tile_size is smaller than radius, workers is below 1, an entry in attributes is not a dimension of the source file, or the file has no CRS and none is supplied.

Source code in src/geoai3d/stream.py

def geometric_features_stream(
    source: str | os.PathLike[str],
    destination: str | os.PathLike[str],
    *,
    radius: float,
    tile_size: float,
    crs: object | None = None,
    attributes: Sequence[str] | None = None,
    chunk_size: int = 1_000_000,
    workers: int = 1,
    verbose: bool = False,
    max_points: int | None = None,
) -> None:
    """Compute fixed-radius geometric features on a LAS/LAZ file out-of-core.

    Streams the file in two passes (partition to per-tile temp files, then
    per-tile feature computation) so the whole cloud is never resident. The
    output matches :func:`~geoai3d.geometric_features` with the same radius,
    bit for bit, and carries every source point attribute through unchanged.

    Args:
        source: Path to the input ``.las`` or ``.laz`` file.
        destination: Path to the output ``.parquet`` file.
        radius: Neighbourhood radius and halo width, in the file's CRS units.
        tile_size: Edge length of the square (x, y) tiles. Choose it so the
            number of tiles is modest (dozens to a few hundred); it must be at
            least ``radius``.
        crs: CRS to assign, overriding the file header. Required if the file
            has none.
        attributes: Source point dimensions to carry through to the output.
            ``None`` (default) carries every dimension the file stores
            (intensity, returns, classification, GPS time, and so on). Pass a
            sequence of dimension names to carry only those and keep the output
            smaller.
        chunk_size: Number of points read per streaming chunk in pass one.
        workers: Number of processes for the per-tile feature pass. 1 runs
            serially.
        verbose: If true, print the time spent in the partition pass and the
            feature pass, and the tile count.
        max_points: If set, process only about this many points (rounded up to
            a chunk). Mainly for benchmarking a slice of a very large file.

    Raises:
        ValueError: If ``radius`` or ``tile_size`` is not positive,
            ``tile_size`` is smaller than ``radius``, ``workers`` is below 1,
            an entry in ``attributes`` is not a dimension of the source file, or
            the file has no CRS and none is supplied.
    """
    if radius <= 0:
        msg = "radius must be a positive number."
        raise ValueError(msg)
    if tile_size <= 0:
        msg = "tile_size must be a positive number."
        raise ValueError(msg)
    if tile_size < radius:
        msg = "tile_size must be at least radius, or halos would span whole tiles."
        raise ValueError(msg)
    if workers < 1:
        msg = "workers must be at least 1."
        raise ValueError(msg)

    with laspy.open(str(source)) as reader:
        header = reader.header
        file_crs = header.parse_crs()
        resolved_crs = coerce_crs(crs) if crs is not None else file_crs
        if resolved_crs is None:
            msg = (
                f"{str(source)!r} has no coordinate reference system and none "
                "was supplied. Pass crs= to set it explicitly."
            )
            raise ValueError(msg)
        attribute_names = _attribute_dimensions(header.point_format, attributes)
        origin_x = float(header.mins[0])
        origin_y = float(header.mins[1])

        temp_dir = Path(tempfile.mkdtemp(prefix="geoai3d_tiles_"))
        try:
            partition_start = time.perf_counter()
            paths, temp_schema = _partition(
                reader,
                temp_dir,
                origin_x,
                origin_y,
                tile_size,
                radius,
                chunk_size,
                attribute_names,
                max_points,
            )
            partition_seconds = time.perf_counter() - partition_start
            if verbose:
                print(
                    f"[geoai3d] partition pass: {partition_seconds:.1f}s "
                    f"into {len(paths)} tiles",
                    flush=True,
                )
            feature_start = time.perf_counter()
            out_schema = pa.schema(
                [
                    *temp_schema,
                    *[pa.field(name, pa.float64()) for name in FEATURE_NAMES],
                ]
            ).with_metadata(
                _output_metadata(
                    resolved_crs, str(source), radius, tile_size, attribute_names
                )
            )
            items = [
                (str(path), tile_i, tile_j, origin_x, origin_y, tile_size, radius)
                for (tile_i, tile_j), path in paths.items()
            ]
            writer = pq.ParquetWriter(str(destination), out_schema)
            try:
                if workers > 1:
                    with ProcessPoolExecutor(max_workers=workers) as pool:
                        for result in pool.map(_tile_worker, items):
                            _write_result(writer, out_schema, result)
                else:
                    for item in items:
                        _write_result(writer, out_schema, _tile_worker(item))
            finally:
                writer.close()
            if verbose:
                feature_seconds = time.perf_counter() - feature_start
                print(
                    f"[geoai3d] feature pass: {feature_seconds:.1f}s "
                    f"({workers} worker(s))",
                    flush=True,
                )
        finally:
            shutil.rmtree(temp_dir, ignore_errors=True)

geometric_features_tiled

geometric_features_tiled(cloud: PointCloud, *, radius: float, tile_size: float) -> PointCloud

Compute fixed-radius geometric features tile by tile with a halo.

The result is identical to :func:geometric_features with the same radius on the whole cloud -- bit for bit -- but each tile is processed with a bounded working set (the tile plus a one-radius halo), so memory stays low for clouds much larger than RAM.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to describe.	required
radius	float	Neighbourhood radius. Points with fewer than three neighbours within it get NaN features. Also the halo width.	required
tile_size	float	Edge length of the square (x, y) tiles, in CRS units.	required

Returns:

Type	Description
PointCloud	A new :class:PointCloud with the geometric feature columns added.

Raises:

Type	Description
ValueError	If radius or tile_size is not positive.

Example

import numpy as np from geoai3d import PointCloud, geometric_features_tiled pts = np.random.default_rng(0).uniform(0, 10, (2000, 3)) cloud = geometric_features_tiled(PointCloud(pts, crs=28992), ... radius=0.6, tile_size=2.5) "planarity" in cloud.attribute_names True

Source code in src/geoai3d/outofcore.py

def geometric_features_tiled(
    cloud: PointCloud,
    *,
    radius: float,
    tile_size: float,
) -> PointCloud:
    """Compute fixed-radius geometric features tile by tile with a halo.

    The result is identical to :func:`geometric_features` with the same
    ``radius`` on the whole cloud -- bit for bit -- but each tile is processed
    with a bounded working set (the tile plus a one-radius halo), so memory
    stays low for clouds much larger than RAM.

    Args:
        cloud: The cloud to describe.
        radius: Neighbourhood radius. Points with fewer than three neighbours
            within it get NaN features. Also the halo width.
        tile_size: Edge length of the square (x, y) tiles, in CRS units.

    Returns:
        A new :class:`PointCloud` with the geometric feature columns added.

    Raises:
        ValueError: If ``radius`` or ``tile_size`` is not positive.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, geometric_features_tiled
        >>> pts = np.random.default_rng(0).uniform(0, 10, (2000, 3))
        >>> cloud = geometric_features_tiled(PointCloud(pts, crs=28992),
        ...                                  radius=0.6, tile_size=2.5)
        >>> "planarity" in cloud.attribute_names
        True
    """
    if radius <= 0:
        msg = "radius must be a positive number."
        raise ValueError(msg)
    if tile_size <= 0:
        msg = "tile_size must be a positive number."
        raise ValueError(msg)

    xyz = cloud.xyz
    n_points = len(cloud)
    global_index = np.arange(n_points, dtype=np.intp)
    origin = xyz[:, :2].min(axis=0)
    tile_ij = np.floor((xyz[:, :2] - origin) / tile_size).astype(np.int64)

    results: dict[str, NDArray[Any]] = {
        name: np.full(n_points, np.nan) for name in FEATURE_NAMES
    }
    for tile_i, tile_j in np.unique(tile_ij, axis=0):
        core_mask = (tile_ij[:, 0] == tile_i) & (tile_ij[:, 1] == tile_j)
        x_min = origin[0] + tile_i * tile_size
        x_max = x_min + tile_size
        y_min = origin[1] + tile_j * tile_size
        y_max = y_min + tile_size
        pool_mask = (
            (xyz[:, 0] >= x_min - radius)
            & (xyz[:, 0] < x_max + radius)
            & (xyz[:, 1] >= y_min - radius)
            & (xyz[:, 1] < y_max + radius)
        )
        pool_positions = np.flatnonzero(pool_mask)
        pool_xyz = xyz[pool_positions]
        pool_global = global_index[pool_positions]
        tree = cKDTree(pool_xyz)

        core_global = np.flatnonzero(core_mask)
        core_in_pool = np.searchsorted(pool_positions, core_global)
        eigenvalues, normals = radius_eigen(
            pool_xyz, pool_global, tree, core_in_pool, radius
        )
        tile_features = features_from_eigen(eigenvalues, normals)
        for name in FEATURE_NAMES:
            results[name][core_global] = tile_features[name]

    return attach_attributes(
        cloud,
        results,
        "geometric_features_tiled",
        {"radius": radius, "tile_size": tile_size},
    )

load_classifier

load_classifier(source: str | PathLike[str]) -> Classifier

Load a classifier saved by :func:save_classifier.

Only load models from a source you trust: the file is unpickled.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to a saved classifier.	required

Returns:

Type	Description
Classifier	The loaded :class:Classifier.

Raises:

Type	Description
ValueError	If the file does not contain a GEOAI_3D classifier.

Source code in src/geoai3d/classify.py

def load_classifier(source: str | os.PathLike[str]) -> Classifier:
    """Load a classifier saved by :func:`save_classifier`.

    Only load models from a source you trust: the file is unpickled.

    Args:
        source: Path to a saved classifier.

    Returns:
        The loaded :class:`Classifier`.

    Raises:
        ValueError: If the file does not contain a GEOAI_3D classifier.
    """
    with open(source, "rb") as handle:
        loaded = pickle.load(handle)
    if not isinstance(loaded, Classifier):
        msg = f"{str(source)!r} does not contain a GEOAI_3D Classifier."
        raise ValueError(msg)
    return loaded

multiscale_features

multiscale_features(cloud: PointCloud, *, ks: tuple[int, ...] = (10, 20, 30, 40, 50)) -> PointCloud

Add geometric features computed at each point's optimal neighbourhood.

Features are computed at several neighbourhood sizes and, per point, the size that minimises eigenentropy is chosen (the dimensionality-based scale selection of Demantke et al.). The chosen size is stored as optimal_k.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to describe.	required
ks	tuple[int, ...]	Candidate neighbour counts to try. Each must be at least 3.	(10, 20, 30, 40, 50)

Returns:

Type	Description
PointCloud	A new :class:PointCloud with the geometric feature columns at the
PointCloud	per-point optimal scale, plus an optimal_k column.

Raises:

Type	Description
ValueError	If any candidate in ks is below 3, or ks is empty.

Example

import numpy as np from geoai3d import PointCloud, multiscale_features pts = np.random.default_rng(0).random((400, 3)) cloud = multiscale_features(PointCloud(pts, crs=28992), ks=(8, 16, 24)) "optimal_k" in cloud.attribute_names True

Source code in src/geoai3d/features.py

def multiscale_features(
    cloud: PointCloud,
    *,
    ks: tuple[int, ...] = (10, 20, 30, 40, 50),
) -> PointCloud:
    """Add geometric features computed at each point's optimal neighbourhood.

    Features are computed at several neighbourhood sizes and, per point, the
    size that minimises eigenentropy is chosen (the dimensionality-based scale
    selection of Demantke et al.). The chosen size is stored as ``optimal_k``.

    Args:
        cloud: The cloud to describe.
        ks: Candidate neighbour counts to try. Each must be at least 3.

    Returns:
        A new :class:`PointCloud` with the geometric feature columns at the
        per-point optimal scale, plus an ``optimal_k`` column.

    Raises:
        ValueError: If any candidate in ``ks`` is below 3, or ``ks`` is empty.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, multiscale_features
        >>> pts = np.random.default_rng(0).random((400, 3))
        >>> cloud = multiscale_features(PointCloud(pts, crs=28992), ks=(8, 16, 24))
        >>> "optimal_k" in cloud.attribute_names
        True
    """
    candidates = tuple(sorted({int(value) for value in ks}))
    if not candidates:
        msg = "ks must contain at least one neighbour count."
        raise ValueError(msg)
    if candidates[0] < 3:
        msg = "every candidate in ks must be at least 3."
        raise ValueError(msg)

    xyz = cloud.xyz
    n_points = len(cloud)
    tree = cKDTree(xyz)
    per_scale: list[dict[str, NDArray[Any]]] = []
    entropies = np.empty((n_points, len(candidates)))
    for column, k in enumerate(candidates):
        _, neighbor_indices = tree.query(xyz, k=min(k, n_points), workers=-1)
        scale_features = features_from_eigen(*knn_eigen(xyz[neighbor_indices]))
        per_scale.append(scale_features)
        entropies[:, column] = scale_features["eigenentropy"]

    filled = np.where(np.isnan(entropies), np.inf, entropies)
    best = np.argmin(filled, axis=1)

    chosen: dict[str, NDArray[Any]] = {
        name: np.empty(n_points) for name in FEATURE_NAMES
    }
    for column in range(len(candidates)):
        mask = best == column
        for name in FEATURE_NAMES:
            chosen[name][mask] = per_scale[column][name][mask]
    chosen["optimal_k"] = np.array([candidates[b] for b in best], dtype=np.float64)

    return attach_attributes(
        cloud, chosen, "multiscale_features", {"ks": list(candidates)}
    )

normalize_intensity

normalize_intensity(cloud: PointCloud, *, sensor_position: NDArray[Any] | Sequence[float] | Sequence[Sequence[float]], intensity_attribute: str = 'intensity', reference_range: float | None = None, correct_range: bool = True, correct_incidence_angle: bool = False, max_incidence_angle: float = 80.0, output_attribute: str = 'normalized_intensity') -> PointCloud

Normalise LiDAR intensity for range and, optionally, incidence angle.

Corrects intensity so it better reflects target reflectance and is comparable across a survey, using the sensor geometry (see the module docstring for the model). The result is added as a new float64 column and the CRS and provenance are carried through.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to normalise. Must carry the intensity attribute and, for incidence-angle correction, surface normals (nx, ny, nz from :func:~geoai3d.estimate_normals).	required
sensor_position	NDArray[Any] | Sequence[float] | Sequence[Sequence[float]]	The scanner position, in the cloud's CRS. A (3,) array for a single origin (terrestrial), or an (N, 3) array of per-point positions (airborne or mobile trajectory).	required
intensity_attribute	str	Name of the raw intensity attribute to read.	'intensity'
reference_range	float | None	Range R_ref the range correction normalises to. If omitted, the median sensor-to-point range of the cloud is used.	None
correct_range	bool	If true, apply the (R / R_ref) ** 2 range correction.	True
correct_incidence_angle	bool	If true, also divide by cos(theta) using the surface normals. Off by default because it needs normals.	False
max_incidence_angle	float	Grazing-angle clamp in degrees for the incidence correction; incidence angles beyond it are treated as this angle so the 1 / cos(theta) factor stays bounded.	80.0
output_attribute	str	Name of the normalised intensity column to add.	'normalized_intensity'

Returns:

Type	Description
PointCloud	A new :class:~geoai3d.PointCloud with the normalised intensity column
PointCloud	added.

Raises:

Type	Description
ValueError	If neither correction is enabled; if the intensity or (for incidence correction) the normal attributes are missing; if sensor_position is not (3,) or (N, 3); if any point coincides with the sensor (zero range); if reference_range is not positive; or if max_incidence_angle is not in (0, 90).

Example

import numpy as np from geoai3d import PointCloud, normalize_intensity pts = np.zeros((100, 3)) pts[:, 0] = np.linspace(0.0, 10.0, 100) cloud = PointCloud( ... pts, attributes={"intensity": np.full(100, 100.0)}, crs=28992 ... ) out = normalize_intensity(cloud, sensor_position=[5.0, 0.0, 50.0]) "normalized_intensity" in out.attribute_names True

Source code in src/geoai3d/intensity.py

def normalize_intensity(
    cloud: PointCloud,
    *,
    sensor_position: NDArray[Any] | Sequence[float] | Sequence[Sequence[float]],
    intensity_attribute: str = "intensity",
    reference_range: float | None = None,
    correct_range: bool = True,
    correct_incidence_angle: bool = False,
    max_incidence_angle: float = 80.0,
    output_attribute: str = "normalized_intensity",
) -> PointCloud:
    """Normalise LiDAR intensity for range and, optionally, incidence angle.

    Corrects intensity so it better reflects target reflectance and is
    comparable across a survey, using the sensor geometry (see the module
    docstring for the model). The result is added as a new ``float64`` column
    and the CRS and provenance are carried through.

    Args:
        cloud: The cloud to normalise. Must carry the intensity attribute and,
            for incidence-angle correction, surface normals (``nx``, ``ny``,
            ``nz`` from :func:`~geoai3d.estimate_normals`).
        sensor_position: The scanner position, in the cloud's CRS. A ``(3,)``
            array for a single origin (terrestrial), or an ``(N, 3)`` array of
            per-point positions (airborne or mobile trajectory).
        intensity_attribute: Name of the raw intensity attribute to read.
        reference_range: Range ``R_ref`` the range correction normalises to. If
            omitted, the median sensor-to-point range of the cloud is used.
        correct_range: If true, apply the ``(R / R_ref) ** 2`` range correction.
        correct_incidence_angle: If true, also divide by ``cos(theta)`` using
            the surface normals. Off by default because it needs normals.
        max_incidence_angle: Grazing-angle clamp in degrees for the incidence
            correction; incidence angles beyond it are treated as this angle so
            the ``1 / cos(theta)`` factor stays bounded.
        output_attribute: Name of the normalised intensity column to add.

    Returns:
        A new :class:`~geoai3d.PointCloud` with the normalised intensity column
        added.

    Raises:
        ValueError: If neither correction is enabled; if the intensity or (for
            incidence correction) the normal attributes are missing;
            if ``sensor_position`` is not ``(3,)`` or ``(N, 3)``; if any point
            coincides with the sensor (zero range); if ``reference_range`` is
            not positive; or if ``max_incidence_angle`` is not in ``(0, 90)``.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, normalize_intensity
        >>> pts = np.zeros((100, 3))
        >>> pts[:, 0] = np.linspace(0.0, 10.0, 100)
        >>> cloud = PointCloud(
        ...     pts, attributes={"intensity": np.full(100, 100.0)}, crs=28992
        ... )
        >>> out = normalize_intensity(cloud, sensor_position=[5.0, 0.0, 50.0])
        >>> "normalized_intensity" in out.attribute_names
        True
    """
    if not correct_range and not correct_incidence_angle:
        msg = (
            "Nothing to normalise: enable correct_range and/or correct_incidence_angle."
        )
        raise ValueError(msg)
    if intensity_attribute not in cloud.attribute_names:
        msg = (
            f"Cloud has no {intensity_attribute!r} attribute to normalise. "
            f"Available attributes: {cloud.attribute_names}. Pass "
            "intensity_attribute= if it is stored under another name."
        )
        raise ValueError(msg)

    n_points = len(cloud)
    sensor = np.asarray(sensor_position, dtype=np.float64)
    if sensor.shape == (3,):
        sensor = np.broadcast_to(sensor, (n_points, 3))
    elif sensor.shape != (n_points, 3):
        msg = (
            "sensor_position must have shape (3,) for a single scanner origin "
            f"or ({n_points}, 3) for a per-point trajectory; got {sensor.shape}."
        )
        raise ValueError(msg)

    beam = cloud.xyz - sensor
    ranges = np.linalg.norm(beam, axis=1)
    if np.any(ranges <= 0.0):
        msg = (
            "Some points coincide with the sensor position (zero range), so "
            "intensity cannot be normalised there. Check sensor_position."
        )
        raise ValueError(msg)

    if reference_range is None:
        resolved_reference = float(np.median(ranges))
    elif reference_range <= 0.0:
        msg = "reference_range must be a positive distance."
        raise ValueError(msg)
    else:
        resolved_reference = float(reference_range)

    corrected = cloud.attribute(intensity_attribute).astype(np.float64)
    if correct_range:
        corrected = corrected * (ranges / resolved_reference) ** 2

    if correct_incidence_angle:
        if not 0.0 < max_incidence_angle < 90.0:
            msg = "max_incidence_angle must be between 0 and 90 degrees."
            raise ValueError(msg)
        missing = [n for n in _NORMAL_ATTRIBUTES if n not in cloud.attribute_names]
        if missing:
            msg = (
                "Incidence-angle correction needs surface normals "
                f"{list(_NORMAL_ATTRIBUTES)}; missing {missing}. Run "
                "estimate_normals on the cloud first."
            )
            raise ValueError(msg)
        normals = np.column_stack(
            [cloud.attribute(name).astype(np.float64) for name in _NORMAL_ATTRIBUTES]
        )
        normal_lengths = np.linalg.norm(normals, axis=1)
        normal_lengths = np.where(normal_lengths == 0.0, 1.0, normal_lengths)
        unit_normals = normals / normal_lengths[:, None]
        unit_beam = beam / ranges[:, None]
        # abs(): a surface normal's sign is arbitrary, so fold the beam onto the
        # front face rather than letting an upward normal give a negative cosine.
        cos_incidence = np.abs(np.einsum("ij,ij->i", unit_beam, unit_normals))
        cos_floor = float(np.cos(np.radians(max_incidence_angle)))
        cos_incidence = np.maximum(cos_incidence, cos_floor)
        corrected = corrected / cos_incidence

    return attach_attributes(
        cloud,
        {output_attribute: corrected},
        "normalize_intensity",
        {
            "intensity_attribute": intensity_attribute,
            "reference_range": resolved_reference,
            "correct_range": correct_range,
            "correct_incidence_angle": correct_incidence_angle,
            "max_incidence_angle": (
                max_incidence_angle if correct_incidence_angle else None
            ),
            "sensor_position": (
                "per_point" if np.asarray(sensor_position).ndim == 2 else "origin"
            ),
        },
    )

read_geotiff

read_geotiff(source: str | PathLike[str], *, band: int = 1) -> Raster

Read one band of a GeoTIFF into a :class:~geoai3d.Raster.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to a .tif file.	required
band	int	1-based band index to read.	1

Returns:

Name	Type	Description
A	Raster	class:~geoai3d.Raster with the band's values, transform, CRS, and
	Raster	nodata value.

Raises:

Type	Description
ValueError	If rasterio (the [gis] extra) is not installed, or the file carries no CRS.

Example

import numpy as np import os, tempfile from geoai3d import Raster, to_geotiff, read_geotiff path = os.path.join(tempfile.mkdtemp(), "out.tif") to_geotiff( ... Raster(np.ones((4, 4)), (1.0, 0.0, 0.0, 0.0, -1.0, 4.0), 28992), path ... ) read_geotiff(path).shape (4, 4)

Source code in src/geoai3d/io/geotiff.py

def read_geotiff(
    source: str | os.PathLike[str],
    *,
    band: int = 1,
) -> Raster:
    """Read one band of a GeoTIFF into a :class:`~geoai3d.Raster`.

    Args:
        source: Path to a ``.tif`` file.
        band: 1-based band index to read.

    Returns:
        A :class:`~geoai3d.Raster` with the band's values, transform, CRS, and
        nodata value.

    Raises:
        ValueError: If rasterio (the ``[gis]`` extra) is not installed, or the
            file carries no CRS.

    Example:
        >>> import numpy as np
        >>> import os, tempfile
        >>> from geoai3d import Raster, to_geotiff, read_geotiff
        >>> path = os.path.join(tempfile.mkdtemp(), "out.tif")
        >>> to_geotiff(
        ...     Raster(np.ones((4, 4)), (1.0, 0.0, 0.0, 0.0, -1.0, 4.0), 28992), path
        ... )
        >>> read_geotiff(path).shape
        (4, 4)
    """
    try:
        import rasterio
    except ImportError as exc:
        raise ValueError(_GIS_HINT) from exc

    with rasterio.open(str(source)) as dataset:
        data = dataset.read(band)
        transform = tuple(dataset.transform)[:6]
        file_crs = dataset.crs
        nodata = dataset.nodata
    if file_crs is None:
        msg = (
            f"{str(source)!r} has no coordinate reference system, so it cannot "
            "be read as a georeferenced raster."
        )
        raise ValueError(msg)
    resolved_nodata = float("nan") if nodata is None else float(nodata)
    return Raster(data, transform, file_crs.to_wkt(), nodata=resolved_nodata)

read_las

read_las(source: str | PathLike[str], *, crs: object | None = None) -> PointCloud

Read a LAS or LAZ file into a :class:PointCloud.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to a .las or .laz file.	required
crs	object | None	Coordinate reference system to assign, overriding whatever the file declares. If the file has no CRS and none is given here, an error is raised rather than returning an unreferenced cloud.	None

Returns:

Name	Type	Description
A	PointCloud	class:PointCloud with scaled float64 coordinates, every
	PointCloud	non-coordinate dimension as an attribute, the resolved CRS, and a
	PointCloud	provenance record noting the read.

Raises:

Type	Description
ValueError	If the file has no CRS and crs is not supplied, or if a supplied crs cannot be interpreted.

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, read_las, to_las cloud = PointCloud(np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]), crs=28992) path = os.path.join(tempfile.mkdtemp(), "example.las") to_las(cloud, path) read_las(path).crs.to_epsg() 28992

Source code in src/geoai3d/io/las.py

def read_las(
    source: str | os.PathLike[str],
    *,
    crs: object | None = None,
) -> PointCloud:
    """Read a LAS or LAZ file into a :class:`PointCloud`.

    Args:
        source: Path to a ``.las`` or ``.laz`` file.
        crs: Coordinate reference system to assign, overriding whatever the
            file declares. If the file has no CRS and none is given here, an
            error is raised rather than returning an unreferenced cloud.

    Returns:
        A :class:`PointCloud` with scaled ``float64`` coordinates, every
        non-coordinate dimension as an attribute, the resolved CRS, and a
        provenance record noting the read.

    Raises:
        ValueError: If the file has no CRS and ``crs`` is not supplied, or if a
            supplied ``crs`` cannot be interpreted.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, read_las, to_las
        >>> cloud = PointCloud(np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]), crs=28992)
        >>> path = os.path.join(tempfile.mkdtemp(), "example.las")
        >>> to_las(cloud, path)
        >>> read_las(path).crs.to_epsg()
        28992
    """
    las = laspy.read(str(source))
    header = las.header
    coordinates: NDArray[np.float64] = np.column_stack(
        (np.asarray(las.x), np.asarray(las.y), np.asarray(las.z))
    ).astype(np.float64, copy=False)

    file_crs = header.parse_crs()
    if crs is not None:
        resolved_crs: object = _coerce_crs(crs)
    elif file_crs is not None:
        resolved_crs = file_crs
    else:
        msg = (
            f"{str(source)!r} has no coordinate reference system and none was "
            "supplied. Pass crs= (an EPSG code, WKT string, or pyproj.CRS) to "
            "set it explicitly."
        )
        raise ValueError(msg)

    attributes: dict[str, NDArray[Any]] = {}
    for name in list(las.point_format.dimension_names):
        if name in _COORDINATE_DIMENSIONS:
            continue
        attributes[name] = np.asarray(las[name]).copy()

    provenance = Provenance(source=str(source))
    provenance.add_step(
        "read LAS/LAZ",
        {
            "path": str(source),
            "point_format": int(las.point_format.id),
            "version": str(header.version),
            "point_count": int(header.point_count),
        },
    )
    return PointCloud(
        coordinates,
        attributes=attributes,
        crs=resolved_crs,
        provenance=provenance,
    )

read_lidar

read_lidar(source: str | PathLike[str], *, crs: object | None = None) -> PointCloud

Read a LiDAR or point-cloud file, choosing the reader by extension.

Dispatches on the file suffix: .las and .laz are read with :func:~geoai3d.read_las, and whitespace-delimited text (.xyz, .txt, .asc, .pts) with :func:~geoai3d.read_xyz using its defaults. For delimited text, non-default column layouts, or to force a particular format, call the format-specific reader directly.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to the point-cloud file.	required
crs	object | None	Coordinate reference system to assign, overriding any the file carries. An EPSG code, WKT/PROJ string, or pyproj.CRS. Required when the file (and, for text, its sibling .prj) has none.	None

Returns:

Name	Type	Description
A	PointCloud	class:~geoai3d.PointCloud with the resolved CRS and a provenance
	PointCloud	record noting the read.

Raises:

Type	Description
ValueError	If the extension is .parquet (read it with :func:~geoai3d.read_parquet), or is not a recognised point-cloud format.

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, read_lidar, to_las cloud = PointCloud(np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]), crs=28992) path = os.path.join(tempfile.mkdtemp(), "scan.las") to_las(cloud, path) read_lidar(path).crs.to_epsg() 28992

Source code in src/geoai3d/io/dispatch.py

def read_lidar(
    source: str | os.PathLike[str],
    *,
    crs: object | None = None,
) -> PointCloud:
    """Read a LiDAR or point-cloud file, choosing the reader by extension.

    Dispatches on the file suffix: ``.las`` and ``.laz`` are read with
    :func:`~geoai3d.read_las`, and whitespace-delimited text
    (``.xyz``, ``.txt``, ``.asc``, ``.pts``) with :func:`~geoai3d.read_xyz`
    using its defaults. For delimited text, non-default column layouts, or to
    force a particular format, call the format-specific reader directly.

    Args:
        source: Path to the point-cloud file.
        crs: Coordinate reference system to assign, overriding any the file
            carries. An EPSG code, WKT/PROJ string, or ``pyproj.CRS``. Required
            when the file (and, for text, its sibling ``.prj``) has none.

    Returns:
        A :class:`~geoai3d.PointCloud` with the resolved CRS and a provenance
        record noting the read.

    Raises:
        ValueError: If the extension is ``.parquet`` (read it with
            :func:`~geoai3d.read_parquet`), or is not a recognised
            point-cloud format.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, read_lidar, to_las
        >>> cloud = PointCloud(np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]), crs=28992)
        >>> path = os.path.join(tempfile.mkdtemp(), "scan.las")
        >>> to_las(cloud, path)
        >>> read_lidar(path).crs.to_epsg()
        28992
    """
    suffix = Path(source).suffix.lower()
    if suffix in _LAS_SUFFIXES:
        return read_las(source, crs=crs)
    if suffix in _TEXT_SUFFIXES:
        return read_xyz(source, crs=crs)
    if suffix == ".parquet":
        msg = (
            "Parquet is GEOAI_3D's serialization format, not an acquisition "
            "format; read it with read_parquet() rather than read_lidar()."
        )
        raise ValueError(msg)
    supported = ", ".join(_LAS_SUFFIXES + _TEXT_SUFFIXES)
    msg = (
        f"Unsupported point-cloud format {suffix!r} for {str(source)!r}. "
        f"read_lidar handles: {supported}. For Parquet use read_parquet()."
    )
    raise ValueError(msg)

read_parquet

read_parquet(source: str | PathLike[str], *, crs: object | None = None) -> PointCloud

Read a Parquet file written by :func:to_parquet into a cloud.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to a .parquet file with x, y, z columns.	required
crs	object | None	Coordinate reference system to assign, overriding the file's metadata. If the file carries no CRS and none is given here, an error is raised rather than returning an unreferenced cloud.	None

Returns:

Name	Type	Description
A	PointCloud	class:PointCloud with the coordinates, every other column as an
	PointCloud	attribute, the resolved CRS, and the stored provenance record.

Raises:

Type	Description
ValueError	If the file lacks x/y/z columns, or has no CRS and crs is not supplied.

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, to_parquet, read_parquet path = os.path.join(tempfile.mkdtemp(), "cloud.parquet") to_parquet(PointCloud(np.ones((2, 3)), crs=28992), path) read_parquet(path).crs.to_epsg() 28992

Source code in src/geoai3d/io/parquet.py

def read_parquet(
    source: str | os.PathLike[str],
    *,
    crs: object | None = None,
) -> PointCloud:
    """Read a Parquet file written by :func:`to_parquet` into a cloud.

    Args:
        source: Path to a ``.parquet`` file with ``x``, ``y``, ``z`` columns.
        crs: Coordinate reference system to assign, overriding the file's
            metadata. If the file carries no CRS and none is given here, an
            error is raised rather than returning an unreferenced cloud.

    Returns:
        A :class:`PointCloud` with the coordinates, every other column as an
        attribute, the resolved CRS, and the stored provenance record.

    Raises:
        ValueError: If the file lacks ``x``/``y``/``z`` columns, or has no CRS
            and ``crs`` is not supplied.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_parquet, read_parquet
        >>> path = os.path.join(tempfile.mkdtemp(), "cloud.parquet")
        >>> to_parquet(PointCloud(np.ones((2, 3)), crs=28992), path)
        >>> read_parquet(path).crs.to_epsg()
        28992
    """
    table = pq.read_table(str(source))
    missing = [name for name in _COORDINATE_COLUMNS if name not in table.column_names]
    if missing:
        msg = (
            f"{str(source)!r} is missing coordinate columns {missing}; it does "
            "not look like a geoai3d Parquet file."
        )
        raise ValueError(msg)

    coordinates = np.column_stack(
        [np.asarray(table.column(name)) for name in _COORDINATE_COLUMNS]
    ).astype(np.float64, copy=False)
    attributes: dict[str, NDArray[Any]] = {
        name: np.asarray(table.column(name))
        for name in table.column_names
        if name not in _COORDINATE_COLUMNS
    }

    metadata = table.schema.metadata or {}
    if crs is not None:
        resolved_crs: object = coerce_crs(crs)
    elif _CRS_KEY in metadata:
        resolved_crs = pyproj.CRS.from_wkt(metadata[_CRS_KEY].decode("utf-8"))
    else:
        msg = (
            f"{str(source)!r} has no coordinate reference system in its "
            "metadata and none was supplied. Pass crs= to set it explicitly."
        )
        raise ValueError(msg)

    provenance = None
    if _PROVENANCE_KEY in metadata:
        provenance = Provenance.from_dict(
            json.loads(metadata[_PROVENANCE_KEY].decode("utf-8"))
        )

    return PointCloud(
        coordinates,
        attributes=attributes,
        crs=resolved_crs,
        provenance=provenance,
    )

read_xyz

read_xyz(source: str | PathLike[str], *, crs: object | None = None, columns: tuple[str, ...] = _COORDINATE_NAMES, delimiter: str | None = None, comments: str = '#') -> PointCloud

Read a whitespace- or delimiter-separated XYZ text file.

Parameters:

Name	Type	Description	Default
source	str | PathLike[str]	Path to the text file.	required
crs	object | None	Coordinate reference system to assign. If omitted, a sibling .prj file is used; if neither is available, an error is raised.	None
columns	tuple[str, ...]	Name of each column in order. Must include x, y, and z; any other names become attributes.	_COORDINATE_NAMES
delimiter	str | None	Column delimiter. None (default) splits on any run of whitespace.	None
comments	str	Character marking comment lines to skip.	'#'

Returns:

Name	Type	Description
A	PointCloud	class:PointCloud with the resolved CRS and a provenance record.

Raises:

Type	Description
ValueError	If columns omits a coordinate axis, the column count does not match, or no CRS can be resolved.

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, read_xyz, to_xyz path = os.path.join(tempfile.mkdtemp(), "points.xyz") to_xyz(PointCloud(np.zeros((3, 3)), crs=28992), path) len(read_xyz(path)) 3

Source code in src/geoai3d/io/xyz.py

def read_xyz(
    source: str | os.PathLike[str],
    *,
    crs: object | None = None,
    columns: tuple[str, ...] = _COORDINATE_NAMES,
    delimiter: str | None = None,
    comments: str = "#",
) -> PointCloud:
    """Read a whitespace- or delimiter-separated XYZ text file.

    Args:
        source: Path to the text file.
        crs: Coordinate reference system to assign. If omitted, a sibling
            ``.prj`` file is used; if neither is available, an error is raised.
        columns: Name of each column in order. Must include ``x``, ``y``, and
            ``z``; any other names become attributes.
        delimiter: Column delimiter. ``None`` (default) splits on any run of
            whitespace.
        comments: Character marking comment lines to skip.

    Returns:
        A :class:`PointCloud` with the resolved CRS and a provenance record.

    Raises:
        ValueError: If ``columns`` omits a coordinate axis, the column count
            does not match, or no CRS can be resolved.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, read_xyz, to_xyz
        >>> path = os.path.join(tempfile.mkdtemp(), "points.xyz")
        >>> to_xyz(PointCloud(np.zeros((3, 3)), crs=28992), path)
        >>> len(read_xyz(path))
        3
    """
    names = tuple(columns)
    for axis in _COORDINATE_NAMES:
        if axis not in names:
            msg = f"columns must include {axis!r}; got {names}."
            raise ValueError(msg)
    data: NDArray[np.float64] = np.loadtxt(
        str(source), delimiter=delimiter, comments=comments, dtype=np.float64, ndmin=2
    )
    if data.shape[1] != len(names):
        msg = (
            f"{str(source)!r} has {data.shape[1]} columns but {len(names)} "
            f"column names were given: {names}."
        )
        raise ValueError(msg)

    resolved_crs = _resolve_crs(source, crs)
    by_name = {name: data[:, index] for index, name in enumerate(names)}
    coordinates = np.column_stack([by_name["x"], by_name["y"], by_name["z"]])
    attributes: dict[str, NDArray[Any]] = {
        name: by_name[name] for name in names if name not in _COORDINATE_NAMES
    }

    provenance = Provenance(source=str(source))
    provenance.add_step("read XYZ", {"path": str(source), "columns": list(names)})
    return PointCloud(
        coordinates,
        attributes=attributes,
        crs=resolved_crs,
        provenance=provenance,
    )

region_growing

region_growing(cloud: PointCloud, *, k: int = 30, smoothness_degrees: float = 15.0, curvature_attribute: str = 'surface_variation', curvature_threshold: float | None = None, min_size: int = 10, output_attribute: str = 'segment') -> PointCloud

Segment a cloud into smooth regions by growing from seeds.

Starting from the smoothest points, a region absorbs a neighbour when the angle between their surface normals is below smoothness_degrees; the neighbour then continues the growth unless a curvature limit stops it. This is the classical normal-based region growing (after Rabbani et al.), which separates a scene by surface orientation -- ground from walls from roofs -- with no training data.

The cloud must carry surface normals (nx, ny, nz from :func:~geoai3d.estimate_normals). If it also carries a curvature attribute (for example surface_variation from :func:~geoai3d.geometric_features) the smoothest points seed first, which gives cleaner regions.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to segment. Must carry surface normals.	required
k	int	Number of nearest neighbours considered for growth.	30
smoothness_degrees	float	Maximum normal-to-normal angle, in degrees, for a neighbour to join a region.	15.0
curvature_attribute	str	Attribute used to order seeds (smoothest first) and, with curvature_threshold, to limit propagation. Ignored if the cloud does not carry it.	'surface_variation'
curvature_threshold	float | None	If set, only points with curvature below it continue growing a region. None lets every joined point continue.	None
min_size	int	Minimum points for a region to be kept; smaller regions are labelled -1.	10
output_attribute	str	Name of the integer segment label column to add.	'segment'

Returns:

Type	Description
PointCloud	A new :class:~geoai3d.PointCloud with the integer segment label column
PointCloud	added (contiguous ids, -1 for unassigned points), the CRS and
PointCloud	provenance carried through.

Raises:

Type	Description
ValueError	If the cloud is empty, lacks the normal attributes, k or min_size is below 1, or smoothness_degrees is not in (0, 180).

Example

import numpy as np from geoai3d import PointCloud, estimate_normals, region_growing rng = np.random.default_rng(0) flat = rng.uniform(0.0, 5.0, (400, 3)) flat[:, 2] = 0.0 cloud = estimate_normals(PointCloud(flat, crs=28992), k=12) segmented = region_growing(cloud, k=12, min_size=20) int(segmented.attribute("segment").max()) 0

Source code in src/geoai3d/cluster.py

def region_growing(
    cloud: PointCloud,
    *,
    k: int = 30,
    smoothness_degrees: float = 15.0,
    curvature_attribute: str = "surface_variation",
    curvature_threshold: float | None = None,
    min_size: int = 10,
    output_attribute: str = "segment",
) -> PointCloud:
    """Segment a cloud into smooth regions by growing from seeds.

    Starting from the smoothest points, a region absorbs a neighbour when the
    angle between their surface normals is below ``smoothness_degrees``; the
    neighbour then continues the growth unless a curvature limit stops it. This
    is the classical normal-based region growing (after Rabbani et al.), which
    separates a scene by surface orientation -- ground from walls from roofs --
    with no training data.

    The cloud must carry surface normals (``nx``, ``ny``, ``nz`` from
    :func:`~geoai3d.estimate_normals`). If it also carries a curvature attribute
    (for example ``surface_variation`` from :func:`~geoai3d.geometric_features`)
    the smoothest points seed first, which gives cleaner regions.

    Args:
        cloud: The cloud to segment. Must carry surface normals.
        k: Number of nearest neighbours considered for growth.
        smoothness_degrees: Maximum normal-to-normal angle, in degrees, for a
            neighbour to join a region.
        curvature_attribute: Attribute used to order seeds (smoothest first) and,
            with ``curvature_threshold``, to limit propagation. Ignored if the
            cloud does not carry it.
        curvature_threshold: If set, only points with curvature below it continue
            growing a region. ``None`` lets every joined point continue.
        min_size: Minimum points for a region to be kept; smaller regions are
            labelled ``-1``.
        output_attribute: Name of the integer segment label column to add.

    Returns:
        A new :class:`~geoai3d.PointCloud` with the integer segment label column
        added (contiguous ids, ``-1`` for unassigned points), the CRS and
        provenance carried through.

    Raises:
        ValueError: If the cloud is empty, lacks the normal attributes, ``k`` or
            ``min_size`` is below 1, or ``smoothness_degrees`` is not in
            ``(0, 180)``.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, estimate_normals, region_growing
        >>> rng = np.random.default_rng(0)
        >>> flat = rng.uniform(0.0, 5.0, (400, 3))
        >>> flat[:, 2] = 0.0
        >>> cloud = estimate_normals(PointCloud(flat, crs=28992), k=12)
        >>> segmented = region_growing(cloud, k=12, min_size=20)
        >>> int(segmented.attribute("segment").max())
        0
    """
    n_points = len(cloud)
    if n_points == 0:
        msg = "Cannot segment an empty cloud."
        raise ValueError(msg)
    if k < 1:
        msg = "k must be at least 1."
        raise ValueError(msg)
    if min_size < 1:
        msg = "min_size must be at least 1."
        raise ValueError(msg)
    if not 0.0 < smoothness_degrees < 180.0:
        msg = "smoothness_degrees must be between 0 and 180."
        raise ValueError(msg)
    missing = [n for n in _NORMAL_ATTRIBUTES if n not in cloud.attribute_names]
    if missing:
        msg = (
            "region_growing needs surface normals "
            f"{list(_NORMAL_ATTRIBUTES)}; missing {missing}. Run estimate_normals "
            "on the cloud first."
        )
        raise ValueError(msg)

    normals = np.column_stack(
        [cloud.attribute(name).astype(np.float64) for name in _NORMAL_ATTRIBUTES]
    )
    lengths = np.linalg.norm(normals, axis=1)
    lengths = np.where(lengths == 0.0, 1.0, lengths)
    normals = normals / lengths[:, None]

    tree = cKDTree(cloud.xyz)
    neighbor_count = min(k + 1, n_points)
    _, neighbors = tree.query(cloud.xyz, k=neighbor_count, workers=-1)
    neighbors = np.atleast_2d(neighbors)[:, 1:]  # drop each point's self match

    has_curvature = curvature_attribute in cloud.attribute_names
    if has_curvature:
        curvature = cloud.attribute(curvature_attribute).astype(np.float64)
        seed_order = np.argsort(curvature)
    else:
        curvature = None
        seed_order = np.arange(n_points)

    cos_threshold = float(np.cos(np.radians(smoothness_degrees)))
    labels = np.full(n_points, -1, dtype=np.int64)
    visited = np.zeros(n_points, dtype=bool)
    next_label = 0
    for seed in seed_order:
        if visited[seed]:
            continue
        queue = [int(seed)]
        visited[seed] = True
        region = [int(seed)]
        while queue:
            point = queue.pop()
            for neighbor in neighbors[point]:
                neighbor = int(neighbor)
                if visited[neighbor]:
                    continue
                if abs(float(normals[point] @ normals[neighbor])) < cos_threshold:
                    continue
                visited[neighbor] = True
                region.append(neighbor)
                propagate = curvature_threshold is None or (
                    curvature is not None and curvature[neighbor] < curvature_threshold
                )
                if propagate:
                    queue.append(neighbor)
        if len(region) >= min_size:
            labels[region] = next_label
            next_label += 1

    return attach_attributes(
        cloud,
        {output_attribute: labels},
        "region_growing",
        {
            "k": k,
            "smoothness_degrees": smoothness_degrees,
            "curvature_threshold": curvature_threshold,
            "min_size": min_size,
        },
    )

remove_statistical_outliers

remove_statistical_outliers(cloud: PointCloud, *, k: int = 16, std_ratio: float = 2.0) -> PointCloud

Remove points whose neighbours are unusually far away.

For each point, the mean distance to its k nearest neighbours is computed. Points whose mean distance exceeds the global mean by more than std_ratio standard deviations are dropped.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to filter.	required
k	int	Number of nearest neighbours per point.	16
std_ratio	float	Multiplier of the global standard deviation above which a point is treated as an outlier. Smaller values remove more points.	2.0

Returns:

Type	Description
PointCloud	A new :class:PointCloud with the outliers removed.

Raises:

Type	Description
ValueError	If k is less than 1.

Example

import numpy as np from geoai3d import PointCloud, remove_statistical_outliers pts = np.random.default_rng(0).random((500, 3)) pts[0] = [100.0, 100.0, 100.0] # a clear outlier cleaned = remove_statistical_outliers(PointCloud(pts, crs=28992)) len(cleaned) < 500 True

Source code in src/geoai3d/geometry.py

def remove_statistical_outliers(
    cloud: PointCloud,
    *,
    k: int = 16,
    std_ratio: float = 2.0,
) -> PointCloud:
    """Remove points whose neighbours are unusually far away.

    For each point, the mean distance to its ``k`` nearest neighbours is
    computed. Points whose mean distance exceeds the global mean by more than
    ``std_ratio`` standard deviations are dropped.

    Args:
        cloud: The cloud to filter.
        k: Number of nearest neighbours per point.
        std_ratio: Multiplier of the global standard deviation above which a
            point is treated as an outlier. Smaller values remove more points.

    Returns:
        A new :class:`PointCloud` with the outliers removed.

    Raises:
        ValueError: If ``k`` is less than 1.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, remove_statistical_outliers
        >>> pts = np.random.default_rng(0).random((500, 3))
        >>> pts[0] = [100.0, 100.0, 100.0]  # a clear outlier
        >>> cleaned = remove_statistical_outliers(PointCloud(pts, crs=28992))
        >>> len(cleaned) < 500
        True
    """
    if k < 1:
        msg = "k must be at least 1."
        raise ValueError(msg)
    xyz = cloud.xyz
    k_effective = min(k, len(cloud) - 1)
    tree = cKDTree(xyz)
    distances, _ = tree.query(xyz, k=k_effective + 1, workers=-1)
    # Column 0 is the point itself (distance 0); average the rest.
    mean_distance = distances[:, 1:].mean(axis=1)
    threshold = mean_distance.mean() + std_ratio * mean_distance.std()
    keep = np.flatnonzero(mean_distance <= threshold)
    return select(
        cloud, keep, "remove_statistical_outliers", {"k": k, "std_ratio": std_ratio}
    )

reproject

reproject(cloud: PointCloud, target_crs: object) -> PointCloud

Reproject a cloud's horizontal coordinates to another CRS.

Only the x and y coordinates are transformed; heights are carried through unchanged. For a transform that also converts heights between vertical datums, use :func:reproject_3d.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to reproject. Its CRS must be set.	required
target_crs	object	Target CRS as an EPSG code, WKT string, or pyproj.CRS.	required

Returns:

Type	Description
PointCloud	A new :class:PointCloud in target_crs, with attributes carried
PointCloud	through and a provenance step appended.

Raises:

Type	Description
ValueError	If the cloud has no CRS, or a CRS cannot be interpreted.

Example

import numpy as np from geoai3d import PointCloud, reproject cloud = PointCloud(np.array([[155000.0, 463000.0, 0.0]]), crs=28992) reproject(cloud, 4326).crs.to_epsg() 4326

Source code in src/geoai3d/crs.py

def reproject(cloud: PointCloud, target_crs: object) -> PointCloud:
    """Reproject a cloud's horizontal coordinates to another CRS.

    Only the ``x`` and ``y`` coordinates are transformed; heights are carried
    through unchanged. For a transform that also converts heights between
    vertical datums, use :func:`reproject_3d`.

    Args:
        cloud: The cloud to reproject. Its CRS must be set.
        target_crs: Target CRS as an EPSG code, WKT string, or ``pyproj.CRS``.

    Returns:
        A new :class:`PointCloud` in ``target_crs``, with attributes carried
        through and a provenance step appended.

    Raises:
        ValueError: If the cloud has no CRS, or a CRS cannot be interpreted.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, reproject
        >>> cloud = PointCloud(np.array([[155000.0, 463000.0, 0.0]]), crs=28992)
        >>> reproject(cloud, 4326).crs.to_epsg()
        4326
    """
    if cloud.crs is None:
        msg = "Cannot reproject a cloud with no CRS. Set the cloud's CRS first."
        raise ValueError(msg)
    source = coerce_crs(cloud.crs)
    target = coerce_crs(target_crs)
    transformer = pyproj.Transformer.from_crs(source, target, always_xy=True)

    x, y = transformer.transform(cloud.xyz[:, 0], cloud.xyz[:, 1])
    coordinates = np.column_stack([np.asarray(x), np.asarray(y), cloud.xyz[:, 2]])

    provenance = _provenance_for(cloud)
    provenance.add_step(
        "reproject (horizontal)",
        {"source_crs": source.to_string(), "target_crs": target.to_string()},
    )
    return PointCloud(
        coordinates,
        attributes=_attributes_of(cloud),
        crs=target,
        provenance=provenance,
    )

reproject_3d

reproject_3d(cloud: PointCloud, target_crs: object, *, allow_network: bool = False) -> PointCloud

Reproject a cloud in three dimensions, including the vertical datum.

Transforms x, y, and z together, so heights are converted between the source and target vertical datums (for example ellipsoidal to NAP orthometric) using the appropriate geoid grid. Both the source and target CRS must carry a vertical axis.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to reproject. Its CRS must be set and three-dimensional.	required
target_crs	object	A 3D or compound target CRS (for example 4979 for WGS84 ellipsoidal, or 7415 for RD New + NAP orthometric).	required
allow_network	bool	If true, let pyproj download any missing transform grid from the PROJ CDN. If false and a grid is missing, an error is raised rather than returning silently wrong heights.	False

Returns:

Type	Description
PointCloud	A new :class:PointCloud in target_crs, with attributes carried
PointCloud	through and a provenance step appended.

Raises:

Type	Description
ValueError	If the cloud has no CRS, either CRS lacks a vertical axis, or the transform produced non-finite coordinates (usually a missing geoid grid).

Example

import numpy as np from geoai3d import PointCloud, reproject_3d cloud = PointCloud(np.array([[5.0, 52.0, 43.0]]), crs=4979) reproject_3d(cloud, 4937).crs.to_epsg() 4937

Source code in src/geoai3d/crs.py

def reproject_3d(
    cloud: PointCloud,
    target_crs: object,
    *,
    allow_network: bool = False,
) -> PointCloud:
    """Reproject a cloud in three dimensions, including the vertical datum.

    Transforms ``x``, ``y``, and ``z`` together, so heights are converted
    between the source and target vertical datums (for example ellipsoidal to
    NAP orthometric) using the appropriate geoid grid. Both the source and
    target CRS must carry a vertical axis.

    Args:
        cloud: The cloud to reproject. Its CRS must be set and three-dimensional.
        target_crs: A 3D or compound target CRS (for example ``4979`` for
            WGS84 ellipsoidal, or ``7415`` for RD New + NAP orthometric).
        allow_network: If true, let pyproj download any missing transform grid
            from the PROJ CDN. If false and a grid is missing, an error is
            raised rather than returning silently wrong heights.

    Returns:
        A new :class:`PointCloud` in ``target_crs``, with attributes carried
        through and a provenance step appended.

    Raises:
        ValueError: If the cloud has no CRS, either CRS lacks a vertical axis,
            or the transform produced non-finite coordinates (usually a missing
            geoid grid).

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, reproject_3d
        >>> cloud = PointCloud(np.array([[5.0, 52.0, 43.0]]), crs=4979)
        >>> reproject_3d(cloud, 4937).crs.to_epsg()
        4937
    """
    if cloud.crs is None:
        msg = "Cannot reproject a cloud with no CRS. Set the cloud's CRS first."
        raise ValueError(msg)
    source = coerce_crs(cloud.crs)
    target = coerce_crs(target_crs)
    for role, crs_obj in (("source", source), ("target", target)):
        if not _has_vertical_axis(crs_obj):
            msg = (
                f"The {role} CRS {crs_obj.to_string()!r} has no vertical axis, "
                "so a 3D transform cannot convert heights. Use reproject() for "
                "a horizontal transform, or supply a 3D/compound CRS such as "
                "EPSG:4979 (ellipsoidal) or EPSG:7415 (RD New + NAP)."
            )
            raise ValueError(msg)

    if allow_network:
        # pyproj's type stubs do not export this runtime function.
        pyproj.network.set_network_enabled(active=True)  # type: ignore[attr-defined]
    transformer = pyproj.Transformer.from_crs(source, target, always_xy=True)

    x, y, z = transformer.transform(cloud.xyz[:, 0], cloud.xyz[:, 1], cloud.xyz[:, 2])
    coordinates = np.column_stack([np.asarray(x), np.asarray(y), np.asarray(z)])
    if not bool(np.all(np.isfinite(coordinates))):
        msg = (
            "The 3D transform produced non-finite coordinates, usually because "
            "the geoid/datum grid it needs is not available offline. Retry with "
            "allow_network=True, or install the PROJ grid locally."
        )
        raise ValueError(msg)

    provenance = _provenance_for(cloud)
    provenance.add_step(
        "reproject (3D)",
        {"source_crs": source.to_string(), "target_crs": target.to_string()},
    )
    return PointCloud(
        coordinates,
        attributes=_attributes_of(cloud),
        crs=target,
        provenance=provenance,
    )

save_classifier

save_classifier(model: Classifier, destination: str | PathLike[str]) -> None

Persist a trained classifier to disk (pickle).

Parameters:

Name	Type	Description	Default
model	Classifier	The classifier to save.	required
destination	str | PathLike[str]	Output path.	required

Source code in src/geoai3d/classify.py

def save_classifier(model: Classifier, destination: str | os.PathLike[str]) -> None:
    """Persist a trained classifier to disk (pickle).

    Args:
        model: The classifier to save.
        destination: Output path.
    """
    with open(destination, "wb") as handle:
        pickle.dump(model, handle)

spatial_block_split

spatial_block_split(cloud: PointCloud, *, block_size: float, n_folds: int = 5, seed: int = 0) -> list[tuple[NDArray[Any], NDArray[Any]]]

Split a cloud into spatially-blocked cross-validation folds.

Points are grouped into square block_size blocks in the horizontal plane, the blocks are shuffled and dealt round-robin into n_folds folds, and each fold's test set is the points of its blocks with the rest as the training set. Because a whole block goes to one fold, test points never neighbour their own training points.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to split.	required
block_size	float	Edge length of the square blocks, in the cloud's CRS units. Make it several times the feature neighbourhood so blocks are genuinely independent.	required
n_folds	int	Number of folds.	5
seed	int	Seed for shuffling blocks into folds.	0

Returns:

Type	Description
list[tuple[NDArray[Any], NDArray[Any]]]	A list of n_folds (train_indices, test_indices) pairs of integer
list[tuple[NDArray[Any], NDArray[Any]]]	index arrays into the cloud.

Raises:

Type	Description
ValueError	If the cloud is empty, block_size is not positive, n_folds is below 2, or there are fewer blocks than folds.

Example

import numpy as np from geoai3d import PointCloud, spatial_block_split pts = np.random.default_rng(0).uniform(0.0, 100.0, (2000, 3)) folds = spatial_block_split( ... PointCloud(pts, crs=28992), block_size=10.0, n_folds=5 ... ) len(folds) 5

Source code in src/geoai3d/crossval.py

def spatial_block_split(
    cloud: PointCloud,
    *,
    block_size: float,
    n_folds: int = 5,
    seed: int = 0,
) -> list[tuple[NDArray[Any], NDArray[Any]]]:
    """Split a cloud into spatially-blocked cross-validation folds.

    Points are grouped into square ``block_size`` blocks in the horizontal
    plane, the blocks are shuffled and dealt round-robin into ``n_folds`` folds,
    and each fold's test set is the points of its blocks with the rest as the
    training set. Because a whole block goes to one fold, test points never
    neighbour their own training points.

    Args:
        cloud: The cloud to split.
        block_size: Edge length of the square blocks, in the cloud's CRS units.
            Make it several times the feature neighbourhood so blocks are
            genuinely independent.
        n_folds: Number of folds.
        seed: Seed for shuffling blocks into folds.

    Returns:
        A list of ``n_folds`` ``(train_indices, test_indices)`` pairs of integer
        index arrays into the cloud.

    Raises:
        ValueError: If the cloud is empty, ``block_size`` is not positive,
            ``n_folds`` is below 2, or there are fewer blocks than folds.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, spatial_block_split
        >>> pts = np.random.default_rng(0).uniform(0.0, 100.0, (2000, 3))
        >>> folds = spatial_block_split(
        ...     PointCloud(pts, crs=28992), block_size=10.0, n_folds=5
        ... )
        >>> len(folds)
        5
    """
    n_points = len(cloud)
    if n_points == 0:
        msg = "Cannot split an empty cloud."
        raise ValueError(msg)
    if block_size <= 0.0:
        msg = "block_size must be a positive number."
        raise ValueError(msg)
    if n_folds < 2:
        msg = "n_folds must be at least 2."
        raise ValueError(msg)

    xyz = cloud.xyz
    block_col = np.floor((xyz[:, 0] - xyz[:, 0].min()) / block_size).astype(np.int64)
    block_row = np.floor((xyz[:, 1] - xyz[:, 1].min()) / block_size).astype(np.int64)
    # A single integer key per block; np.unique gives each point its block index.
    _, block_index = np.unique(
        np.column_stack([block_col, block_row]), axis=0, return_inverse=True
    )
    block_index = np.asarray(block_index).ravel()
    n_blocks = int(block_index.max()) + 1
    if n_blocks < n_folds:
        msg = (
            f"Only {n_blocks} spatial blocks for {n_folds} folds. Use a smaller "
            "block_size or fewer folds."
        )
        raise ValueError(msg)

    rng = np.random.default_rng(seed)
    fold_of_block = rng.permutation(n_blocks) % n_folds
    point_fold = fold_of_block[block_index]

    folds: list[tuple[NDArray[Any], NDArray[Any]]] = []
    for fold in range(n_folds):
        test_indices = np.flatnonzero(point_fold == fold)
        train_indices = np.flatnonzero(point_fold != fold)
        folds.append((train_indices, test_indices))
    return folds

to_dsm

to_dsm(cloud: PointCloud, *, resolution: float, bounds: tuple[float, float, float, float] | None = None, nodata: float = float('nan')) -> Raster

Rasterise a cloud to a Digital Surface Model (the top surface).

Uses the highest z per cell over every point, capturing the top of vegetation and buildings.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to rasterise. Its CRS must be set.	required
resolution	float	Cell size in the cloud's CRS units.	required
bounds	tuple[float, float, float, float] | None	Grid extent (min_x, min_y, max_x, max_y). Defaults to the cloud's full horizontal extent.	None
nodata	float	Value for empty cells. Defaults to NaN.	float('nan')

Returns:

Name	Type	Description
A	Raster	class:~geoai3d.Raster of surface heights. Write it to disk with
	Raster	func:~geoai3d.to_geotiff.

Raises:

Type	Description
ValueError	If the cloud is empty or has no CRS, or resolution is not positive.

Example

import numpy as np from geoai3d import PointCloud, to_dsm pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3)) dsm = to_dsm(PointCloud(pts, crs=28992), resolution=1.0) dsm.shape (10, 10)

Source code in src/geoai3d/dem.py

def to_dsm(
    cloud: PointCloud,
    *,
    resolution: float,
    bounds: tuple[float, float, float, float] | None = None,
    nodata: float = float("nan"),
) -> Raster:
    """Rasterise a cloud to a Digital Surface Model (the top surface).

    Uses the highest z per cell over every point, capturing the top of
    vegetation and buildings.

    Args:
        cloud: The cloud to rasterise. Its CRS must be set.
        resolution: Cell size in the cloud's CRS units.
        bounds: Grid extent ``(min_x, min_y, max_x, max_y)``. Defaults to the
            cloud's full horizontal extent.
        nodata: Value for empty cells. Defaults to ``NaN``.

    Returns:
        A :class:`~geoai3d.Raster` of surface heights. Write it to disk with
        :func:`~geoai3d.to_geotiff`.

    Raises:
        ValueError: If the cloud is empty or has no CRS, or ``resolution`` is
            not positive.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_dsm
        >>> pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3))
        >>> dsm = to_dsm(PointCloud(pts, crs=28992), resolution=1.0)
        >>> dsm.shape
        (10, 10)
    """
    _check_inputs(cloud, resolution)
    grid_bounds = bounds if bounds is not None else _cloud_bounds(cloud)
    return _rasterize(cloud.xyz, grid_bounds, cloud.crs, resolution, "max", nodata)

to_dtm

to_dtm(cloud: PointCloud, *, resolution: float, ground_attribute: str = 'is_ground', bounds: tuple[float, float, float, float] | None = None, nodata: float = float('nan')) -> Raster

Rasterise a cloud to a Digital Terrain Model (bare earth).

If the cloud carries a boolean ground attribute (from :func:~geoai3d.ground), only those points are used; otherwise every point is used and the lowest z per cell approximates the bare earth.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to rasterise. Its CRS must be set.	required
resolution	float	Cell size in the cloud's CRS units.	required
ground_attribute	str	Boolean attribute selecting ground points, if present.	'is_ground'
bounds	tuple[float, float, float, float] | None	Grid extent (min_x, min_y, max_x, max_y). Defaults to the cloud's full horizontal extent, so a DTM and DSM from the same cloud align for differencing.	None
nodata	float	Value for empty cells. Defaults to NaN.	float('nan')

Returns:

Name	Type	Description
A	Raster	class:~geoai3d.Raster of ground heights. Write it to disk with
	Raster	func:~geoai3d.to_geotiff.

Raises:

Type	Description
ValueError	If the cloud is empty or has no CRS, resolution is not positive, or the ground attribute is present but all false.

Example

import numpy as np from geoai3d import PointCloud, to_dtm pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3)) dtm = to_dtm(PointCloud(pts, crs=28992), resolution=1.0) dtm.shape (10, 10)

Source code in src/geoai3d/dem.py

def to_dtm(
    cloud: PointCloud,
    *,
    resolution: float,
    ground_attribute: str = "is_ground",
    bounds: tuple[float, float, float, float] | None = None,
    nodata: float = float("nan"),
) -> Raster:
    """Rasterise a cloud to a Digital Terrain Model (bare earth).

    If the cloud carries a boolean ground attribute (from
    :func:`~geoai3d.ground`), only those points are used; otherwise every point
    is used and the lowest z per cell approximates the bare earth.

    Args:
        cloud: The cloud to rasterise. Its CRS must be set.
        resolution: Cell size in the cloud's CRS units.
        ground_attribute: Boolean attribute selecting ground points, if present.
        bounds: Grid extent ``(min_x, min_y, max_x, max_y)``. Defaults to the
            cloud's full horizontal extent, so a DTM and DSM from the same cloud
            align for differencing.
        nodata: Value for empty cells. Defaults to ``NaN``.

    Returns:
        A :class:`~geoai3d.Raster` of ground heights. Write it to disk with
        :func:`~geoai3d.to_geotiff`.

    Raises:
        ValueError: If the cloud is empty or has no CRS, ``resolution`` is not
            positive, or the ground attribute is present but all false.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_dtm
        >>> pts = np.random.default_rng(0).uniform(0.0, 10.0, (2000, 3))
        >>> dtm = to_dtm(PointCloud(pts, crs=28992), resolution=1.0)
        >>> dtm.shape
        (10, 10)
    """
    _check_inputs(cloud, resolution)
    if ground_attribute in cloud.attribute_names:
        mask = np.asarray(cloud.attribute(ground_attribute)).astype(bool)
        xyz = cloud.xyz[mask]
        if xyz.shape[0] == 0:
            msg = (
                f"No ground points to rasterise: attribute {ground_attribute!r} "
                "is all false."
            )
            raise ValueError(msg)
    else:
        xyz = cloud.xyz
    grid_bounds = bounds if bounds is not None else _cloud_bounds(cloud)
    return _rasterize(xyz, grid_bounds, cloud.crs, resolution, "min", nodata)

to_geopackage

to_geopackage(cloud: PointCloud, destination: str | PathLike[str], *, attributes: Sequence[str] | None = None, layer: str = 'points', include_z: bool = True) -> None

Write a point cloud to a GeoPackage vector layer.

Each point becomes a (3D) point feature; the chosen attributes become columns. The coordinate reference system is written into the file. If the cloud carries a provenance record, it is written to a sibling <destination>.provenance.json.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to write. Its CRS must be set.	required
destination	str | PathLike[str]	Output .gpkg path.	required
attributes	Sequence[str] | None	Attribute columns to include. If omitted, every attribute is written.	None
layer	str	Layer name inside the GeoPackage.	'points'
include_z	bool	If true, write 3D PointZ geometry carrying elevation; if false, write 2D points.	True

Raises:

Type	Description
ValueError	If geopandas (the [gis] extra) is not installed, or the cloud has no CRS.

Example

import numpy as np import os, tempfile from geoai3d import PointCloud, to_geopackage cloud = PointCloud( ... np.zeros((3, 3)), ... attributes={"classification": np.array([2, 2, 6])}, ... crs=28992, ... ) to_geopackage(cloud, os.path.join(tempfile.mkdtemp(), "points.gpkg"))

Source code in src/geoai3d/io/vector.py

def to_geopackage(
    cloud: PointCloud,
    destination: str | os.PathLike[str],
    *,
    attributes: Sequence[str] | None = None,
    layer: str = "points",
    include_z: bool = True,
) -> None:
    """Write a point cloud to a GeoPackage vector layer.

    Each point becomes a (3D) point feature; the chosen attributes become
    columns. The coordinate reference system is written into the file. If the
    cloud carries a provenance record, it is written to a sibling
    ``<destination>.provenance.json``.

    Args:
        cloud: The cloud to write. Its CRS must be set.
        destination: Output ``.gpkg`` path.
        attributes: Attribute columns to include. If omitted, every attribute is
            written.
        layer: Layer name inside the GeoPackage.
        include_z: If true, write 3D ``PointZ`` geometry carrying elevation; if
            false, write 2D points.

    Raises:
        ValueError: If geopandas (the ``[gis]`` extra) is not installed, or the
            cloud has no CRS.

    Example:
        >>> import numpy as np
        >>> import os, tempfile
        >>> from geoai3d import PointCloud, to_geopackage
        >>> cloud = PointCloud(
        ...     np.zeros((3, 3)),
        ...     attributes={"classification": np.array([2, 2, 6])},
        ...     crs=28992,
        ... )
        >>> to_geopackage(cloud, os.path.join(tempfile.mkdtemp(), "points.gpkg"))
    """
    try:
        import geopandas
        import shapely
    except ImportError as exc:
        raise ValueError(_GIS_HINT) from exc

    if cloud.crs is None:
        msg = (
            "Cannot write a GeoPackage without a coordinate reference system. "
            "Set the cloud's CRS first."
        )
        raise ValueError(msg)

    names = list(cloud.attribute_names) if attributes is None else list(attributes)
    missing = [name for name in names if name not in cloud.attribute_names]
    if missing:
        msg = (
            f"Attributes {missing} are not present. Available: {cloud.attribute_names}."
        )
        raise ValueError(msg)

    coordinates = cloud.xyz
    if include_z:
        geometry = shapely.points(
            coordinates[:, 0], coordinates[:, 1], coordinates[:, 2]
        )
    else:
        geometry = shapely.points(coordinates[:, 0], coordinates[:, 1])

    columns = {name: cloud.attribute(name) for name in names}
    frame = geopandas.GeoDataFrame(columns, geometry=geometry, crs=cloud.crs)
    frame.to_file(str(destination), driver="GPKG", layer=layer, engine="pyogrio")

    if cloud.provenance is not None:
        sidecar = Path(f"{destination}.provenance.json")
        sidecar.write_text(
            json.dumps(cloud.provenance.to_dict(), indent=2), encoding="utf-8"
        )

to_geotiff

to_geotiff(raster: Raster, destination: str | PathLike[str]) -> None

Write a :class:~geoai3d.Raster to a single-band GeoTIFF.

Parameters:

Name	Type	Description	Default
raster	Raster	The raster to write. Its CRS, transform, and nodata value are stored in the file.	required
destination	str | PathLike[str]	Output .tif path.	required

Raises:

Type	Description
ValueError	If rasterio (the [gis] extra) is not installed.

Example

import numpy as np import os, tempfile from geoai3d import Raster, to_geotiff raster = Raster(np.ones((4, 4)), (1.0, 0.0, 0.0, 0.0, -1.0, 4.0), 28992) to_geotiff(raster, os.path.join(tempfile.mkdtemp(), "out.tif"))

Source code in src/geoai3d/io/geotiff.py

def to_geotiff(raster: Raster, destination: str | os.PathLike[str]) -> None:
    """Write a :class:`~geoai3d.Raster` to a single-band GeoTIFF.

    Args:
        raster: The raster to write. Its CRS, transform, and nodata value are
            stored in the file.
        destination: Output ``.tif`` path.

    Raises:
        ValueError: If rasterio (the ``[gis]`` extra) is not installed.

    Example:
        >>> import numpy as np
        >>> import os, tempfile
        >>> from geoai3d import Raster, to_geotiff
        >>> raster = Raster(np.ones((4, 4)), (1.0, 0.0, 0.0, 0.0, -1.0, 4.0), 28992)
        >>> to_geotiff(raster, os.path.join(tempfile.mkdtemp(), "out.tif"))
    """
    try:
        import rasterio
    except ImportError as exc:
        raise ValueError(_GIS_HINT) from exc

    data = np.asarray(raster.data)
    with rasterio.open(
        str(destination),
        "w",
        driver="GTiff",
        height=raster.height,
        width=raster.width,
        count=1,
        dtype=str(data.dtype),
        crs=raster.crs.to_wkt(),
        transform=rasterio.Affine(*raster.transform),
        nodata=raster.nodata,
    ) as dataset:
        dataset.write(data, 1)

to_las

to_las(cloud: PointCloud, destination: str | PathLike[str], *, point_format: int = _DEFAULT_POINT_FORMAT, file_version: str = _DEFAULT_FILE_VERSION, scales: NDArray[float64] | Sequence[float] | None = None, offsets: NDArray[float64] | Sequence[float] | None = None) -> None

Write a :class:PointCloud to a LAS or LAZ file.

Attributes whose names are standard dimensions of point_format are written as those fields; any other attribute is written as a user-defined extra dimension, preserving its values and dtype.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to write. Its CRS must be set.	required
destination	str | PathLike[str]	Output path. A .laz suffix triggers compression and requires the lazrs backend (pip install geoai3d[laz]).	required
point_format	int	LAS point format id to write. The default, 6, carries coordinates, intensity, returns, classification, scan angle, GPS time, and point source id.	_DEFAULT_POINT_FORMAT
file_version	str	LAS version string. The default, "1.4", stores the CRS as WKT, which supports arbitrary coordinate systems.	_DEFAULT_FILE_VERSION
scales	NDArray[float64] | Sequence[float] | None	Per-axis coordinate scales. Defaults to one millimetre on each axis.	None
offsets	NDArray[float64] | Sequence[float] | None	Per-axis coordinate offsets. Defaults to the minimum coordinate on each axis.	None

Raises:

Type	Description
ValueError	If the cloud has no CRS, or the CRS cannot be interpreted.

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, to_las cloud = PointCloud(np.zeros((5, 3)), crs=28992) to_las(cloud, os.path.join(tempfile.mkdtemp(), "out.laz"))

Source code in src/geoai3d/io/las.py

def to_las(
    cloud: PointCloud,
    destination: str | os.PathLike[str],
    *,
    point_format: int = _DEFAULT_POINT_FORMAT,
    file_version: str = _DEFAULT_FILE_VERSION,
    scales: NDArray[np.float64] | Sequence[float] | None = None,
    offsets: NDArray[np.float64] | Sequence[float] | None = None,
) -> None:
    """Write a :class:`PointCloud` to a LAS or LAZ file.

    Attributes whose names are standard dimensions of ``point_format`` are
    written as those fields; any other attribute is written as a user-defined
    extra dimension, preserving its values and dtype.

    Args:
        cloud: The cloud to write. Its CRS must be set.
        destination: Output path. A ``.laz`` suffix triggers compression and
            requires the ``lazrs`` backend (``pip install geoai3d[laz]``).
        point_format: LAS point format id to write. The default, 6, carries
            coordinates, intensity, returns, classification, scan angle,
            GPS time, and point source id.
        file_version: LAS version string. The default, ``"1.4"``, stores the
            CRS as WKT, which supports arbitrary coordinate systems.
        scales: Per-axis coordinate scales. Defaults to one millimetre on each
            axis.
        offsets: Per-axis coordinate offsets. Defaults to the minimum
            coordinate on each axis.

    Raises:
        ValueError: If the cloud has no CRS, or the CRS cannot be interpreted.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_las
        >>> cloud = PointCloud(np.zeros((5, 3)), crs=28992)
        >>> to_las(cloud, os.path.join(tempfile.mkdtemp(), "out.laz"))
    """
    if cloud.crs is None:
        msg = (
            "Cannot write a LAS/LAZ file without a coordinate reference "
            "system. Set the cloud's CRS first, for example by reading it "
            "with read_las(..., crs=...)."
        )
        raise ValueError(msg)
    resolved_crs = _coerce_crs(cloud.crs)

    coordinates = cloud.xyz
    if offsets is None:
        offsets_array = (
            np.min(coordinates, axis=0) if len(cloud) else np.zeros(3, dtype=np.float64)
        )
    else:
        offsets_array = np.asarray(offsets, dtype=np.float64)
    scales_array = (
        np.full(3, _DEFAULT_SCALE)
        if scales is None
        else np.asarray(scales, dtype=np.float64)
    )

    header = laspy.LasHeader(version=file_version, point_format=point_format)
    header.offsets = offsets_array
    header.scales = scales_array

    standard_names = set(header.point_format.standard_dimension_names)
    for name in cloud.attribute_names:
        if name in standard_names or name in _COORDINATE_DIMENSIONS:
            continue
        values = cloud.attribute(name)
        header.add_extra_dim(laspy.ExtraBytesParams(name=name, type=values.dtype))

    header.add_crs(resolved_crs)

    las = laspy.LasData(header)
    # Assigning the coordinates first sizes the point record; other
    # dimensions are then assigned into the allocated record.
    las.x = coordinates[:, 0]
    las.y = coordinates[:, 1]
    las.z = coordinates[:, 2]
    for name in cloud.attribute_names:
        if name in _COORDINATE_DIMENSIONS:
            continue
        las[name] = cloud.attribute(name)

    las.write(str(destination))

to_parquet

to_parquet(cloud: PointCloud, destination: str | PathLike[str]) -> None

Write a :class:PointCloud to a Parquet file.

Coordinates are stored as x, y, z columns and each attribute as a column of its own dtype. The CRS and provenance record are written into the schema metadata.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to write. Its CRS must be set.	required
destination	str | PathLike[str]	Output .parquet path.	required

Raises:

Type	Description
ValueError	If the cloud has no CRS, or an attribute is named x, y, or z (reserved for coordinates).

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, to_parquet, read_parquet cloud = PointCloud(np.zeros((4, 3)), crs=28992) path = os.path.join(tempfile.mkdtemp(), "cloud.parquet") to_parquet(cloud, path) len(read_parquet(path)) 4

Source code in src/geoai3d/io/parquet.py

def to_parquet(cloud: PointCloud, destination: str | os.PathLike[str]) -> None:
    """Write a :class:`PointCloud` to a Parquet file.

    Coordinates are stored as ``x``, ``y``, ``z`` columns and each attribute as
    a column of its own dtype. The CRS and provenance record are written into
    the schema metadata.

    Args:
        cloud: The cloud to write. Its CRS must be set.
        destination: Output ``.parquet`` path.

    Raises:
        ValueError: If the cloud has no CRS, or an attribute is named ``x``,
            ``y``, or ``z`` (reserved for coordinates).

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_parquet, read_parquet
        >>> cloud = PointCloud(np.zeros((4, 3)), crs=28992)
        >>> path = os.path.join(tempfile.mkdtemp(), "cloud.parquet")
        >>> to_parquet(cloud, path)
        >>> len(read_parquet(path))
        4
    """
    if cloud.crs is None:
        msg = (
            "Cannot write a Parquet file without a coordinate reference "
            "system. Set the cloud's CRS first."
        )
        raise ValueError(msg)
    reserved = set(_COORDINATE_COLUMNS) & set(cloud.attribute_names)
    if reserved:
        msg = (
            f"Attributes {sorted(reserved)} clash with the reserved coordinate "
            "column names 'x', 'y', 'z'. Rename them before writing."
        )
        raise ValueError(msg)

    coordinates = cloud.xyz
    columns: dict[str, NDArray[Any]] = {
        "x": np.ascontiguousarray(coordinates[:, 0]),
        "y": np.ascontiguousarray(coordinates[:, 1]),
        "z": np.ascontiguousarray(coordinates[:, 2]),
    }
    for name in cloud.attribute_names:
        columns[name] = cloud.attribute(name)

    metadata: dict[bytes, bytes] = {
        _SCHEMA_VERSION_KEY: _SCHEMA_VERSION,
        _CRS_KEY: coerce_crs(cloud.crs).to_wkt().encode("utf-8"),
    }
    if cloud.provenance is not None:
        metadata[_PROVENANCE_KEY] = json.dumps(cloud.provenance.to_dict()).encode(
            "utf-8"
        )

    table = pa.table(columns).replace_schema_metadata(metadata)
    pq.write_table(table, str(destination))

to_xyz

to_xyz(cloud: PointCloud, destination: str | PathLike[str], *, include_attributes: bool = True, delimiter: str = ' ', fmt: str = '%.6f', write_header: bool = False, write_prj: bool = True) -> None

Write a :class:PointCloud to an XYZ text file.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to write.	required
destination	str | PathLike[str]	Output path.	required
include_attributes	bool	If true, write attribute columns after x y z.	True
delimiter	str	Column delimiter.	' '
fmt	str	numpy.savetxt format applied to every column.	'%.6f'
write_header	bool	If true, write a commented header line of column names.	False
write_prj	bool	If true and the cloud has a CRS, write it as WKT to a sibling .prj file so read_xyz can recover it.	True

Example

import os, tempfile import numpy as np from geoai3d import PointCloud, to_xyz cloud = PointCloud(np.zeros((2, 3)), crs=28992) to_xyz(cloud, os.path.join(tempfile.mkdtemp(), "out.xyz"))

Source code in src/geoai3d/io/xyz.py

def to_xyz(
    cloud: PointCloud,
    destination: str | os.PathLike[str],
    *,
    include_attributes: bool = True,
    delimiter: str = " ",
    fmt: str = "%.6f",
    write_header: bool = False,
    write_prj: bool = True,
) -> None:
    """Write a :class:`PointCloud` to an XYZ text file.

    Args:
        cloud: The cloud to write.
        destination: Output path.
        include_attributes: If true, write attribute columns after ``x y z``.
        delimiter: Column delimiter.
        fmt: ``numpy.savetxt`` format applied to every column.
        write_header: If true, write a commented header line of column names.
        write_prj: If true and the cloud has a CRS, write it as WKT to a
            sibling ``.prj`` file so ``read_xyz`` can recover it.

    Example:
        >>> import os, tempfile
        >>> import numpy as np
        >>> from geoai3d import PointCloud, to_xyz
        >>> cloud = PointCloud(np.zeros((2, 3)), crs=28992)
        >>> to_xyz(cloud, os.path.join(tempfile.mkdtemp(), "out.xyz"))
    """
    coordinates = cloud.xyz
    arrays: list[NDArray[Any]] = [
        coordinates[:, 0],
        coordinates[:, 1],
        coordinates[:, 2],
    ]
    names = ["x", "y", "z"]
    if include_attributes:
        for name in cloud.attribute_names:
            arrays.append(np.asarray(cloud.attribute(name), dtype=np.float64))
            names.append(name)

    table = np.column_stack(arrays)
    header = delimiter.join(names) if write_header else ""
    np.savetxt(
        str(destination),
        table,
        delimiter=delimiter,
        fmt=fmt,
        header=header,
        comments="# ",
    )
    if write_prj and cloud.crs is not None:
        _prj_path(destination).write_text(
            coerce_crs(cloud.crs).to_wkt(), encoding="utf-8"
        )

train_classifier

train_classifier(cloud: PointCloud, *, label_attribute: str = 'classification', feature_names: tuple[str, ...] | list[str] | None = None, n_estimators: int = 200, max_depth: int | None = None, random_state: int = 0) -> Classifier

Train a random-forest classifier on a labelled cloud.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The labelled training cloud.	required
label_attribute	str	Attribute holding the integer class label per point.	'classification'
feature_names	tuple[str, ...] | list[str] | None	Feature attributes to train on. If omitted, resolved as in :func:feature_matrix.	None
n_estimators	int	Number of trees in the forest.	200
max_depth	int | None	Maximum tree depth, or None for unlimited.	None
random_state	int	Seed for reproducible training.	0

Returns:

Name	Type	Description
A	Classifier	class:Classifier wrapping the fitted forest and its feature names.

Raises:

Type	Description
ValueError	If the label attribute is absent, or features cannot be resolved.

Example

import numpy as np from geoai3d import PointCloud, train_classifier, classify rng = np.random.default_rng(0) planarity = np.concatenate([rng.normal(0.9, 0.05, 100), ... rng.normal(0.1, 0.05, 100)]) labels = np.concatenate([np.zeros(100), np.ones(100)]).astype(int) cloud = PointCloud( ... np.zeros((200, 3)), ... attributes={"planarity": planarity, "classification": labels}, ... crs=28992, ... ) model = train_classifier(cloud, feature_names=["planarity"]) model.feature_names ('planarity',)

Source code in src/geoai3d/classify.py

def train_classifier(
    cloud: PointCloud,
    *,
    label_attribute: str = "classification",
    feature_names: tuple[str, ...] | list[str] | None = None,
    n_estimators: int = 200,
    max_depth: int | None = None,
    random_state: int = 0,
) -> Classifier:
    """Train a random-forest classifier on a labelled cloud.

    Args:
        cloud: The labelled training cloud.
        label_attribute: Attribute holding the integer class label per point.
        feature_names: Feature attributes to train on. If omitted, resolved as
            in :func:`feature_matrix`.
        n_estimators: Number of trees in the forest.
        max_depth: Maximum tree depth, or ``None`` for unlimited.
        random_state: Seed for reproducible training.

    Returns:
        A :class:`Classifier` wrapping the fitted forest and its feature names.

    Raises:
        ValueError: If the label attribute is absent, or features cannot be
            resolved.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud, train_classifier, classify
        >>> rng = np.random.default_rng(0)
        >>> planarity = np.concatenate([rng.normal(0.9, 0.05, 100),
        ...                             rng.normal(0.1, 0.05, 100)])
        >>> labels = np.concatenate([np.zeros(100), np.ones(100)]).astype(int)
        >>> cloud = PointCloud(
        ...     np.zeros((200, 3)),
        ...     attributes={"planarity": planarity, "classification": labels},
        ...     crs=28992,
        ... )
        >>> model = train_classifier(cloud, feature_names=["planarity"])
        >>> model.feature_names
        ('planarity',)
    """
    if label_attribute not in cloud.attribute_names:
        msg = (
            f"Cloud has no {label_attribute!r} label attribute to train on. "
            "Pass label_attribute= to name the training labels."
        )
        raise ValueError(msg)
    names = _resolve_features(cloud, feature_names)
    features = feature_matrix(cloud, names)
    labels = np.asarray(cloud.attribute(label_attribute))

    from sklearn.ensemble import RandomForestClassifier

    estimator = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        random_state=random_state,
        n_jobs=-1,
    )
    estimator.fit(features, labels)
    return Classifier(estimator=estimator, feature_names=tuple(names))

view

view(cloud: PointCloud, *, color_by: str | None = None, max_points: int = _DEFAULT_MAX_POINTS, point_size: float = 2.0, seed: int = 0) -> go.Figure

Render a point cloud as an interactive 3D plot.

Parameters:

Name	Type	Description	Default
cloud	PointCloud	The cloud to display.	required
color_by	str | None	Attribute name to colour by, or "z"/"height", or "rgb" for true colour. Defaults to RGB when red/green/blue attributes are present, otherwise height.	None
max_points	int	Randomly thin the cloud to at most this many points for display, so the browser stays responsive.	_DEFAULT_MAX_POINTS
point_size	float	Marker size.	2.0
seed	int	Seed for the display thinning, for a reproducible view.	0

Returns:

Type	Description
Figure	A plotly.graph_objects.Figure. In a notebook it renders inline;
Figure	elsewhere call .show() on it.

Raises:

Type	Description
ImportError	If plotly is not installed.
ValueError	If color_by names something that cannot be used.

Example

import numpy as np from geoai3d import PointCloud from geoai3d.viz import view cloud = PointCloud(np.random.default_rng(0).random((500, 3)), crs=28992) figure = view(cloud) figure.data[0].type 'scatter3d'

Source code in src/geoai3d/viz.py

def view(
    cloud: PointCloud,
    *,
    color_by: str | None = None,
    max_points: int = _DEFAULT_MAX_POINTS,
    point_size: float = 2.0,
    seed: int = 0,
) -> go.Figure:
    """Render a point cloud as an interactive 3D plot.

    Args:
        cloud: The cloud to display.
        color_by: Attribute name to colour by, or ``"z"``/``"height"``, or
            ``"rgb"`` for true colour. Defaults to RGB when red/green/blue
            attributes are present, otherwise height.
        max_points: Randomly thin the cloud to at most this many points for
            display, so the browser stays responsive.
        point_size: Marker size.
        seed: Seed for the display thinning, for a reproducible view.

    Returns:
        A ``plotly.graph_objects.Figure``. In a notebook it renders inline;
        elsewhere call ``.show()`` on it.

    Raises:
        ImportError: If plotly is not installed.
        ValueError: If ``color_by`` names something that cannot be used.

    Example:
        >>> import numpy as np
        >>> from geoai3d import PointCloud
        >>> from geoai3d.viz import view
        >>> cloud = PointCloud(np.random.default_rng(0).random((500, 3)), crs=28992)
        >>> figure = view(cloud)
        >>> figure.data[0].type
        'scatter3d'
    """
    try:
        import plotly.graph_objects as go
    except ImportError as exc:  # pragma: no cover - exercised without the extra
        msg = (
            "view() needs plotly. Install the viewer extra with "
            "'pip install geoai3d[viz]'."
        )
        raise ImportError(msg) from exc

    display = cloud
    if len(cloud) > max_points:
        display = subsample(cloud, method="random", count=max_points, seed=seed)

    xyz = display.xyz
    color, colorscale, colorbar_title = _resolve_color(display, color_by)
    if colorscale is None:
        marker = {"size": point_size, "color": color}
    else:
        marker = {
            "size": point_size,
            "color": color,
            "colorscale": colorscale,
            "colorbar": {"title": colorbar_title},
        }

    figure = go.Figure(
        data=[
            go.Scatter3d(
                x=xyz[:, 0],
                y=xyz[:, 1],
                z=xyz[:, 2],
                mode="markers",
                marker=marker,
            )
        ]
    )
    figure.update_layout(
        scene={"aspectmode": "data"},
        margin={"l": 0, "r": 0, "t": 0, "b": 0},
    )
    return figure

volume

volume(raster: Raster, *, base: float = 0.0) -> dict[str, float]

Compute cut and fill volumes of a raster relative to a base level.

Interprets the raster as a height field (typically a difference raster) and integrates it against base: cells above base contribute fill, cells below contribute cut, each weighted by the cell area. Nodata cells are ignored.

Parameters:

Name	Type	Description	Default
raster	Raster	The height raster to integrate.	required
base	float	The reference level to measure against.	0.0

Returns:

Type	Description
dict[str, float]	A dict with "cut", "fill", and "net" (fill minus cut)
dict[str, float]	volumes, in cubic CRS units.

Example

import numpy as np from geoai3d import Raster, volume data = np.array([[1.0, 1.0], [1.0, 1.0]]) # 1 m over four 2 m cells raster = Raster(data, (2.0, 0.0, 0.0, 0.0, -2.0, 4.0), 28992) volume(raster)["fill"] 16.0

Source code in src/geoai3d/dem.py

def volume(raster: Raster, *, base: float = 0.0) -> dict[str, float]:
    """Compute cut and fill volumes of a raster relative to a base level.

    Interprets the raster as a height field (typically a difference raster) and
    integrates it against ``base``: cells above ``base`` contribute fill, cells
    below contribute cut, each weighted by the cell area. Nodata cells are
    ignored.

    Args:
        raster: The height raster to integrate.
        base: The reference level to measure against.

    Returns:
        A dict with ``"cut"``, ``"fill"``, and ``"net"`` (fill minus cut)
        volumes, in cubic CRS units.

    Example:
        >>> import numpy as np
        >>> from geoai3d import Raster, volume
        >>> data = np.array([[1.0, 1.0], [1.0, 1.0]])  # 1 m over four 2 m cells
        >>> raster = Raster(data, (2.0, 0.0, 0.0, 0.0, -2.0, 4.0), 28992)
        >>> volume(raster)["fill"]
        16.0
    """
    data = _finite(raster)
    valid = ~np.isnan(data)
    difference_from_base = data - base
    res_x, res_y = raster.resolution
    cell_area = res_x * res_y
    fill = float(
        np.sum(difference_from_base[valid & (difference_from_base > 0.0)]) * cell_area
    )
    cut = float(
        -np.sum(difference_from_base[valid & (difference_from_base < 0.0)]) * cell_area
    )
    return {"cut": cut, "fill": fill, "net": fill - cut}