Skip to content

Point Cloud API

obia.pointcloud

obia.pointcloud

Optional point-cloud feature extraction for OBIA segment objects.

add_pointcloud_features(segments, pointcloud, *, metrics=DEFAULT_METRICS, crs=None, reader='auto', predicate='within')

Add point-cloud features to segment objects.

Parameters

segments: Segment polygons. The output preserves the same rows and appends pc_* feature columns. pointcloud: LAS/LAZ path or in-memory point cloud. Path inputs require optional pointcloud dependencies. metrics: Feature groups to calculate: "height", "intensity", and "density". crs: CRS assigned to point clouds that do not already have one. reader: Backend for path inputs: "auto", "pdal", or "pyforestscan". predicate: Spatial join predicate used to attach points to polygons.

Source code in obia/pointcloud/segment_features.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def add_pointcloud_features(
    segments: gpd.GeoDataFrame,
    pointcloud: str | Path | np.ndarray | pd.DataFrame | gpd.GeoDataFrame,
    *,
    metrics: Iterable[str] = DEFAULT_METRICS,
    crs=None,
    reader: str = "auto",
    predicate: str = "within",
) -> gpd.GeoDataFrame:
    """Add point-cloud features to segment objects.

    Parameters
    ----------
    segments:
        Segment polygons. The output preserves the same rows and appends
        ``pc_*`` feature columns.
    pointcloud:
        LAS/LAZ path or in-memory point cloud. Path inputs require optional
        pointcloud dependencies.
    metrics:
        Feature groups to calculate: ``"height"``, ``"intensity"``, and
        ``"density"``.
    crs:
        CRS assigned to point clouds that do not already have one.
    reader:
        Backend for path inputs: ``"auto"``, ``"pdal"``, or ``"pyforestscan"``.
    predicate:
        Spatial join predicate used to attach points to polygons.
    """
    if segments.empty:
        return segments.copy()
    if segments.geometry.name is None:
        raise ValueError("segments must have an active geometry column")

    points = read_pointcloud(pointcloud, crs=crs or segments.crs, reader=reader)
    if points.empty:
        return _append_empty_features(segments, metrics)

    if segments.crs is not None:
        if points.crs is None:
            points = points.set_crs(segments.crs)
        elif points.crs != segments.crs:
            points = points.to_crs(segments.crs)

    segment_index_name = "__obia_segment_index"
    point_index_name = "__obia_point_index"

    segment_table = segments[[segments.geometry.name]].copy()
    segment_table[segment_index_name] = segments.index

    point_table = points.copy()
    point_table[point_index_name] = np.arange(len(point_table))

    joined = gpd.sjoin(
        point_table,
        segment_table,
        how="inner",
        predicate=predicate,
    )

    features_by_segment = {}
    for segment_index, group in joined.groupby(segment_index_name):
        points_in_segment = point_table.iloc[group[point_index_name].to_numpy()]
        area = float(segments.loc[segment_index].geometry.area)
        features_by_segment[segment_index] = calculate_pointcloud_features(
            points_in_segment,
            metrics=metrics,
            area=area,
        )

    empty_features = calculate_pointcloud_features(pd.DataFrame(), metrics=metrics, area=np.nan)
    output = segments.copy()
    for column in empty_features:
        output[column] = np.nan

    for segment_index, features in features_by_segment.items():
        for column, value in features.items():
            output.loc[segment_index, column] = value

    output["pc_point_count"] = output["pc_point_count"].fillna(0).astype(int)
    return output

calculate_pointcloud_features(points, *, metrics=DEFAULT_METRICS, area=None, height_column='z', intensity_column='intensity')

Calculate point-cloud features for one segment.

Parameters

points: Points inside one segment. The table may contain z and intensity columns. metrics: Feature groups to calculate. Supported values are "height", "intensity", and "density". area: Segment area in CRS units squared. Required for density. height_column: Column used for height metrics. intensity_column: Column used for intensity metrics.

Source code in obia/pointcloud/features.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def calculate_pointcloud_features(
    points: pd.DataFrame,
    *,
    metrics: Iterable[str] = DEFAULT_METRICS,
    area: float | None = None,
    height_column: str = "z",
    intensity_column: str = "intensity",
) -> dict[str, float]:
    """Calculate point-cloud features for one segment.

    Parameters
    ----------
    points:
        Points inside one segment. The table may contain ``z`` and
        ``intensity`` columns.
    metrics:
        Feature groups to calculate. Supported values are ``"height"``,
        ``"intensity"``, and ``"density"``.
    area:
        Segment area in CRS units squared. Required for density.
    height_column:
        Column used for height metrics.
    intensity_column:
        Column used for intensity metrics.
    """
    metric_set = set(metrics)
    unknown = metric_set.difference(DEFAULT_METRICS)
    if unknown:
        raise ValueError(f"Unknown point-cloud metrics: {sorted(unknown)}")

    features: dict[str, float] = {"pc_point_count": int(len(points))}

    if "density" in metric_set:
        if area is None or area <= 0:
            features["pc_density"] = np.nan
        else:
            features["pc_density"] = float(len(points) / area)

    if "height" in metric_set:
        features.update(_numeric_series_features(points, height_column, "pc_z", HEIGHT_PERCENTILES))

    if "intensity" in metric_set:
        features.update(
            _numeric_series_features(points, intensity_column, "pc_intensity", INTENSITY_PERCENTILES)
        )

    return features

read_pointcloud(pointcloud, *, crs=None, reader='auto')

Return point-cloud data as a point GeoDataFrame.

Parameters

pointcloud: LAS/LAZ path, structured NumPy array, (n, >=3) NumPy array, pandas DataFrame, or GeoDataFrame. crs: CRS assigned to in-memory point clouds when one is not already present. CRS alignment with segment polygons is checked later by add_pointcloud_features. reader: Reader used for path inputs. "auto" tries PDAL first and then pyforestscan. "pdal" and "pyforestscan" force one backend.

Source code in obia/pointcloud/io.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def read_pointcloud(
    pointcloud: str | Path | np.ndarray | pd.DataFrame | gpd.GeoDataFrame,
    *,
    crs: Any = None,
    reader: str = "auto",
) -> gpd.GeoDataFrame:
    """Return point-cloud data as a point GeoDataFrame.

    Parameters
    ----------
    pointcloud:
        LAS/LAZ path, structured NumPy array, ``(n, >=3)`` NumPy array,
        pandas DataFrame, or GeoDataFrame.
    crs:
        CRS assigned to in-memory point clouds when one is not already present.
        CRS alignment with segment polygons is checked later by
        ``add_pointcloud_features``.
    reader:
        Reader used for path inputs. ``"auto"`` tries PDAL first and then
        pyforestscan. ``"pdal"`` and ``"pyforestscan"`` force one backend.
    """
    if isinstance(pointcloud, gpd.GeoDataFrame):
        gdf = pointcloud.copy()
        if gdf.crs is None and crs is not None:
            gdf = gdf.set_crs(crs)
        return _normalize_point_geodataframe(gdf)

    if isinstance(pointcloud, pd.DataFrame):
        return _dataframe_to_points(pointcloud.copy(), crs=crs)

    if isinstance(pointcloud, np.ndarray):
        return _array_to_points(pointcloud, crs=crs)

    if isinstance(pointcloud, (str, Path)):
        return _read_pointcloud_path(Path(pointcloud), crs=crs, reader=reader)

    raise TypeError(
        "pointcloud must be a path, NumPy array, pandas DataFrame, or GeoDataFrame"
    )

obia.pointcloud.io

obia.pointcloud.io

Point-cloud input normalization.

Path-based LAS/LAZ reading is optional and requires either PDAL or pyforestscan. In-memory arrays and GeoDataFrames are supported without those dependencies, which keeps the core package importable in lightweight environments.

read_pointcloud(pointcloud, *, crs=None, reader='auto')

Return point-cloud data as a point GeoDataFrame.

Parameters

pointcloud: LAS/LAZ path, structured NumPy array, (n, >=3) NumPy array, pandas DataFrame, or GeoDataFrame. crs: CRS assigned to in-memory point clouds when one is not already present. CRS alignment with segment polygons is checked later by add_pointcloud_features. reader: Reader used for path inputs. "auto" tries PDAL first and then pyforestscan. "pdal" and "pyforestscan" force one backend.

Source code in obia/pointcloud/io.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def read_pointcloud(
    pointcloud: str | Path | np.ndarray | pd.DataFrame | gpd.GeoDataFrame,
    *,
    crs: Any = None,
    reader: str = "auto",
) -> gpd.GeoDataFrame:
    """Return point-cloud data as a point GeoDataFrame.

    Parameters
    ----------
    pointcloud:
        LAS/LAZ path, structured NumPy array, ``(n, >=3)`` NumPy array,
        pandas DataFrame, or GeoDataFrame.
    crs:
        CRS assigned to in-memory point clouds when one is not already present.
        CRS alignment with segment polygons is checked later by
        ``add_pointcloud_features``.
    reader:
        Reader used for path inputs. ``"auto"`` tries PDAL first and then
        pyforestscan. ``"pdal"`` and ``"pyforestscan"`` force one backend.
    """
    if isinstance(pointcloud, gpd.GeoDataFrame):
        gdf = pointcloud.copy()
        if gdf.crs is None and crs is not None:
            gdf = gdf.set_crs(crs)
        return _normalize_point_geodataframe(gdf)

    if isinstance(pointcloud, pd.DataFrame):
        return _dataframe_to_points(pointcloud.copy(), crs=crs)

    if isinstance(pointcloud, np.ndarray):
        return _array_to_points(pointcloud, crs=crs)

    if isinstance(pointcloud, (str, Path)):
        return _read_pointcloud_path(Path(pointcloud), crs=crs, reader=reader)

    raise TypeError(
        "pointcloud must be a path, NumPy array, pandas DataFrame, or GeoDataFrame"
    )

obia.pointcloud.features

obia.pointcloud.features

Point-cloud feature calculations.

calculate_pointcloud_features(points, *, metrics=DEFAULT_METRICS, area=None, height_column='z', intensity_column='intensity')

Calculate point-cloud features for one segment.

Parameters

points: Points inside one segment. The table may contain z and intensity columns. metrics: Feature groups to calculate. Supported values are "height", "intensity", and "density". area: Segment area in CRS units squared. Required for density. height_column: Column used for height metrics. intensity_column: Column used for intensity metrics.

Source code in obia/pointcloud/features.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def calculate_pointcloud_features(
    points: pd.DataFrame,
    *,
    metrics: Iterable[str] = DEFAULT_METRICS,
    area: float | None = None,
    height_column: str = "z",
    intensity_column: str = "intensity",
) -> dict[str, float]:
    """Calculate point-cloud features for one segment.

    Parameters
    ----------
    points:
        Points inside one segment. The table may contain ``z`` and
        ``intensity`` columns.
    metrics:
        Feature groups to calculate. Supported values are ``"height"``,
        ``"intensity"``, and ``"density"``.
    area:
        Segment area in CRS units squared. Required for density.
    height_column:
        Column used for height metrics.
    intensity_column:
        Column used for intensity metrics.
    """
    metric_set = set(metrics)
    unknown = metric_set.difference(DEFAULT_METRICS)
    if unknown:
        raise ValueError(f"Unknown point-cloud metrics: {sorted(unknown)}")

    features: dict[str, float] = {"pc_point_count": int(len(points))}

    if "density" in metric_set:
        if area is None or area <= 0:
            features["pc_density"] = np.nan
        else:
            features["pc_density"] = float(len(points) / area)

    if "height" in metric_set:
        features.update(_numeric_series_features(points, height_column, "pc_z", HEIGHT_PERCENTILES))

    if "intensity" in metric_set:
        features.update(
            _numeric_series_features(points, intensity_column, "pc_intensity", INTENSITY_PERCENTILES)
        )

    return features

empty_pointcloud_features(*, metrics=DEFAULT_METRICS)

Return a feature dictionary filled with empty values.

Source code in obia/pointcloud/features.py 
65
66
67
68
69
70
def empty_pointcloud_features(
    *,
    metrics: Iterable[str] = DEFAULT_METRICS,
) -> dict[str, float]:
    """Return a feature dictionary filled with empty values."""
    return calculate_pointcloud_features(pd.DataFrame(), metrics=metrics, area=np.nan)

obia.pointcloud.segment_features

obia.pointcloud.segment_features

Join point-cloud feature groups onto segment objects.

add_pointcloud_features(segments, pointcloud, *, metrics=DEFAULT_METRICS, crs=None, reader='auto', predicate='within')

Add point-cloud features to segment objects.

Parameters

segments: Segment polygons. The output preserves the same rows and appends pc_* feature columns. pointcloud: LAS/LAZ path or in-memory point cloud. Path inputs require optional pointcloud dependencies. metrics: Feature groups to calculate: "height", "intensity", and "density". crs: CRS assigned to point clouds that do not already have one. reader: Backend for path inputs: "auto", "pdal", or "pyforestscan". predicate: Spatial join predicate used to attach points to polygons.

Source code in obia/pointcloud/segment_features.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def add_pointcloud_features(
    segments: gpd.GeoDataFrame,
    pointcloud: str | Path | np.ndarray | pd.DataFrame | gpd.GeoDataFrame,
    *,
    metrics: Iterable[str] = DEFAULT_METRICS,
    crs=None,
    reader: str = "auto",
    predicate: str = "within",
) -> gpd.GeoDataFrame:
    """Add point-cloud features to segment objects.

    Parameters
    ----------
    segments:
        Segment polygons. The output preserves the same rows and appends
        ``pc_*`` feature columns.
    pointcloud:
        LAS/LAZ path or in-memory point cloud. Path inputs require optional
        pointcloud dependencies.
    metrics:
        Feature groups to calculate: ``"height"``, ``"intensity"``, and
        ``"density"``.
    crs:
        CRS assigned to point clouds that do not already have one.
    reader:
        Backend for path inputs: ``"auto"``, ``"pdal"``, or ``"pyforestscan"``.
    predicate:
        Spatial join predicate used to attach points to polygons.
    """
    if segments.empty:
        return segments.copy()
    if segments.geometry.name is None:
        raise ValueError("segments must have an active geometry column")

    points = read_pointcloud(pointcloud, crs=crs or segments.crs, reader=reader)
    if points.empty:
        return _append_empty_features(segments, metrics)

    if segments.crs is not None:
        if points.crs is None:
            points = points.set_crs(segments.crs)
        elif points.crs != segments.crs:
            points = points.to_crs(segments.crs)

    segment_index_name = "__obia_segment_index"
    point_index_name = "__obia_point_index"

    segment_table = segments[[segments.geometry.name]].copy()
    segment_table[segment_index_name] = segments.index

    point_table = points.copy()
    point_table[point_index_name] = np.arange(len(point_table))

    joined = gpd.sjoin(
        point_table,
        segment_table,
        how="inner",
        predicate=predicate,
    )

    features_by_segment = {}
    for segment_index, group in joined.groupby(segment_index_name):
        points_in_segment = point_table.iloc[group[point_index_name].to_numpy()]
        area = float(segments.loc[segment_index].geometry.area)
        features_by_segment[segment_index] = calculate_pointcloud_features(
            points_in_segment,
            metrics=metrics,
            area=area,
        )

    empty_features = calculate_pointcloud_features(pd.DataFrame(), metrics=metrics, area=np.nan)
    output = segments.copy()
    for column in empty_features:
        output[column] = np.nan

    for segment_index, features in features_by_segment.items():
        for column, value in features.items():
            output.loc[segment_index, column] = value

    output["pc_point_count"] = output["pc_point_count"].fillna(0).astype(int)
    return output