Repository: single-cell-data/TileDB-SOMA (language: python)

Issue:
[spatial] Add PointCloud class to main branch
[[sc-50281]](https://app.shortcut.com/tiledb-inc/story/50281/add-pointcloud-class-to-main-branch)



Discussion / hints:

Fix patch:
diff --git a/apis/python/setup.py b/apis/python/setup.py
index 5e9102f256..11c373a70c 100644
--- a/apis/python/setup.py
+++ b/apis/python/setup.py
@@ -307,6 +307,7 @@ def run(self):
                 "src/tiledbsoma/soma_array.cc",
                 "src/tiledbsoma/soma_object.cc",
                 "src/tiledbsoma/soma_dataframe.cc",
+                "src/tiledbsoma/soma_point_cloud_dataframe.cc",
                 "src/tiledbsoma/soma_dense_ndarray.cc",
                 "src/tiledbsoma/soma_sparse_ndarray.cc",
                 "src/tiledbsoma/soma_group.cc",
diff --git a/apis/python/src/tiledbsoma/__init__.py b/apis/python/src/tiledbsoma/__init__.py
index 308bf19012..569a785e44 100644
--- a/apis/python/src/tiledbsoma/__init__.py
+++ b/apis/python/src/tiledbsoma/__init__.py
@@ -135,7 +135,17 @@
         # Otherwise try loading by name only.
         ctypes.CDLL(libtiledbsoma_name)
 
-from somacore import AxisColumnNames, AxisQuery, ExperimentAxisQuery
+from somacore import (
+    Axis,
+    CoordinateSpace,
+    AffineTransform,
+    ScaleTransform,
+    IdentityTransform,
+    UniformScaleTransform,
+    AxisColumnNames,
+    AxisQuery,
+    ExperimentAxisQuery,
+)
 from somacore.options import ResultOrder
 
 # This is important since we need to do the above dll/dylib/so business
@@ -165,6 +175,7 @@
 )
 from ._indexer import IntIndexer, tiledbsoma_build_index
 from ._measurement import Measurement
+from ._point_cloud_dataframe import PointCloudDataFrame
 from ._sparse_nd_array import SparseNDArray, SparseNDArrayRead
 from .options import SOMATileDBContext, TileDBCreateOptions, TileDBWriteOptions
 from .pytiledbsoma import (
@@ -182,9 +193,11 @@
 
 __all__ = [
     "AlreadyExistsError",
+    "Axis",
     "AxisColumnNames",
     "AxisQuery",
     "Collection",
+    "CoordinateSpace",
     "DataFrame",
     "DenseNDArray",
     "DoesNotExistError",
@@ -198,6 +211,7 @@
     "Measurement",
     "NotCreateableError",
     "open",
+    "PointCloudDataFrame",
     "ResultOrder",
     "show_package_versions",
     "SOMA_JOINID",
diff --git a/apis/python/src/tiledbsoma/_constants.py b/apis/python/src/tiledbsoma/_constants.py
index 724b4f91dd..0e29438e17 100644
--- a/apis/python/src/tiledbsoma/_constants.py
+++ b/apis/python/src/tiledbsoma/_constants.py
@@ -1,5 +1,5 @@
-# Copyright (c) 2021-2023 The Chan Zuckerberg Initiative Foundation
-# Copyright (c) 2021-2023 TileDB, Inc.
+# Copyright (c) 2021-2024 The Chan Zuckerberg Initiative Foundation
+# Copyright (c) 2021-2024 TileDB, Inc.
 #
 # Licensed under the MIT License.
 
@@ -8,6 +8,8 @@
 
 SOMA_JOINID = "soma_joinid"
 SOMA_GEOMETRY = "soma_geometry"
+SOMA_COORDINATE_SPACE_METADATA_KEY = "soma_coordinate_space"
+SOMA_MULTISCALE_IMAGE_SCHAME = "soma_multiscale_image_schema"
 SOMA_OBJECT_TYPE_METADATA_KEY = "soma_object_type"
 SOMA_ENCODING_VERSION_METADATA_KEY = "soma_encoding_version"
 SOMA_ENCODING_VERSION = "1.1.0"
diff --git a/apis/python/src/tiledbsoma/_factory.py b/apis/python/src/tiledbsoma/_factory.py
index 040c3bf330..a04befdc00 100644
--- a/apis/python/src/tiledbsoma/_factory.py
+++ b/apis/python/src/tiledbsoma/_factory.py
@@ -28,6 +28,7 @@
     _dense_nd_array,
     _experiment,
     _measurement,
+    _point_cloud_dataframe,
     _soma_object,
     _sparse_nd_array,
     _tdb_handles,
@@ -216,6 +217,7 @@ def _type_name_to_cls(type_name: str) -> Type[AnySOMAObject]:
             _experiment.Experiment,
             _measurement.Measurement,
             _sparse_nd_array.SparseNDArray,
+            _point_cloud_dataframe.PointCloudDataFrame,
         )
     }
     try:
diff --git a/apis/python/src/tiledbsoma/_point_cloud_dataframe.py b/apis/python/src/tiledbsoma/_point_cloud_dataframe.py
index 660e5b9ea7..2a621971d0 100644
--- a/apis/python/src/tiledbsoma/_point_cloud_dataframe.py
+++ b/apis/python/src/tiledbsoma/_point_cloud_dataframe.py
@@ -3,40 +3,67 @@
 #
 # Licensed under the MIT License.
 """
-Implementation of a SOMA Point Cloud
+Implementation of a SOMA Point Cloud DataFrame
 """
 
 import warnings
-from typing import Any, Optional, Sequence, Tuple, Union
+from typing import Any, Optional, Sequence, Tuple, Union, cast
 
 import pyarrow as pa
 import somacore
-from somacore import CoordinateSpace, CoordinateTransform, options
+from somacore import Axis, CoordinateSpace, CoordinateTransform, options
 from typing_extensions import Self
 
-from ._constants import SOMA_JOINID, SPATIAL_DISCLAIMER
-from ._dataframe import Domain
+from . import _arrow_types, _util
+from . import pytiledbsoma as clib
+from ._constants import (
+    SOMA_COORDINATE_SPACE_METADATA_KEY,
+    SOMA_JOINID,
+    SPATIAL_DISCLAIMER,
+)
+from ._dataframe import (
+    Domain,
+    _canonicalize_schema,
+    _fill_out_slot_soma_domain,
+    _find_extent_for_domain,
+    _revise_domain_for_extent,
+)
+from ._exception import SOMAError, map_exception_for_create
+from ._flags import NEW_SHAPE_FEATURE_FLAG_ENABLED
+from ._query_condition import QueryCondition
 from ._read_iters import TableReadIter
+from ._spatial_dataframe import SpatialDataFrame
+from ._spatial_util import (
+    coordinate_space_from_json,
+    coordinate_space_to_json,
+    process_spatial_df_region,
+)
+from ._tdb_handles import PointCloudDataFrameWrapper
 from ._types import OpenTimestamp
 from .options import SOMATileDBContext
+from .options._soma_tiledb_context import _validate_soma_tiledb_context
+from .options._tiledb_create_write_options import (
+    TileDBCreateOptions,
+    TileDBWriteOptions,
+)
 
 _UNBATCHED = options.BatchSize()
 
 
-class PointCloudDataFrame(somacore.PointCloudDataFrame):
+class PointCloudDataFrame(SpatialDataFrame, somacore.PointCloudDataFrame):
     """A specialized SOMA DataFrame for storing collections of points in
     multi-dimensional space.
 
-    The ``PointCloudDataFrame`` class is designed to efficiently store and query point
-    data, where each point is represented by coordinates in one or more spatial
-    dimensions (e.g., x, y, z) and may have additional columns for associated
-    attributes.
+    The ``PointCloudDataFrame`` class is designed to efficiently store and query point data,
+    where each point is represented by coordinates in one or more spatial dimensions
+    (e.g., x, y, z) and may have additional columns for associated attributes.
 
     Lifecycle:
         Experimental.
     """
 
-    __slots__ = ()
+    __slots__ = ("_coord_space",)
+    _wrapper_type = PointCloudDataFrameWrapper
 
     @classmethod
     def create(
@@ -87,10 +114,179 @@ def create(
             Experimental.
         """
         warnings.warn(SPATIAL_DISCLAIMER)
-        raise NotImplementedError()
+
+        axis_dtype: Optional[pa.DataType] = None
+        for column_name in axis_names:
+            if column_name not in index_column_names:
+                raise ValueError(f"Spatial column '{column_name}' must an index column")
+            # Check axis column type is valid and all axis columns have the same type.
+            if axis_dtype is None:
+                axis_dtype = schema.field(column_name).type
+                if not (
+                    pa.types.is_integer(axis_dtype) or pa.types.is_floating(axis_dtype)
+                ):
+                    raise ValueError(
+                        f"Spatial column '{column_name}' must have an integer or "
+                        f"floating-point type. Column type is {axis_dtype!r}"
+                    )
+            else:
+                column_dtype = schema.field(column_name).type
+                if column_dtype != axis_dtype:
+                    raise ValueError("All spatial axes must have the same datatype.")
+
+        # mypy false positive https://github.com/python/mypy/issues/5313
+        coord_space = CoordinateSpace(
+            tuple(Axis(axis_name) for axis_name in axis_names)  # type: ignore
+        )
+        context = _validate_soma_tiledb_context(context)
+        schema = _canonicalize_schema(schema, index_column_names)
+
+        # SOMA-to-core mappings:
+        #
+        # Before the current-domain feature was enabled (possible after core 2.25):
+        #
+        # * SOMA domain <-> core domain, AKA "max domain" which is a name we'll use for clarity
+        # * core current domain did not exist
+        #
+        # After the current-domain feature was enabled:
+        #
+        # * SOMA max_domain <-> core domain
+        # * SOMA domain <-> core current domain
+        #
+        # As far as the user is concerned, the SOMA-level domain is the only
+        # thing they see and care about. Before 2.25 support, it was immutable
+        # (since it was implemented by core domain). After 2.25 support, it is
+        # mutable/up-resizeable (since it is implemented by core current domain).
+
+        # At this point shift from API terminology "domain" to specifying a soma_ or core_
+        # prefix for these variables. This is crucial to avoid developer confusion.
+        soma_domain = domain
+        domain = None
+
+        if soma_domain is None:
+            soma_domain = tuple(None for _ in index_column_names)
+        else:
+            ndom = len(soma_domain)
+            nidx = len(index_column_names)
+            if ndom != nidx:
+                raise ValueError(
+                    f"if domain is specified, it must have the same length as index_column_names; got {ndom} != {nidx}"
+                )
+
+        index_column_schema = []
+        index_column_data = {}
+
+        for index_column_name, slot_soma_domain in zip(index_column_names, soma_domain):
+            pa_field = schema.field(index_column_name)
+            dtype = _arrow_types.tiledb_type_from_arrow_type(
+                pa_field.type, is_indexed_column=True
+            )
+
+            (slot_core_current_domain, saturated_cd) = _fill_out_slot_soma_domain(
+                slot_soma_domain, index_column_name, pa_field.type, dtype
+            )
+            (slot_core_max_domain, saturated_md) = _fill_out_slot_soma_domain(
+                None, index_column_name, pa_field.type, dtype
+            )
+
+            extent = _find_extent_for_domain(
+                index_column_name,
+                TileDBCreateOptions.from_platform_config(platform_config),
+                dtype,
+                slot_core_current_domain,
+            )
+
+            # Necessary to avoid core array-creation error "Reduce domain max by
+            # 1 tile extent to allow for expansion."
+            slot_core_current_domain = _revise_domain_for_extent(
+                slot_core_current_domain, extent, saturated_cd
+            )
+            slot_core_max_domain = _revise_domain_for_extent(
+                slot_core_max_domain, extent, saturated_md
+            )
+
+            # Here is our Arrow data API for communicating schema info between
+            # Python/R and C++ libtiledbsoma:
+            #
+            # [0] core max domain lo
+            # [1] core max domain hi
+            # [2] core extent parameter
+            # If present, these next two signal to use the current-domain feature:
+            # [3] core current domain lo
+            # [4] core current domain hi
+
+            index_column_schema.append(pa_field)
+            if NEW_SHAPE_FEATURE_FLAG_ENABLED:
+
+                index_column_data[pa_field.name] = [
+                    *slot_core_max_domain,
+                    extent,
+                    *slot_core_current_domain,
+                ]
+
+            else:
+                index_column_data[pa_field.name] = [*slot_core_current_domain, extent]
+
+        index_column_info = pa.RecordBatch.from_pydict(
+            index_column_data, schema=pa.schema(index_column_schema)
+        )
+
+        plt_cfg = _util.build_clib_platform_config(platform_config)
+        timestamp_ms = context._open_timestamp_ms(tiledb_timestamp)
+        try:
+            clib.SOMAPointCloudDataFrame.create(
+                uri,
+                schema=schema,
+                index_column_info=index_column_info,
+                ctx=context.native_context,
+                platform_config=plt_cfg,
+                timestamp=(0, timestamp_ms),
+            )
+        except SOMAError as e:
+            raise map_exception_for_create(e, uri) from None
+
+        handle = cls._wrapper_type.open(uri, "w", context, tiledb_timestamp)
+        handle.meta[SOMA_COORDINATE_SPACE_METADATA_KEY] = coordinate_space_to_json(
+            coord_space
+        )
+        return cls(
+            handle,
+            _dont_call_this_use_create_or_open_instead="tiledbsoma-internal-code",
+        )
+
+    def __init__(
+        self,
+        handle: PointCloudDataFrameWrapper,
+        **kwargs: Any,
+    ):
+        super().__init__(handle, **kwargs)
+
+        # Get and validate coordinate space.
+        try:
+            coord_space = self.metadata[SOMA_COORDINATE_SPACE_METADATA_KEY]
+        except KeyError as ke:
+            raise SOMAError("Missing coordinate space metadata") from ke
+        self._coord_space = coordinate_space_from_json(coord_space)
+        for name in self._coord_space.axis_names:
+            if name not in self.index_column_names:
+                raise SOMAError(
+                    f"Point cloud dataframe axis '{name}' does not match any of the "
+                    f"index column names."
+                )
 
     # Data operations
 
+    def __len__(self) -> int:
+        """Returns the number of rows in the point cound dataframe."""
+        return self.count
+
+    @property
+    def count(self) -> int:
+        """Returns the number of rows in the dataframe."""
+        self._check_open_read()
+        # if is it in read open mode, then it is a PointCloudDataFrameWrapper
+        return cast(PointCloudDataFrameWrapper, self._handle).count
+
     def read(
         self,
         coords: options.SparseDFCoords = (),
@@ -123,7 +319,34 @@ def read(
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        del batch_size  # Currently unused.
+        _util.check_unpartitioned(partitions)
+        self._check_open_read()
+
+        handle = self._handle._handle
+
+        context = handle.context()
+        if platform_config is not None:
+            config = context.tiledb_config.copy()
+            config.update(platform_config)
+            context = clib.SOMAContext(config)
+
+        sr = clib.SOMAPointCloudDataFrame.open(
+            uri=handle.uri,
+            mode=clib.OpenMode.read,
+            context=context,
+            column_names=column_names or [],
+            result_order=_util.to_clib_result_order(result_order),
+            timestamp=handle.timestamp and (0, handle.timestamp),
+        )
+
+        if value_filter is not None:
+            sr.set_condition(QueryCondition(value_filter), handle.schema)
+
+        self._set_reader_coords(sr, coords)
+
+        # # TODO: batch_size
+        return TableReadIter(sr)
 
     def read_spatial_region(
         self,
@@ -141,7 +364,6 @@ def read_spatial_region(
         """Reads data intersecting an user-defined region of space into a
         :class:`SpatialRead` with data in Arrow tables.
 
-
         Args:
             region: The region to query. May be a box in the form
                 [x_min, y_min, x_max, y_max] (for 2D images), a box in the form
@@ -171,7 +393,60 @@ def read_spatial_region(
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        # Set/check transform and region coordinate space.
+        if region_transform is None:
+            region_transform = somacore.IdentityTransform(
+                self.axis_names, self.axis_names
+            )
+            if region_coord_space is not None:
+                raise ValueError(
+                    "Cannot specify the output coordinate space when region transform i"
+                    "is ``None``."
+                )
+            region_coord_space = self._coord_space
+        else:
+            if region_coord_space is None:
+                # mypy false positive https://github.com/python/mypy/issues/5313
+                region_coord_space = CoordinateSpace(
+                    tuple(Axis(axis_name) for axis_name in region_transform.input_axes)  # type: ignore
+                )
+            elif region_transform.input_axes != region_coord_space.axis_names:
+                raise ValueError(
+                    f"The input axes '{region_transform.input_axes}' of the region "
+                    f"transform must match the axes '{region_coord_space.axis_names}' "
+                    f"of the coordinate space the requested region is defined in."
+                )
+            if region_transform.output_axes != self._coord_space.axis_names:
+                raise ValueError(
+                    f"The output axes of '{region_transform.output_axes}' of the "
+                    f"transform must match the axes '{self._coord_space.axis_names}' "
+                    f"of the coordinate space of this point cloud dataframe."
+                )
+
+        # Process the user provided region.
+        coords, data_region, inv_transform = process_spatial_df_region(
+            region,
+            region_transform,
+            dict(),  #  Move index value_filters into this dict to optimize queries
+            self._tiledb_dim_names(),
+            self._coord_space.axis_names,
+            self._handle.schema,
+        )
+
+        return somacore.SpatialRead(
+            self.read(
+                coords,
+                column_names,
+                result_order=result_order,
+                value_filter=value_filter,
+                batch_size=batch_size,
+                partitions=partitions,
+                platform_config=platform_config,
+            ),
+            self.coordinate_space,
+            region_coord_space,
+            inv_transform,
+        )
 
     def write(
         self,
@@ -194,45 +469,57 @@ def write(
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        _util.check_type("values", values, (pa.Table,))
 
-    # Metadata operations
+        write_options: Union[TileDBCreateOptions, TileDBWriteOptions]
+        sort_coords = None
+        if isinstance(platform_config, TileDBCreateOptions):
+            raise ValueError(
+                "As of TileDB-SOMA 1.13, the write method takes "
+                "TileDBWriteOptions instead of TileDBCreateOptions"
+            )
+        write_options = TileDBWriteOptions.from_platform_config(platform_config)
+        sort_coords = write_options.sort_coords
 
-    @property
-    def schema(self) -> pa.Schema:
-        """The schema of the data in this dataframe.
+        clib_dataframe = self._handle._handle
 
-        Lifecycle:
-            Experimental.
-        """
-        raise NotImplementedError()
+        for batch in values.to_batches():
+            clib_dataframe.write(batch, sort_coords or False)
 
-    @property
-    def index_column_names(self) -> Tuple[str, ...]:
-        """The names of the index (dimension) columns.
+        if write_options.consolidate_and_vacuum:
+            clib_dataframe.consolidate_and_vacuum()
 
-        Lifecycle:
-            Experimental.
-        """
-        raise NotImplementedError()
+        return self
+
+    # Metadata operations
 
     @property
     def coordinate_space(self) -> CoordinateSpace:
-        """Coordinate space for this point cloud.
+        """Coordinate space for this point cloud dataframe.
 
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return self._coord_space
 
     @coordinate_space.setter
     def coordinate_space(self, value: CoordinateSpace) -> None:
-        """Coordinate space for this point cloud.
+        """Coordinate space for this point cloud dataframe.
 
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        if self._coord_space is not None:
+            if value.axis_names != self._coord_space.axis_names:
+                raise ValueError(
+                    f"Cannot change axis names of a point cloud dataframe. Existing "
+                    f"axis names are {self._coord_space.axis_names}. New coordinate "
+                    f"space has axis names {self._coord_space.axis_names}."
+                )
+        self.metadata[SOMA_COORDINATE_SPACE_METADATA_KEY] = coordinate_space_to_json(
+            value
+        )
+        self._coord_space = value
 
     @property
     def axis_names(self) -> Tuple[str, ...]:
@@ -241,16 +528,4 @@ def axis_names(self) -> Tuple[str, ...]:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
-
-    @property
-    def domain(self) -> Tuple[Tuple[Any, Any], ...]:
-        """The allowable range of values in each index column.
-
-        Returns: a tuple of minimum and maximum values, inclusive,
-            storable on each index column of the dataframe.
-
-        Lifecycle:
-            Experimental.
-        """
-        raise NotImplementedError()
+        return self._coord_space.axis_names
diff --git a/apis/python/src/tiledbsoma/_scene.py b/apis/python/src/tiledbsoma/_scene.py
index 5e10db2cc5..e4b5b254be 100644
--- a/apis/python/src/tiledbsoma/_scene.py
+++ b/apis/python/src/tiledbsoma/_scene.py
@@ -110,7 +110,7 @@ def add_multiscale_image(
         axis_types: Sequence[str] = ("channel", "height", "width"),
     ) -> MultiscaleImage:
         """Adds a ``MultiscaleImage`` to the scene and sets a coordinate transform
-        between the scene and the dataframe.
+        between the scene and the multiscale image.
 
         Parameters are as in :meth:`spatial.MultiscaleImage.create`.
         See :meth:`add_new_collection` for details about child URIs.
@@ -176,7 +176,7 @@ def set_transform_to_geometry_dataframe(
         to set a transformation for geometry dataframe named  "transcripts" in the
         "var/RNA" collection::
 
-            scene.set_transfrom_for_geometry_dataframe(
+            scene.set_transfrom_to_geometry_dataframe(
                 'transcripts', transform, subcollection=['var', 'RNA'],
             )
 
@@ -226,7 +226,7 @@ def set_transform_to_multiscale_image(
         """
         raise NotImplementedError()
 
-    def set_transform_to_point_cloud(
+    def set_transform_to_point_cloud_dataframe(
         self,
         key: str,
         transform: CoordinateTransform,
@@ -235,14 +235,14 @@ def set_transform_to_point_cloud(
         coordinate_space: Optional[CoordinateSpace] = None,
     ) -> PointCloudDataFrame:
         """Adds the coordinate transform for the scene coordinate space to
-        a point cloud stored in the scene.
+        a point cloud dataframe stored in the scene.
 
         If the subcollection the point cloud is inside of is more than one
         layer deep, the input should be provided as a sequence of names. For example,
         to set a transform for  a point named `transcripts` in the `var/RNA`
         collection::
 
-            scene.set_transformation_for_point_cloud(
+            scene.set_transformation_to_point_cloud_dataframe(
                 'transcripts', transform, subcollection=['var', 'RNA'],
             )
 
@@ -305,19 +305,19 @@ def get_transform_from_multiscale_image(
         """
         raise NotImplementedError()
 
-    def get_transform_from_point_cloud(
+    def get_transform_from_point_cloud_dataframe(
         self, key: str, *, subcollection: str = "obsl"
     ) -> CoordinateTransform:
-        """Returns the coordinate transformation from the requested point cloud to
-        the scene.
+        """Returns the coordinate transformation from the requested point cloud
+        dataframe to the scene.
 
         Args:
-            key: The name of the point cloud.
+            key: The name of the point cloud dataframe.
             subcollection: The name, or sequence of names, of the subcollection the
-                point cloud is stored in. Defaults to ``'obsl'``.
+                dataframe is stored in. Defaults to ``'obsl'``.
 
         Returns:
-            Coordinate transform from the scene to the point cloud.
+            Coordinate transform from the scene to the point cloud dataframe.
 
         Lifecycle: experimental
         """
@@ -366,7 +366,7 @@ def get_transform_to_multiscale_image(
         """
         raise NotImplementedError()
 
-    def get_transform_to_point_cloud(
+    def get_transform_to_point_cloud_dataframe(
         self, key: str, *, subcollection: str = "obsl"
     ) -> CoordinateTransform:
         """Returns the coordinate transformation from the scene to a requested
@@ -375,7 +375,7 @@ def get_transform_to_point_cloud(
         Args:
             key: The name of the point cloud.
             subcollection: The name, or sequence of names, of the subcollection the
-                point cloud is stored in. Defaults to ``'obsl'``.
+                point cloud dataframe is stored in. Defaults to ``'obsl'``.
 
         Returns:
             Coordinate transform from the scene to the requested point cloud.
diff --git a/apis/python/src/tiledbsoma/_spatial_dataframe.py b/apis/python/src/tiledbsoma/_spatial_dataframe.py
new file mode 100644
index 0000000000..571b5d116c
--- /dev/null
+++ b/apis/python/src/tiledbsoma/_spatial_dataframe.py
@@ -0,0 +1,309 @@
+# Copyright (c) 2024 The Chan Zuckerberg Initiative Foundation
+# Copyright (c) 2024 TileDB, Inc.
+#
+# Licensed under the MIT License.
+
+"""
+Implementation of a base class shared between GeometryDataFrame and PointCloudDataFrame
+"""
+from typing import Any, Optional, Sequence, Tuple, Type, Union
+
+import numpy as np
+import pyarrow as pa
+import somacore
+from somacore import CoordinateSpace, CoordinateTransform, options
+from typing_extensions import Self
+
+from . import _util
+from . import pytiledbsoma as clib
+from ._read_iters import TableReadIter
+from ._soma_array import SOMAArray
+from ._types import Slice, is_slice_of
+
+_UNBATCHED = options.BatchSize()
+
+
+class SpatialDataFrame(SOMAArray):
+
+    __slots__ = ()
+
+    def keys(self) -> Tuple[str, ...]:
+        """Returns the names of the columns when read back as a spatial dataframe.
+
+        Examples:
+            >>> with tiledbsoma.open("a_dataframe") as soma_df:
+            ...     k = soma_df.keys()
+            ...
+            >>> k
+            ('soma_joinid', 'col1')
+
+        Lifecycle:
+            Experimental.
+        """
+        return self._tiledb_array_keys()
+
+    @property
+    def index_column_names(self) -> Tuple[str, ...]:
+        """Returns index (dimension) column names.
+
+        Lifecycle:
+            Experimental.
+        """
+        return self._tiledb_dim_names()
+
+    @property
+    def axis_names(self) -> Tuple[str, ...]:
+        """The names of the axes of the coordinate space the data is defined on.
+
+        Lifecycle: experimental
+        """
+        raise NotImplementedError("Must be implemented by the child class")
+
+    @property
+    def domain(self) -> Tuple[Tuple[Any, Any], ...]:
+        """Returns a tuple of minimum and maximum values, inclusive, storable
+        on each index column of the dataframe.
+
+        Lifecycle:
+            Experimental.
+        """
+        return self._domain()
+
+    def read(
+        self,
+        coords: options.SparseDFCoords = (),
+        column_names: Optional[Sequence[str]] = None,
+        *,
+        result_order: options.ResultOrderStr = options.ResultOrder.AUTO,
+        value_filter: Optional[str] = None,
+        batch_size: options.BatchSize = _UNBATCHED,
+        partitions: Optional[options.ReadPartitions] = None,
+        platform_config: Optional[options.PlatformConfig] = None,
+    ) -> TableReadIter:
+        """Reads a user-defined slice of data into Arrow tables.
+
+        Args:
+            coords: for each index dimension, which rows to read.
+                Defaults to ``()``, meaning no constraint -- all IDs.
+            column_names: the named columns to read and return.
+                Defaults to ``None``, meaning no constraint -- all column names.
+            partitions: If present, specifies that this is part of
+                a partitioned read, and which part of the data to include.
+            result_order: the order to return results, specified as a
+                :class:`~options.ResultOrder` or its string value.
+            value_filter: an optional value filter to apply to the results.
+                The default of ``None`` represents no filter. Value filter
+                syntax is implementation-defined; see the documentation
+                for the particular SOMA implementation for details.
+        Returns:
+            A :class:`ReadIter` of :class:`pa.Table`s.
+
+        Lifecycle:
+            Experimental.
+        """
+        raise NotImplementedError("must be implemented by child class")
+
+    def read_spatial_region(
+        self,
+        region: Optional[options.SpatialRegion] = None,
+        column_names: Optional[Sequence[str]] = None,
+        *,
+        region_transform: Optional[CoordinateTransform] = None,
+        region_coord_space: Optional[CoordinateSpace] = None,
+        batch_size: options.BatchSize = _UNBATCHED,
+        partitions: Optional[options.ReadPartitions] = None,
+        result_order: options.ResultOrderStr = options.ResultOrder.AUTO,
+        value_filter: Optional[str] = None,
+        platform_config: Optional[options.PlatformConfig] = None,
+    ) -> somacore.SpatialRead[somacore.ReadIter[pa.Table]]:
+        """Reads data intersecting an user-defined region of space into a
+        :class:`SpatialRead` with data in Arrow tables.
+
+
+        Args:
+            region: The region to query. May be a box in the form
+                [x_min, y_min, x_max, y_max] (for 2D images), a box in the form
+                [x_min, y_min, z_min, x_max, y_max, z_max] (for 3D images), or
+                a shapely Geometry.
+            column_names: The named columns to read and return.
+                Defaults to ``None``, meaning no constraint -- all column names.
+            region_transform: An optional coordinate transform from the read region to
+                the coordinate system of the spatial dataframe.
+                Defaults to ``None``, meaning an identity transform.
+            region_coord_space: An optional coordinate space for the region being read.
+                Defaults to ``None``, coordinate space will be inferred from transform.
+            batch_size: The size of batched reads.
+                Defaults to `unbatched`.
+            partitions: If present, specifies that this is part of a partitioned read,
+                and which part of the data to include.
+            result_order: the order to return results, specified as a
+                :class:`~options.ResultOrder` or its string value.
+            value_filter: an optional value filter to apply to the results.
+                The default of ``None`` represents no filter. Value filter
+                syntax is implementation-defined; see the documentation
+                for the particular SOMA implementation for details.
+
+        Returns:
+            A :class:`SpatialRead` with :class:`ReadIter` of :class:`pa.Table`s data.
+
+        Lifecycle:
+            Experimental.
+        """
+        raise NotImplementedError("must be implemented by child class")
+
+    def write(
+        self,
+        values: Union[pa.RecordBatch, pa.Table],
+        *,
+        platform_config: Optional[options.PlatformConfig] = None,
+    ) -> Self:
+        """Writes the data from an Arrow table to the persistent object.
+
+        As duplicate index values are not allowed, index values already present
+        in the object are overwritten and new index values are added.
+
+        Args:
+            values: An Arrow table containing all columns, including
+                the index columns. The schema for the values must match
+                the schema for the ``DataFrame``.
+
+        Returns: ``self``, to enable method chaining.
+
+        Lifecycle:
+            Experimental.
+        """
+        raise NotImplementedError("must be implemented by child class")
+
+    def _set_reader_coord(
+        self,
+        sr: clib.SOMAArray,
+        dim_idx: int,
+        dim: pa.Field,
+        coord: object,
+    ) -> bool:
+        if coord is None:
+            return True  # No constraint; select all in this dimension
+
+        if isinstance(coord, (str, bytes)):
+            sr.set_dim_points_string_or_bytes(dim.name, [coord])
+            return True
+
+        if isinstance(coord, (pa.Array, pa.ChunkedArray)):
+            # sr.set_dim_points_arrow does type disambiguation based on array schema -- so we do
+            # not.
+            sr.set_dim_points_arrow(dim.name, coord)
+            return True
+
+        if isinstance(coord, (Sequence, np.ndarray)):
+            if self._set_reader_coord_by_py_seq_or_np_array(sr, dim_idx, dim, coord):
+                return True
+
+        if isinstance(coord, slice):
+            _util.validate_slice(coord)
+            if coord.start is None and coord.stop is None:
+                return True
+
+        if isinstance(coord, slice):
+            _util.validate_slice(coord)
+            if self._set_reader_coord_by_numeric_slice(sr, dim_idx, dim, coord):
+                return True
+
+        domain = self.domain[dim_idx]
+
+        # Note: slice(None, None) matches the is_slice_of part, unless we also check the dim-type
+        # part.
+        if (
+            is_slice_of(coord, str) or is_slice_of(coord, bytes)
+        ) and _util.pa_types_is_string_or_bytes(dim.type):
+            _util.validate_slice(coord)
+            # Figure out which one.
+            dim_type: Union[Type[str], Type[bytes]] = type(domain[0])
+            # A ``None`` or empty start is always equivalent to empty str/bytes.
+            start = coord.start or dim_type()
+            if coord.stop is None:
+                # There's no way to specify "to infinity" for strings.
+                # We have to get the nonempty domain and use that as the end.
+                ned = self._handle.non_empty_domain()
+                _, stop = ned[dim_idx]
+            else:
+                stop = coord.stop
+            sr.set_dim_ranges_string_or_bytes(dim.name, [(start, stop)])
+            return True
+
+        # Note: slice(None, None) matches the is_slice_of part, unless we also check the dim-type
+        # part.
+        if is_slice_of(coord, np.datetime64) and pa.types.is_timestamp(dim.type):
+            _util.validate_slice(coord)
+            # These timestamp types are stored in Arrow as well as TileDB as 64-bit integers (with
+            # distinguishing metadata of course). For purposes of the query logic they're just
+            # int64.
+            istart = coord.start or domain[0]
+            istart = int(istart.astype("int64"))
+            istop = coord.stop or domain[1]
+            istop = int(istop.astype("int64"))
+            sr.set_dim_ranges_int64(dim.name, [(istart, istop)])
+            return True
+
+        if super()._set_reader_coord(sr, dim_idx, dim, coord):
+            return True
+
+        return False
+
+    def _set_reader_coord_by_py_seq_or_np_array(
+        self,
+        sr: clib.SOMAArray,
+        dim_idx: int,
+        dim: pa.Field,
+        coord: object,
+    ) -> bool:
+        if isinstance(coord, np.ndarray):
+            if coord.ndim != 1:
+                raise ValueError(
+                    f"only 1D numpy arrays may be used to index; got {coord.ndim}"
+                )
+
+        try:
+            set_dim_points = getattr(sr, f"set_dim_points_{dim.type}")
+        except AttributeError:
+            # We have to handle this type specially below
+            pass
+        else:
+            set_dim_points(dim.name, coord)
+            return True
+
+        if _util.pa_types_is_string_or_bytes(dim.type):
+            sr.set_dim_points_string_or_bytes(dim.name, coord)
+            return True
+        elif pa.types.is_timestamp(dim.type):
+            if not isinstance(coord, (tuple, list, np.ndarray)):
+                raise ValueError(
+                    f"unhandled coord type {type(coord)} for index column named {dim.name}"
+                )
+            icoord = [
+                int(e.astype("int64")) if isinstance(e, np.datetime64) else e
+                for e in coord
+            ]
+            sr.set_dim_points_int64(dim.name, icoord)
+            return True
+
+        # TODO: bool
+
+        raise ValueError(f"unhandled type {dim.type} for index column named {dim.name}")
+
+    def _set_reader_coord_by_numeric_slice(
+        self, sr: clib.SOMAArray, dim_idx: int, dim: pa.Field, coord: Slice[Any]
+    ) -> bool:
+        try:
+            lo_hi = _util.slice_to_numeric_range(coord, self.domain[dim_idx])
+        except _util.NonNumericDimensionError:
+            return False  # We only handle numeric dimensions here.
+
+        if not lo_hi:
+            return True
+
+        try:
+            set_dim_range = getattr(sr, f"set_dim_ranges_{dim.type}")
+            set_dim_range(dim.name, [lo_hi])
+            return True
+        except AttributeError:
+            return False
diff --git a/apis/python/src/tiledbsoma/_spatial_util.py b/apis/python/src/tiledbsoma/_spatial_util.py
new file mode 100644
index 0000000000..1cad36ea8a
--- /dev/null
+++ b/apis/python/src/tiledbsoma/_spatial_util.py
@@ -0,0 +1,152 @@
+import json
+from typing import Dict, Optional, Tuple
+
+import numpy as np
+import pyarrow as pa
+import shapely
+import somacore
+from somacore import options
+
+from ._exception import SOMAError
+
+
+def coordinate_space_from_json(data: str) -> somacore.CoordinateSpace:
+    """Returns a coordinate space from a json string."""
+    # TODO: Needs good, comprehensive error handling.
+    raw = json.loads(data)
+    # mypy false positive https://github.com/python/mypy/issues/5313
+    return somacore.CoordinateSpace(tuple(somacore.Axis(**axis) for axis in raw))  # type: ignore[misc]
+
+
+def coordinate_space_to_json(coord_space: somacore.CoordinateSpace) -> str:
+    """Returns json string representation of the coordinate space."""
+    return json.dumps(
+        tuple({"name": axis.name, "unit": axis.unit} for axis in coord_space.axes)
+    )
+
+
+def transform_from_json(data: str) -> somacore.CoordinateTransform:
+    """Creates a coordinate transform from a json string."""
+    raw = json.loads(data)
+    try:
+        transform_type = raw.pop("transform_type")
+    except KeyError as ke:
+        raise KeyError("Missing required key 'transform_type'.") from ke
+    try:
+        kwargs = raw.pop("transform")
+    except KeyError as ke:
+        raise KeyError("Missing reequired key 'transform'.") from ke
+    if transform_type == "IdentityTransform":
+        return somacore.IdentityTransform(**kwargs)
+    elif transform_type == "AffineTransform":
+        return somacore.AffineTransform(**kwargs)
+    else:
+        raise KeyError(f"Unrecognized transform type '{transform_type}'.")
+
+
+def transform_to_json(transform: somacore.CoordinateTransform) -> str:
+    kwargs = {
+        "input_axes": transform.input_axes,
+        "output_axes": transform.output_axes,
+    }
+    if isinstance(transform, somacore.IdentityTransform):
+        pass
+    elif isinstance(transform, somacore.AffineTransform):
+        kwargs["matrix"] = transform.augmented_matrix.tolist()
+    else:
+        raise TypeError(f"Unrecognized coordinate transform type {type(transform)!r}.")
+
+    transform_type = type(transform).__name__
+    return json.dumps({"transform_type": transform_type, "transform": kwargs})
+
+
+def transform_region(
+    region: options.SpatialRegion,
+    transform: somacore.CoordinateTransform,
+) -> shapely.geometry.base.BaseGeometry:
+    if len(transform.input_axes) != 2:
+        raise NotImplementedError(
+            "Spatial queries are currently only supported for 2D coordinates."
+        )
+    if isinstance(region, shapely.geometry.base.BaseGeometry):
+        if region.has_z:
+            raise ValueError("Only 2d shapely geometries are supported.")
+        # Following check is currently unneeded, but leaving it for reference if
+        # 3D support is added.
+        ndim = 3 if region.has_z else 2
+        if ndim != len(transform.input_axes):
+            raise ValueError(
+                f"Input region must have {len(transform.input_axes)} dimension, "
+                f"but region with {ndim} dimensions provided."
+            )
+    else:
+        if len(region) != 4:
+            raise ValueError(f"Unexpected region with size {len(region)}")
+        region = shapely.box(region[0], region[1], region[2], region[3])
+
+    if not isinstance(transform, somacore.AffineTransform):
+        raise NotImplementedError("Only affine transforms are supported.")
+    aug = transform.augmented_matrix
+    affine = aug[:-1, :-1].flatten("C").tolist() + aug[-1, :-1].tolist()
+    return shapely.affinity.affine_transform(region, affine)
+
+
+def process_spatial_df_region(
+    region: Optional[options.SpatialRegion],
+    transform: somacore.CoordinateTransform,
+    coords_by_name: Dict[str, options.SparseDFCoord],
+    index_columns: Tuple[str, ...],
+    axis_names: Tuple[str, ...],
+    schema: pa.Schema,
+) -> Tuple[
+    options.SparseDFCoords,
+    Optional[options.SpatialRegion],
+    somacore.CoordinateTransform,
+]:
+    # Check provided coords are valid.
+    if not set(axis_names).isdisjoint(coords_by_name):
+        raise KeyError("Extra coords cannot contain a spatial index column.")
+    if not set(index_columns).issuperset(coords_by_name):
+        raise KeyError("Extra coords must be index columns.")
+
+    # Transform the region into the data region and add the spatial coordinates
+    # to the coords_by_name map.
+    if region is None:
+        # Leave spatial coords as None - this will select the entire region.
+        data_region: Optional[options.SpatialRegion] = None
+    else:
+        # Restricted to guarantee data region is a box.
+        if isinstance(region, shapely.GeometryType):
+            raise NotImplementedError(
+                "Support for querying point clouds by geometries is not yet implemented."
+            )
+        if not isinstance(transform, somacore.ScaleTransform):
+            raise NotImplementedError(
+                f"Support for querying point clouds with a transform of type "
+                f"{type(transform)!r} our a bounding box region is not yet supported."
+            )
+        # Note: transform_region currently only supports 2D regions. This code block
+        # operates under the same assumption.
+        data_region = transform_region(region, transform)
+
+        def axis_slice(a: float, b: float, dtype: pa.DataType) -> slice:
+            if pa.types.is_floating(dtype):
+                return slice(a, b)
+            if pa.types.is_integer(dtype):
+                # Round so only points within the requested region are returned.
+                return slice(int(np.ceil(a)), int(np.floor(b)))
+            raise SOMAError(f"Unexpected spatial axis datatype {dtype}")
+
+        # Add the transform region to coords. This gets the bounding box for the
+        # requested region.
+        (x_min, y_min, x_max, y_max) = shapely.bounds(data_region)
+        coords_by_name[axis_names[0]] = axis_slice(
+            x_min, x_max, schema.field(axis_names[0]).type
+        )
+        coords_by_name[axis_names[1]] = axis_slice(
+            y_min, y_max, schema.field(axis_names[1]).type
+        )
+
+    coords = tuple(coords_by_name.get(index_name) for index_name in index_columns)
+    inv_transform = transform.inverse_transform()
+    return (coords, data_region, inv_transform)
diff --git a/apis/python/src/tiledbsoma/_tdb_handles.py b/apis/python/src/tiledbsoma/_tdb_handles.py
index decaab86c3..2284636c82 100644
--- a/apis/python/src/tiledbsoma/_tdb_handles.py
+++ b/apis/python/src/tiledbsoma/_tdb_handles.py
@@ -41,6 +41,7 @@
 RawHandle = Union[
     clib.SOMAArray,
     clib.SOMADataFrame,
+    clib.SOMAPointCloudDataFrame,
     clib.SOMASparseNDArray,
     clib.SOMADenseNDArray,
     clib.SOMAGroup,
@@ -77,6 +78,7 @@ def open(
 
     _type_to_class = {
         "somadataframe": DataFrameWrapper,
+        "somapointclouddataframe": PointCloudDataFrameWrapper,
         "somadensendarray": DenseNDArrayWrapper,
         "somasparsendarray": SparseNDArrayWrapper,
         "somacollection": CollectionWrapper,
@@ -513,6 +515,19 @@ def resize_soma_joinid(self, newshape: int) -> None:
         self._handle.resize_soma_joinid(newshape)
 
 
+class PointCloudDataFrameWrapper(SOMAArrayWrapper[clib.SOMAPointCloudDataFrame]):
+    """Wrapper around a Pybind11 SOMADataFrame handle."""
+
+    _ARRAY_WRAPPED_TYPE = clib.SOMAPointCloudDataFrame
+
+    @property
+    def count(self) -> int:
+        return int(self._handle.count)
+
+    def write(self, values: pa.RecordBatch) -> None:
+        self._handle.write(values)
+
+
 class DenseNDArrayWrapper(SOMAArrayWrapper[clib.SOMADenseNDArray]):
     """Wrapper around a Pybind11 DenseNDArrayWrapper handle."""
 
diff --git a/apis/python/src/tiledbsoma/pytiledbsoma.cc b/apis/python/src/tiledbsoma/pytiledbsoma.cc
index f9d56ff7fe..fa97811f04 100644
--- a/apis/python/src/tiledbsoma/pytiledbsoma.cc
+++ b/apis/python/src/tiledbsoma/pytiledbsoma.cc
@@ -20,6 +20,7 @@ void load_soma_context(py::module&);
 void load_soma_object(py::module&);
 void load_soma_array(py::module&);
 void load_soma_dataframe(py::module&);
+void load_soma_point_cloud_dataframe(py::module&);
 void load_soma_dense_ndarray(py::module&);
 void load_soma_sparse_ndarray(py::module&);
 void load_soma_group(py::module&);
@@ -209,6 +210,7 @@ PYBIND11_MODULE(pytiledbsoma, m) {
     load_soma_dataframe(m);
     load_soma_dense_ndarray(m);
     load_soma_sparse_ndarray(m);
+    load_soma_point_cloud_dataframe(m);
     load_soma_group(m);
     load_soma_collection(m);
     load_query_condition(m);
diff --git a/apis/python/src/tiledbsoma/soma_object.cc b/apis/python/src/tiledbsoma/soma_object.cc
index 44388821c1..61561f0756 100644
--- a/apis/python/src/tiledbsoma/soma_object.cc
+++ b/apis/python/src/tiledbsoma/soma_object.cc
@@ -73,6 +73,9 @@ void load_soma_object(py::module& m) {
                     if (soma_obj_type == "somadataframe")
                         return py::cast(
                             dynamic_cast<SOMADataFrame&>(*soma_obj));
+                    if (soma_obj_type == "somapointclouddataframe")
+                        return py::cast(
+                            dynamic_cast<SOMAPointCloudDataFrame&>(*soma_obj));
                     else if (soma_obj_type == "somasparsendarray")
                         return py::cast(
                             dynamic_cast<SOMASparseNDArray&>(*soma_obj));
diff --git a/apis/python/src/tiledbsoma/soma_point_cloud_dataframe.cc b/apis/python/src/tiledbsoma/soma_point_cloud_dataframe.cc
new file mode 100644
index 0000000000..cbc16bb59e
--- /dev/null
+++ b/apis/python/src/tiledbsoma/soma_point_cloud_dataframe.cc
@@ -0,0 +1,153 @@
+/**
+ * @file   soma_point_cloud_dataframe.cc
+ *
+ * @section LICENSE
+ *
+ * The MIT License
+ *
+ * @copyright Copyright (c) 2024 TileDB, Inc.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ * @section DESCRIPTION
+ *
+ * This file defines the SOMAPointCloudDataFrame bindings.
+ */
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/pytypes.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+#include <tiledbsoma/tiledbsoma>
+
+#include "common.h"
+
+namespace libtiledbsomacpp {
+
+namespace py = pybind11;
+using namespace py::literals;
+using namespace tiledbsoma;
+
+void load_soma_point_cloud_dataframe(py::module& m) {
+    py::class_<SOMAPointCloudDataFrame, SOMAArray, SOMAObject>(
+        m, "SOMAPointCloudDataFrame")
+
+        .def_static(
+            "create",
+            [](std::string_view uri,
+               py::object py_schema,
+               py::object index_column_info,
+               std::shared_ptr<SOMAContext> context,
+               PlatformConfig platform_config,
+               std::optional<std::pair<uint64_t, uint64_t>> timestamp) {
+                ArrowSchema schema;
+                uintptr_t schema_ptr = (uintptr_t)(&schema);
+                py_schema.attr("_export_to_c")(schema_ptr);
+
+                // Please see
+                // https://github.com/single-cell-data/TileDB-SOMA/issues/2869
+                // for the reasoning here.
+                //
+                // TL;DR:
+                // * The table has an `ArrowSchema`; each of its children
+                //   is also an `ArrowSchema`.
+                // * Arrow fields are nullable by default in the user API.
+                // * There is a field-level nullability flag, _and_ users
+                //   can set a "nullable" metadata as well.
+                // * In the absence of metadata, respect the flag we get.
+                // * In the present of metdata with "nullable", let that
+                //   override.
+
+                auto metadata = py_schema.attr("metadata");
+                if (py::hasattr(metadata, "get")) {
+                    for (int64_t i = 0; i < schema.n_children; ++i) {
+                        auto child = schema.children[i];
+                        auto val = metadata.attr("get")(
+                            py::str(child->name).attr("encode")("utf-8"));
+
+                        if (!val.is(py::none()) &&
+                            val.cast<std::string>() == "nullable") {
+                            child->flags |= ARROW_FLAG_NULLABLE;
+                        }
+                    }
+                }
+
+                ArrowSchema index_column_schema;
+                ArrowArray index_column_array;
+                uintptr_t
+                    index_column_schema_ptr = (uintptr_t)(&index_column_schema);
+                uintptr_t
+                    index_column_array_ptr = (uintptr_t)(&index_column_array);
+                index_column_info.attr("_export_to_c")(
+                    index_column_array_ptr, index_column_schema_ptr);
+
+                try {
+                    SOMAPointCloudDataFrame::create(
+                        uri,
+                        std::make_unique<ArrowSchema>(schema),
+                        ArrowTable(
+                            std::make_unique<ArrowArray>(index_column_array),
+                            std::make_unique<ArrowSchema>(index_column_schema)),
+                        context,
+                        platform_config,
+                        timestamp);
+                } catch (const std::out_of_range& e) {
+                    throw py::type_error(e.what());
+                } catch (const std::exception& e) {
+                    TPY_ERROR_LOC(e.what());
+                }
+                schema.release(&schema);
+            },
+            "uri"_a,
+            py::kw_only(),
+            "schema"_a,
+            "index_column_info"_a,
+            "ctx"_a,
+            "platform_config"_a,
+            "timestamp"_a = py::none())
+
+        .def_static(
+            "open",
+            py::overload_cast<
+                std::string_view,
+                OpenMode,
+                std::shared_ptr<SOMAContext>,
+                std::vector<std::string>,
+                ResultOrder,
+                std::optional<std::pair<uint64_t, uint64_t>>>(
+                &SOMAPointCloudDataFrame::open),
+            "uri"_a,
+            "mode"_a,
+            "context"_a,
+            py::kw_only(),
+            "column_names"_a = py::tuple(),
+            "result_order"_a = ResultOrder::automatic,
+            "timestamp"_a = py::none())
+
+        .def_static("exists", &SOMAPointCloudDataFrame::exists)
+        .def_property_readonly(
+            "index_column_names", &SOMAPointCloudDataFrame::index_column_names)
+        .def_property_readonly(
+            "count",
+            &SOMAPointCloudDataFrame::count,
+            py::call_guard<py::gil_scoped_release>());
+}
+}  // namespace libtiledbsomacpp


Test patch:
diff --git a/apis/python/tests/test_coordinates.py b/apis/python/tests/test_coordinates.py
new file mode 100644
index 0000000000..046a3e579e
--- /dev/null
+++ b/apis/python/tests/test_coordinates.py
@@ -0,0 +1,63 @@
+import numpy as np
+import pytest
+from somacore import (
+    AffineTransform,
+    Axis,
+    CoordinateSpace,
+    CoordinateTransform,
+    IdentityTransform,
+)
+
+import tiledbsoma as soma
+
+
+def check_transform_is_equal(
+    actual: CoordinateTransform, desired: CoordinateTransform
+) -> None:
+    assert actual.input_axes == desired.input_axes
+    assert actual.output_axes == desired.output_axes
+    if isinstance(desired, IdentityTransform):
+        assert isinstance(actual, IdentityTransform)
+    elif isinstance(desired, AffineTransform):
+        assert isinstance(actual, AffineTransform)
+        np.testing.assert_array_equal(actual.augmented_matrix, desired.augmented_matrix)
+    else:
+        assert False
+
+
+@pytest.mark.parametrize(
+    "original",
+    [
+        CoordinateSpace((Axis(name="dim0", unit="meter"),)),  # type: ignore[arg-type]
+        CoordinateSpace(
+            (
+                Axis(name="dim0", unit="micrometer"),
+                Axis(name="dim1"),
+                Axis(name="dim2", unit="micrometer"),
+            ),  # type: ignore[arg-type]
+        ),
+    ],
+)
+def test_coordinate_system_json_roundtrip(original: CoordinateSpace):
+    json_blob = soma._spatial_util.coordinate_space_to_json(original)
+    result = soma._spatial_util.coordinate_space_from_json(json_blob)
+    assert len(result) == len(original)
+    for index in range(len(result)):
+        assert result[index] == original[index]
+
+
+@pytest.mark.parametrize(
+    "original",
+    [
+        IdentityTransform(["y1", "x1"], ["y2", "x2"]),
+        AffineTransform(
+            ["x1", "y1"],
+            ["x2", "y2"],
+            np.array([[1.5, 3, 0], [-1.5, 3, 1]], dtype=np.float64),
+        ),
+    ],
+)
+def test_transform_json_roundtrip(original: CoordinateTransform):
+    json_blob = soma._spatial_util.transform_to_json(original)
+    result = soma._spatial_util.transform_from_json(json_blob)
+    check_transform_is_equal(result, original)
diff --git a/apis/python/tests/test_point_cloud_dataframe.py b/apis/python/tests/test_point_cloud_dataframe.py
new file mode 100644
index 0000000000..95d50ed99d
--- /dev/null
+++ b/apis/python/tests/test_point_cloud_dataframe.py
@@ -0,0 +1,719 @@
+from urllib.parse import urljoin
+
+import numpy as np
+import pyarrow as pa
+import pytest
+import shapely
+import typeguard
+
+import tiledbsoma as soma
+
+
+def test_point_cloud_bad_create(tmp_path):
+    baseuri = urljoin(f"{tmp_path.as_uri()}/", "bad_create")
+
+    # axis names must be in index column names
+    asch = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+    with pytest.raises(ValueError):
+        soma.PointCloudDataFrame.create(
+            urljoin(baseuri, "bad_name_subset"), schema=asch, index_column_names="x"
+        )
+
+    # all spatial axis must have the same type
+    asch = pa.schema([("x", pa.float64()), ("y", pa.int64())])
+    with pytest.raises(ValueError):
+        soma.PointCloudDataFrame.create(
+            urljoin(baseuri, "different_types"), schema=asch
+        )
+
+    # type must be integral or floating-point
+    asch = pa.schema([("x", pa.large_string()), ("y", pa.large_string())])
+    with pytest.raises(ValueError):
+        soma.PointCloudDataFrame.create(urljoin(baseuri, "bad_type"), schema=asch)
+
+
+def test_point_cloud_basic_read(tmp_path):
+    baseuri = urljoin(f"{tmp_path.as_uri()}/", "basic_read")
+
+    asch = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    # defaults
+    with soma.PointCloudDataFrame.create(
+        urljoin(baseuri, "default"), schema=asch
+    ) as ptc:
+        pydict = {}
+        pydict["soma_joinid"] = [1, 2, 3, 4, 5]
+        pydict["x"] = [10, 20, 30, 40, 50]
+        pydict["y"] = [4.1, 5.2, 6.3, 7.4, 8.5]
+
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(urljoin(baseuri, "default"), "r") as ptc:
+        assert set(ptc.schema.names) == set(ptc.index_column_names)
+        assert ptc.index_column_names == ("soma_joinid", "x", "y")
+        assert ptc.axis_names == ("x", "y")
+
+        table = ptc.read().concat()
+        assert ptc.count == len(ptc) == table.num_rows == 5
+        assert table.num_columns == 3
+        assert [e.as_py() for e in table["soma_joinid"]] == pydict["soma_joinid"]
+        assert [e.as_py() for e in table["x"]] == pydict["x"]
+        assert [e.as_py() for e in table["y"]] == pydict["y"]
+
+    # with user defined values
+    with soma.PointCloudDataFrame.create(
+        urljoin(baseuri, "user_defined"),
+        schema=asch,
+        index_column_names="x",
+        axis_names="x",
+        domain=((1, 10),),
+    ) as ptc:
+        pydict = {}
+        pydict["soma_joinid"] = [1, 2, 3, 4, 5]
+        pydict["x"] = [1, 2, 3, 4, 5]
+        pydict["y"] = [4.1, 5.2, 6.3, 7.4, 8.5]
+
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(urljoin(baseuri, "user_defined"), "r") as ptc:
+        assert set(ptc.schema.names) == set(["soma_joinid", "x", "y"])
+        assert ptc.index_column_names == ("x",)
+        assert ptc.axis_names == ("x",)
+        assert ptc.domain == ((1, 10),)
+
+        table = ptc.read().concat()
+        assert ptc.count == len(ptc) == table.num_rows == 5
+        assert table.num_columns == 3
+        assert [e.as_py() for e in table["soma_joinid"]] == pydict["soma_joinid"]
+        assert [e.as_py() for e in table["x"]] == pydict["x"]
+        assert [e.as_py() for e in table["y"]] == pydict["y"]
+
+
+def test_point_cloud_coordinate_space(tmp_path):
+    uri = tmp_path.as_uri()
+
+    asch = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=asch) as ptc:
+        assert len(ptc.coordinate_space) == 2
+        assert ptc.coordinate_space.axis_names == ("x", "y")
+        assert ptc.coordinate_space.axes == (soma.Axis(name="x"), soma.Axis(name="y"))
+
+        # Axis names do not match
+        with pytest.raises(ValueError):
+            ptc.coordinate_space = soma.CoordinateSpace(
+                [soma.Axis(name="a"), soma.Axis(name="y")]
+            )
+
+        ptc.coordinate_space = soma.CoordinateSpace(
+            [soma.Axis(name="x", unit="m"), soma.Axis(name="y", unit="in")]
+        )
+        assert ptc.coordinate_space[0] == soma.Axis(name="x", unit="m")
+        assert ptc.coordinate_space[1] == soma.Axis(name="y", unit="in")
+
+
+def test_point_cloud_bad_read_spatial_region(tmp_path):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        # Cannot specify the output coordinate space when transform is None
+        with pytest.raises(ValueError):
+            ptc.read_spatial_region(
+                region_transform=None,
+                region_coord_space=soma.CoordinateSpace(
+                    (soma.Axis(name="x"), soma.Axis(name="y"))
+                ),
+            )
+
+        # The input axes of the transorm must match the axes of the coordinate
+        # space the requested region is defined in
+        with pytest.raises(ValueError):
+            ptc.read_spatial_region(
+                region_transform=soma.IdentityTransform(
+                    input_axes=("x", "y"), output_axes=("x", "y")
+                ),
+                region_coord_space=soma.CoordinateSpace(
+                    [soma.Axis(name="y"), soma.Axis(name="x")]
+                ),
+            )
+
+
+@pytest.mark.parametrize(
+    "name, region, value_filter, expected_output",
+    [
+        (
+            "Bounding box",
+            [10, 4.1, 50, 8.5],
+            None,
+            {
+                "soma_joinid": [1, 2, 3, 4, 5],
+                "x": [10, 20, 30, 40, 50],
+                "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+            },
+        ),
+        (
+            "Full region",
+            [0, 0, 100, 100],
+            None,
+            {
+                "soma_joinid": [1, 2, 3, 4, 5],
+                "x": [10, 20, 30, 40, 50],
+                "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+            },
+        ),
+        (
+            "Partial region",
+            [15, 4, 35, 7],
+            None,
+            {
+                "soma_joinid": [2, 3],
+                "x": [20, 30],
+                "y": [5.2, 6.3],
+            },
+        ),
+        (
+            "Single point",
+            [10, 4.1, 10, 4.1],
+            None,
+            {
+                "soma_joinid": [1],
+                "x": [10],
+                "y": [4.1],
+            },
+        ),
+        (
+            "No points",
+            [60, 10, 80, 20],
+            None,
+            {
+                "soma_joinid": [],
+                "x": [],
+                "y": [],
+            },
+        ),
+        (
+            "Flipped box",
+            [50, 10, 20, 0],
+            None,
+            {
+                "soma_joinid": [2, 3, 4, 5],
+                "x": [20, 30, 40, 50],
+                "y": [5.2, 6.3, 7.4, 8.5],
+            },
+        ),
+        (
+            "Bounding box with value filter",
+            [10, 4.1, 50, 8.5],
+            "soma_joinid in [1, 2, 3, 4, 5]",
+            {
+                "soma_joinid": [1, 2, 3, 4, 5],
+                "x": [10, 20, 30, 40, 50],
+                "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+            },
+        ),
+        (
+            "Partial region with value filter",
+            [15, 4, 35, 7],
+            "soma_joinid in [2, 3]",
+            {
+                "soma_joinid": [2, 3],
+                "x": [20, 30],
+                "y": [5.2, 6.3],
+            },
+        ),
+        (
+            "Single point with value filter",
+            [10, 4.1, 10, 4.1],
+            "soma_joinid == 1",
+            {
+                "soma_joinid": [1],
+                "x": [10],
+                "y": [4.1],
+            },
+        ),
+        (
+            "No points with value filter",
+            [60, 10, 80, 20],
+            "soma_joinid in [1, 2, 3, 4, 5]",
+            {
+                "soma_joinid": [],
+                "x": [],
+                "y": [],
+            },
+        ),
+        (
+            "Flipped box with value filter",
+            [50, 10, 20, 0],
+            "soma_joinid in [2, 3, 4, 5]",
+            {
+                "soma_joinid": [2, 3, 4, 5],
+                "x": [20, 30, 40, 50],
+                "y": [5.2, 6.3, 7.4, 8.5],
+            },
+        ),
+    ],
+)
+def test_point_cloud_read_spatial_region_basic_2d(
+    tmp_path, name, region, value_filter, expected_output
+):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        actual_output = ptc.read_spatial_region(
+            region=region, value_filter=value_filter
+        )
+        assert actual_output.data.concat().to_pydict() == expected_output
+
+        actual_output = ptc.read_spatial_region(
+            region=shapely.box(*region), value_filter=value_filter
+        )
+        assert actual_output.data.concat().to_pydict() == expected_output
+
+        actual_output = ptc.read_spatial_region(
+            region=region, column_names=["x"], value_filter=value_filter
+        )
+        if value_filter:
+            assert actual_output.data.concat().to_pydict() == {
+                "soma_joinid": expected_output["soma_joinid"],
+                "x": expected_output["x"],
+            }
+        else:
+            assert actual_output.data.concat().to_pydict() == {
+                "x": expected_output["x"],
+            }
+
+        actual_output = ptc.read_spatial_region(
+            region=region, column_names=["y"], value_filter=value_filter
+        )
+        if value_filter:
+            assert actual_output.data.concat().to_pydict() == {
+                "soma_joinid": expected_output["soma_joinid"],
+                "y": expected_output["y"],
+            }
+        else:
+            assert actual_output.data.concat().to_pydict() == {
+                "y": expected_output["y"],
+            }
+
+
+@pytest.mark.skip("3D regions not supported yet")
+@pytest.mark.parametrize(
+    "name, region, expected_output",
+    [
+        (
+            "Bounding box",
+            [10, 4.1, 50, 8.5, 0, 100],
+            {
+                "soma_joinid": [1, 2, 3, 4, 5],
+                "x": [10, 20, 30, 40, 50],
+                "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+                "z": [0, 1, 2, 3, 4],
+            },
+        ),
+        (
+            "Full region",
+            [0, 0, 100, 100, 0, 10],
+            {
+                "soma_joinid": [1, 2, 3, 4, 5],
+                "x": [10, 20, 30, 40, 50],
+                "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+                "z": [0, 1, 2, 3, 4],
+            },
+        ),
+        (
+            "Partial region",
+            [15, 4, 35, 7, 0, 5],
+            {
+                "soma_joinid": [2, 3],
+                "x": [20, 30],
+                "y": [5.2, 6.3],
+                "z": [1, 2],
+            },
+        ),
+        (
+            "Single point",
+            [10, 4.1, 10, 4.1, 0, 0],
+            {
+                "soma_joinid": [1],
+                "x": [10],
+                "y": [4.1],
+                "z": [0],
+            },
+        ),
+        (
+            "No points",
+            [60, 10, 80, 20, 0, 1],
+            {
+                "soma_joinid": [],
+                "x": [],
+                "y": [],
+                "z": [],
+            },
+        ),
+        (
+            "Flipped box",
+            [50, 10, 20, 0, 0, 5],
+            {
+                "soma_joinid": [2, 3, 4, 5],
+                "x": [20, 30, 40, 50],
+                "y": [5.2, 6.3, 7.4, 8.5],
+                "z": [1, 2, 3, 4],
+            },
+        ),
+    ],
+)
+def test_point_cloud_read_spatial_region_basic_3d(
+    tmp_path, name, region, expected_output
+):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64()), ("z", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(
+        uri,
+        schema=schema,
+        index_column_names=("soma_joinid", "x", "y", "z"),
+        axis_names=("x", "y", "z"),
+    ) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+            "z": [0, 1, 2, 3, 4],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        actual_output = ptc.read_spatial_region(region=region)
+        assert actual_output.data.concat().to_pydict() == expected_output
+
+
+@pytest.mark.parametrize(
+    "name, region, exc_type",
+    [
+        ("Empty", [], ValueError),
+        ("Too few", [10, 4.1], ValueError),
+        ("Too many", [0, 0, 100, 100, 50], ValueError),
+        ("Bad type", ["a", "b", "c", "d"], typeguard.TypeCheckError),
+    ],
+)
+def test_point_cloud_read_spatial_region_2d_bad(tmp_path, name, region, exc_type):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        with pytest.raises(exc_type):
+            ptc.read_spatial_region(region=region)
+
+
+@pytest.mark.skip("3D regions not supported yet")
+@pytest.mark.parametrize(
+    "name, region, exc_type",
+    [
+        ("Empty", [], ValueError),
+        ("Too few", [10, 4.1, 20, 6.3], ValueError),
+        ("Too many", [0, 0, 100, 100, 50, 10, 0], ValueError),
+        ("Bad type", ["a", "b", "c", "d", "e", "f"], typeguard.TypeCheckError),
+    ],
+)
+def test_point_cloud_read_spatial_region_3d_bad(tmp_path, name, region, exc_type):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64()), ("z", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+            "z": [0, 1, 2, 3, 4],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        with pytest.raises(exc_type):
+            ptc.read_spatial_region(region=region)
+
+
+def point_cloud_read_spatial_region_transform_setup(uri, transform, input_axes, kwargs):
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    output_names = ("x", "y")
+    input_names = tuple(axis.name for axis in input_axes)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        read_spatial_region = ptc.read_spatial_region(
+            region_transform=transform(
+                input_axes=input_names, output_axes=output_names, **kwargs
+            )
+        )
+
+        assert read_spatial_region.data.concat() == ptc.read().concat()
+
+        data_coordinate_space = read_spatial_region.data_coordinate_space
+        assert (
+            data_coordinate_space.axis_names
+            == ptc.coordinate_space.axis_names
+            == output_names
+        )
+        assert len(data_coordinate_space) == len(output_names)
+
+        output_coordinate_space = read_spatial_region.output_coordinate_space
+        assert output_coordinate_space.axes == input_axes
+        assert output_coordinate_space.axis_names == input_names
+        assert len(output_coordinate_space) == len(input_names)
+
+        assert isinstance(read_spatial_region.coordinate_transform, transform)
+
+    return read_spatial_region.coordinate_transform
+
+
+def test_point_cloud_read_spatial_region_identity_transform(tmp_path):
+    coordinate_transform = point_cloud_read_spatial_region_transform_setup(
+        uri=tmp_path.as_uri(),
+        transform=soma.IdentityTransform,
+        input_axes=(soma.Axis(name="x"), soma.Axis(name="y")),
+        kwargs={},
+    )
+
+    # Identity will always have a scale of 1.0
+    assert list(coordinate_transform.scale_factors) == [1, 1]
+    assert coordinate_transform.scale == 1.0
+    assert list(coordinate_transform.inverse_transform().scale_factors) == [1, 1]
+
+
+@pytest.mark.parametrize(
+    "name, scale, expected_scale_factors, expected_scale",
+    [
+        ("Identity", 1, [1, 1], 1.0),
+        ("Scale up", 2, [0.5, 0.5], 0.5),
+        ("Negative scale up", -2, [-0.5, -0.5], -0.5),
+        ("Very large", 1e6, [1e-6, 1e-6], 1e-6),
+        ("Very small", 1e-6, [1e6, 1e6], 1e6),
+        ("Scale down", 0.5, [2, 2], 2.0),
+        ("Negative scale down", -0.5, [-2, -2], -2.0),
+        ("Inversion", -1, [-1, -1], -1.0),
+        ("Fractional scale", 0.75, [4 / 3, 4 / 3], 4 / 3),
+    ],
+)
+def test_point_cloud_read_spatial_region_uniform_scale_transform(
+    tmp_path, name, scale, expected_scale_factors, expected_scale
+):
+    coordinate_transform = point_cloud_read_spatial_region_transform_setup(
+        uri=tmp_path.as_uri(),
+        transform=soma.UniformScaleTransform,
+        input_axes=(soma.Axis(name="x2"), soma.Axis(name="y2")),
+        kwargs={"scale": scale},
+    )
+
+    assert list(coordinate_transform.scale_factors) == expected_scale_factors
+    assert coordinate_transform.scale == expected_scale
+    assert coordinate_transform.inverse_transform().scale == scale
+
+
+def test_point_cloud_read_spatial_region_uniform_scale_transform_bad(tmp_path):
+    # Error out for zero scale
+    with pytest.raises(ZeroDivisionError):
+        point_cloud_read_spatial_region_transform_setup(
+            uri=tmp_path.as_uri(),
+            transform=soma.UniformScaleTransform,
+            input_axes=(soma.Axis(name="x"), soma.Axis(name="y")),
+            kwargs={"scale": 0},
+        )
+
+
+@pytest.mark.parametrize(
+    "name, scale_factors, expected_scale_factors",
+    [
+        ["Identity", (1, 1), (1, 1)],
+        ["Scale up and down", (2, 0.5), (0.5, 2.0)],
+        ["Very large and small scale", (1e6, 1e-6), (1e-6, 1e6)],
+        ["Single zero scale", (0, 2), (float("inf"), 0.5)],
+        ["Double zero scale", (0, 0), (float("inf"), float("inf"))],
+        ["Single inversion", (-1, 2), (-1, 0.5)],
+        ["Double inversion", (-2, -0.5), (-0.5, -2)],
+        ["One dimension inverted", (1, -1), (1, -1)],
+        ["Scale up in both dimensions", (2, 3), (0.5, (1 / 3))],
+        ["Small scale factors", (0.1, 0.1), (10, 10)],
+    ],
+)
+def test_point_cloud_read_spatial_region_scale_transform(
+    tmp_path, name, scale_factors, expected_scale_factors
+):
+    coordinate_transform = point_cloud_read_spatial_region_transform_setup(
+        uri=tmp_path.as_uri(),
+        transform=soma.ScaleTransform,
+        input_axes=(soma.Axis(name="x2"), soma.Axis(name="y2")),
+        kwargs={"scale_factors": scale_factors},
+    )
+
+    assert tuple(coordinate_transform.scale_factors) == expected_scale_factors
+    assert (
+        tuple(coordinate_transform.inverse_transform().scale_factors) == scale_factors
+    )
+
+
+@pytest.mark.parametrize(
+    "name, matrix, expected_aug_matrix",
+    [
+        (
+            "Identity",
+            np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+            np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+        ),
+        (
+            "Translation",
+            np.array([[1, 0, 2], [0, 1, 3], [0, 0, 1]]),
+            np.array([[1, 0, -2], [0, 1, -3], [0, 0, 1]]),
+        ),
+        (
+            "Scaling",
+            np.array([[2, 0, 0], [0, 3, 0], [0, 0, 1]]),
+            np.array([[0.5, 0, 0], [0, 1 / 3, 0], [0, 0, 1]]),
+        ),
+        (
+            "Rotation (90 degrees counter-clockwise)",
+            np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]),
+            np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]]),
+        ),
+        (
+            "Reflection over the x-axis",
+            np.array([[1, 0, 0], [0, -1, 0], [0, 0, 1]]),
+            np.array([[1, 0, 0], [0, -1, 0], [0, 0, 1]]),
+        ),
+        (
+            "Reflection over the y-axis",
+            np.array([[-1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+            np.array([[-1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+        ),
+        (
+            "Shearing (x direction)",
+            np.array([[1, 1, 0], [0, 1, 0], [0, 0, 1]]),
+            np.array([[1, -1, 0], [0, 1, 0], [0, 0, 1]]),
+        ),
+        (
+            "Shearing (y direction)",
+            np.array([[1, 0, 0], [1, 1, 0], [0, 0, 1]]),
+            np.array([[1, 0, 0], [-1, 1, 0], [0, 0, 1]]),
+        ),
+        (
+            "Scaling and translation",
+            np.array([[2, 0, 1], [0, 3, 2], [0, 0, 1]]),
+            np.array([[0.5, 0, -0.5], [0, 1 / 3, -2 / 3], [0, 0, 1]]),
+        ),
+        (
+            "Rotation and scaling",
+            np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]),
+            np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]]),
+        ),
+    ],
+)
+def test_point_cloud_read_spatial_region_affine_transform(
+    tmp_path, name, matrix, expected_aug_matrix
+):
+    coordinate_transform = point_cloud_read_spatial_region_transform_setup(
+        uri=tmp_path.as_uri(),
+        transform=soma.AffineTransform,
+        input_axes=(soma.Axis(name="x2"), soma.Axis(name="y2")),
+        kwargs={"matrix": matrix},
+    )
+
+    aug_matrix = coordinate_transform.augmented_matrix
+    assert np.array_equal(aug_matrix, expected_aug_matrix)
+
+    aug_matrix_inv = coordinate_transform.inverse_transform().augmented_matrix
+    assert np.array_equal(aug_matrix_inv, matrix)
+
+
+def test_point_cloud_read_spatial_region_region_coord_space(tmp_path):
+    uri = tmp_path.as_uri()
+
+    schema = pa.schema([("x", pa.float64()), ("y", pa.float64())])
+
+    with soma.PointCloudDataFrame.create(uri, schema=schema) as ptc:
+        pydict = {
+            "soma_joinid": [1, 2, 3, 4, 5],
+            "x": [10, 20, 30, 40, 50],
+            "y": [4.1, 5.2, 6.3, 7.4, 8.5],
+        }
+        rb = pa.Table.from_pydict(pydict)
+        ptc.write(rb)
+
+    with soma.PointCloudDataFrame.open(uri, "r") as ptc:
+        output = ptc.read_spatial_region()
+        assert output.output_coordinate_space.axis_names == ("x", "y")
+
+        output = ptc.read_spatial_region(
+            region_transform=soma.IdentityTransform(
+                input_axes=("a", "b"), output_axes=("x", "y")
+            ),
+            region_coord_space=soma.CoordinateSpace(
+                [soma.Axis(name="a"), soma.Axis(name="b")]
+            ),
+        )
+        assert output.output_coordinate_space.axis_names == ("a", "b")
+
+        # Cannot specify the output coordinate space when transform is None
+        with pytest.raises(ValueError):
+            ptc.read_spatial_region(
+                region_coord_space=soma.CoordinateSpace(
+                    [soma.Axis(name="x"), soma.Axis(name="y")]
+                )
+            )
+
+        # The input axes of the transform must match the axes of the coordinate
+        # space the requested region is defined in
+        with pytest.raises(ValueError):
+            ptc.read_spatial_region(
+                region_transform=soma.IdentityTransform(
+                    input_axes=("a", "b"), output_axes=("x", "y")
+                ),
+                region_coord_space=soma.CoordinateSpace(
+                    [soma.Axis(name="x"), soma.Axis(name="y")]
+                ),
+            )