diff --git a/apis/python/src/tiledbsoma/__init__.py b/apis/python/src/tiledbsoma/__init__.py
index 569a785e44..cc344ce5fe 100644
--- a/apis/python/src/tiledbsoma/__init__.py
+++ b/apis/python/src/tiledbsoma/__init__.py
@@ -175,6 +175,7 @@
 )
 from ._indexer import IntIndexer, tiledbsoma_build_index
 from ._measurement import Measurement
+from ._multiscale_image import MultiscaleImage
 from ._point_cloud_dataframe import PointCloudDataFrame
 from ._sparse_nd_array import SparseNDArray, SparseNDArrayRead
 from .options import SOMATileDBContext, TileDBCreateOptions, TileDBWriteOptions
@@ -209,6 +210,7 @@
     "get_storage_engine",
     "IntIndexer",
     "Measurement",
+    "MultiscaleImage",
     "NotCreateableError",
     "open",
     "PointCloudDataFrame",
diff --git a/apis/python/src/tiledbsoma/_arrow_types.py b/apis/python/src/tiledbsoma/_arrow_types.py
index 96ef3b68d8..0d6a580933 100644
--- a/apis/python/src/tiledbsoma/_arrow_types.py
+++ b/apis/python/src/tiledbsoma/_arrow_types.py
@@ -74,6 +74,23 @@
     pa.timestamp("ns"): "tsn:",
 }
 
+_CARROW_TO_PYARROW: Dict[pa.DataType, str] = {
+    "c": pa.int8(),
+    "s": pa.int16(),
+    "i": pa.int32(),
+    "l": pa.int64(),
+    "C": pa.uint8(),
+    "S": pa.uint16(),
+    "I": pa.uint32(),
+    "L": pa.uint64(),
+    "f": pa.float32(),
+    "g": pa.float64(),
+    "tss:": pa.timestamp("s"),
+    "tsm:": pa.timestamp("ms"),
+    "tsu:": pa.timestamp("us"),
+    "tsn:": pa.timestamp("ns"),
+}
+
 # Same as _ARROW_TO_TDB_ATTR, but used for DataFrame indexed columns, aka TileDB Dimensions.
 # Any type system differences from the base-case Attr should be added here.
 _ARROW_TO_TDB_DIM: Dict[Any, Union[str, TypeError]] = _ARROW_TO_TDB_ATTR.copy()
@@ -284,3 +301,10 @@ def pyarrow_to_carrow_type(pa_type: pa.DataType) -> str:
         return _PYARROW_TO_CARROW[pa_type]
     except KeyError:
         raise TypeError(f"Invalid pyarrow type {pa_type}") from None
+
+
+def carrow_type_to_pyarrow(ca_type: str) -> pa.DataType:
+    try:
+        return _CARROW_TO_PYARROW[ca_type]
+    except KeyError:
+        raise TypeError(f"Invalid carrrow type {ca_type}") from None
diff --git a/apis/python/src/tiledbsoma/_constants.py b/apis/python/src/tiledbsoma/_constants.py
index 0e29438e17..3e139ea728 100644
--- a/apis/python/src/tiledbsoma/_constants.py
+++ b/apis/python/src/tiledbsoma/_constants.py
@@ -9,7 +9,7 @@
 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_MULTISCALE_IMAGE_SCHEMA = "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 a04befdc00..7bdb4044f3 100644
--- a/apis/python/src/tiledbsoma/_factory.py
+++ b/apis/python/src/tiledbsoma/_factory.py
@@ -28,6 +28,7 @@
     _dense_nd_array,
     _experiment,
     _measurement,
+    _multiscale_image,
     _point_cloud_dataframe,
     _soma_object,
     _sparse_nd_array,
@@ -216,6 +217,7 @@ def _type_name_to_cls(type_name: str) -> Type[AnySOMAObject]:
             _dense_nd_array.DenseNDArray,
             _experiment.Experiment,
             _measurement.Measurement,
+            _multiscale_image.MultiscaleImage,
             _sparse_nd_array.SparseNDArray,
             _point_cloud_dataframe.PointCloudDataFrame,
         )
diff --git a/apis/python/src/tiledbsoma/_multiscale_image.py b/apis/python/src/tiledbsoma/_multiscale_image.py
index 6bd30edd62..2b7e499fde 100644
--- a/apis/python/src/tiledbsoma/_multiscale_image.py
+++ b/apis/python/src/tiledbsoma/_multiscale_image.py
@@ -7,20 +7,45 @@
 Implementation of a SOMA MultiscaleImage.
 """
 
+import json
 import warnings
 from typing import Any, Optional, Sequence, Tuple, Union
 
+import attrs
 import pyarrow as pa
 import somacore
-from somacore import CoordinateSpace, CoordinateTransform, ScaleTransform, options
-from typing_extensions import Self
-
-from ._constants import SPATIAL_DISCLAIMER
+from somacore import (
+    Axis,
+    CoordinateSpace,
+    CoordinateTransform,
+    ScaleTransform,
+    options,
+)
+from typing_extensions import Final, Self
+
+from . import _funcs, _tdb_handles
+from . import pytiledbsoma as clib
+from ._arrow_types import carrow_type_to_pyarrow, pyarrow_to_carrow_type
+from ._constants import (
+    SOMA_COORDINATE_SPACE_METADATA_KEY,
+    SOMA_MULTISCALE_IMAGE_SCHEMA,
+    SPATIAL_DISCLAIMER,
+)
 from ._dense_nd_array import DenseNDArray
+from ._exception import SOMAError, map_exception_for_create
+from ._soma_group import SOMAGroup
 from ._soma_object import AnySOMAObject
+from ._spatial_util import (
+    coordinate_space_from_json,
+    coordinate_space_to_json,
+    process_image_region,
+)
+from ._types import OpenTimestamp
 from .options import SOMATileDBContext
+from .options._soma_tiledb_context import _validate_soma_tiledb_context
 
 
+@attrs.define(frozen=True)
 class ImageProperties:
     """Properties for a single resolution level in a multiscale image.
 
@@ -28,26 +53,57 @@ class ImageProperties:
         Experimental.
     """
 
-    @property
-    def name(self) -> str:
-        """The key for the image.
-
-        Lifecycle:
-            Experimental
-        """
-        raise NotImplementedError()
-
-    @property
-    def shape(self) -> Tuple[int, ...]:
-        """Size of each axis of the image.
-
-        Lifecycle:
-            Experimental.
-        """
-        raise NotImplementedError()
-
-
-class MultiscaleImage(somacore.MultiscaleImage[DenseNDArray, AnySOMAObject]):
+    name: str
+    image_type: str
+    shape: Tuple[int, ...] = attrs.field(converter=tuple)
+    width: int = attrs.field(init=False)
+    height: int = attrs.field(init=False)
+    depth: Optional[int] = attrs.field(init=False)
+    nchannels: Optional[int] = attrs.field(init=False)
+
+    def __attrs_post_init__(self):  # type: ignore[no-untyped-def]
+        if len(self.image_type) != len(set(self.image_type)):
+            raise ValueError(
+                f"Invalid image type '{self.image_type}'. Image type cannot contain "
+                f"repeated values."
+            )
+        if len(self.image_type) != len(self.shape):
+            raise ValueError(
+                f"{len(self.image_type)} axis names must be provided for a multiscale "
+                f"image with image type {self.image_type}."
+            )
+
+        nchannels: Optional[int] = None
+        width: Optional[int] = None
+        height: Optional[int] = None
+        depth: Optional[int] = None
+        for val, size in zip(self.image_type, self.shape):
+            if val == "X":
+                width = size
+            elif val == "Y":
+                height = size
+            elif val == "Z":
+                depth = size
+            elif val == "C":
+                nchannels = size
+            else:
+                raise SOMAError(f"Invalid image type '{self.image_type}'")
+        if width is None or height is None:
+            raise ValueError(
+                f"Invalid image type '{self.image_type}'. Image type must include "
+                f"'X' and 'Y'."
+            )
+
+        object.__setattr__(self, "nchannels", nchannels)
+        object.__setattr__(self, "width", width)
+        object.__setattr__(self, "height", height)
+        object.__setattr__(self, "depth", depth)
+
+
+class MultiscaleImage(  # type: ignore[misc]  # __eq__ false positive
+    SOMAGroup[DenseNDArray],
+    somacore.MultiscaleImage[DenseNDArray, AnySOMAObject],
+):
     """A multiscale image represented as a collection of images at multiple resolution levels.
 
     Each level of the multiscale image must have the following consistent properties:
@@ -59,7 +115,10 @@ class MultiscaleImage(somacore.MultiscaleImage[DenseNDArray, AnySOMAObject]):
         Experimental.
     """
 
-    __slots__ = ()
+    __slots__ = ("_schema", "_coord_space", "_levels")
+    _wrapper_type = _tdb_handles.MultiscaleImageWrapper
+
+    _level_prefix: Final = "soma_level_"
 
     # Lifecycle
 
@@ -74,6 +133,7 @@ def create(
         axis_types: Sequence[str] = ("channel", "height", "width"),
         platform_config: Optional[options.PlatformConfig] = None,
         context: Optional[SOMATileDBContext] = None,
+        tiledb_timestamp: Optional[OpenTimestamp] = None,
     ) -> Self:
         """Creates a new ``MultiscaleImage`` at the given URI.
 
@@ -93,9 +153,109 @@ def create(
         Lifecycle:
             Experimental.
         """
+        # Warn about the experimental nature of the spatial classes.
         warnings.warn(SPATIAL_DISCLAIMER)
-        raise NotImplementedError()
 
+        context = _validate_soma_tiledb_context(context)
+        if len(set(axis_types)) != len(axis_types):
+            raise ValueError(
+                "Invalid axis types {axis_types} - cannot have repeated values."
+            )
+        if len(axis_names) != len(axis_types):
+            raise ValueError("Mismatched lengths for axis names and types.")
+        axis_type_map = {"channel": "C", "height": "Y", "width": "X", "depth": "Z"}
+        image_type = []
+        for val in axis_types:
+            try:
+                image_type.append(axis_type_map[val])
+            except KeyError as ke:
+                raise ValueError("Invalid axis type name '{val}'.") from ke
+        schema = MultiscaleImageSchema(
+            ImageProperties(
+                name="reference_level",
+                image_type="".join(image_type),
+                shape=tuple(reference_level_shape),  # type: ignore
+            ),
+            axis_names=tuple(axis_names),
+            datatype=type,
+        )
+
+        # mypy false positive https://github.com/python/mypy/issues/5313
+        coord_space = CoordinateSpace(
+            tuple(Axis(name) for name in schema.get_coordinate_space_axis_names())  # type: ignore[misc]
+        )
+        schema_str = schema.to_json()
+        coord_space_str = coordinate_space_to_json(coord_space)
+        try:
+            timestamp_ms = context._open_timestamp_ms(tiledb_timestamp)
+            clib.SOMAGroup.create(
+                uri=uri,
+                soma_type=somacore.MultiscaleImage.soma_type,
+                ctx=context.native_context,
+                timestamp=(0, timestamp_ms),
+            )
+            handle = _tdb_handles.MultiscaleImageWrapper.open(
+                uri, "w", context, tiledb_timestamp
+            )
+            handle.metadata[SOMA_MULTISCALE_IMAGE_SCHEMA] = schema_str
+            handle.metadata[SOMA_COORDINATE_SPACE_METADATA_KEY] = coord_space_str
+            return cls(
+                handle,
+                _dont_call_this_use_create_or_open_instead="tiledbsoma-internal-code",
+            )
+        except SOMAError as e:
+            raise map_exception_for_create(e, uri) from None
+
+    def __init__(
+        self,
+        handle: _tdb_handles.SOMAGroupWrapper[Any],
+        **kwargs: Any,
+    ):
+        # Do generic SOMA collection initialization.
+        super().__init__(handle, **kwargs)
+
+        # Get schema for the multiscale image.
+        try:
+            schema_json = self.metadata[SOMA_MULTISCALE_IMAGE_SCHEMA]
+        except KeyError as ke:
+            raise SOMAError("Missing multiscale image schema metadata") from ke
+        if isinstance(schema_json, bytes):
+            schema_json = str(schema_json, "utf-8")
+        if not isinstance(schema_json, str):
+            raise SOMAError(
+                f"Stored '{SOMA_MULTISCALE_IMAGE_SCHEMA}' metadata is unexpected "
+                f"type {type(schema_json)!r}."
+            )
+        self._schema = MultiscaleImageSchema.from_json(schema_json)
+
+        # Get the 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)
+
+        # Check schema and coordinate space have the same axis order
+        schema_axes = self._schema.get_coordinate_space_axis_names()
+        if schema_axes != self._coord_space.axis_names:
+            raise SOMAError(
+                f"Inconsistent axis names stored in metadata. Multiscale schema metadata"
+                f" has coordinate axes '{schema_axes}', but the coordinate space "
+                f"metadata has coordinate axes '{self._coord_space.axis_names}'"
+            )
+
+        # Get the image levels.
+        # TODO: Optimize and push down to C++ level
+        self._levels = [
+            ImageProperties(name=key[len(self._level_prefix) :], **json.loads(val))
+            for key, val in self.metadata.items()
+            if key.startswith(self._level_prefix)
+        ]
+        self._levels.sort(key=lambda level: (-level.width, -level.height, level.name))
+
+    @_funcs.forwards_kwargs_to(
+        DenseNDArray.create, exclude=("context", "shape", "tiledb_timestamp")
+    )
     def add_new_level(
         self,
         key: str,
@@ -113,7 +273,71 @@ def add_new_level(
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        # Check if key already exists in either the collection or level metadata.
+        if key in self:
+            raise KeyError(f"{key!r} already exists in {type(self)}")
+        meta_key = f"{self._level_prefix}{key}"
+        if meta_key in self.metadata:
+            raise KeyError(f"{key!r} already exists in {type(self)} scales")
+
+        # Check if the shape is valid.
+        ref_props = self._schema.reference_level_properties
+        shape = tuple(shape)
+        if len(shape) != len(ref_props.shape):
+            raise ValueError(
+                f"New level must have {len(ref_props.shape)} dimensions, but shape "
+                f"{shape} has {len(shape)} dimensions."
+            )
+
+        # Check, create, and store as metadata the new level image properties.
+        props = ImageProperties(
+            image_type=ref_props.image_type,
+            name=key,
+            shape=shape,  # type: ignore
+        )
+        if ref_props.nchannels is not None and ref_props.nchannels != props.nchannels:
+            raise ValueError(
+                f"New level must have {ref_props.nchannels}, but provided shape has "
+                f"{props.nchannels} channels."
+            )
+
+        props_str = json.dumps(
+            {
+                "image_type": ref_props.image_type,
+                "shape": shape,
+            }
+        )
+        self.metadata[meta_key] = props_str
+
+        # Add the level properties to level list.
+        # Note: The names are guaranteed to be different from the earlier checks.
+        for index, val in enumerate(self._levels):
+            # Note: Name is unique, so guaranteed to be strict ordering.
+            if (-props.width, -props.height, props.name) < (
+                -val.width,
+                -val.height,
+                val.name,
+            ):
+                self._levels.insert(index, props)
+                break
+        else:
+            self._levels.append(props)
+
+        # Create and return new level array.
+
+        return self._add_new_element(
+            key,
+            DenseNDArray,
+            lambda create_uri: DenseNDArray.create(
+                create_uri,
+                context=self.context,
+                tiledb_timestamp=self.tiledb_timestamp_ms,
+                shape=props.shape,
+                type=self._schema.datatype,
+                **kwargs,
+            ),
+            uri,
+        )
 
     # Data operations
 
@@ -158,7 +382,91 @@ def read_spatial_region(
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        # Get reference level. Check image is 2D.
+        if self._schema.reference_level_properties.depth is not None:
+            raise NotImplementedError(
+                "Support for reading the levels of 3D images it not yet implemented."
+            )
+
+        # Check input query region type is supported.
+        if (
+            channel_coords is not None
+            and self._schema.reference_level_properties.nchannels is None
+        ):
+            raise ValueError(
+                "Invalide channel coordinate provided. This image has no channel "
+                "dimension."
+            )
+
+        # Get the transformation for the group and the data coordinate space.
+        # We may want to revisit copying the units for the data coordinate space.
+        group_to_level = self.get_transform_to_level(level)
+        data_coord_space = self.coordinate_space
+
+        # Update transform and set region coordinate space.
+        # - Add transformation from reference coord system to requested image level.
+        # - Create or check the coordinate space for the input data region.
+        if region_transform is None:
+            if region_coord_space is not None:
+                raise ValueError(
+                    "Cannot specify the output coordinate space when region transform "
+                    "is ``None``."
+                )
+            region_transform = group_to_level
+            region_coord_space = data_coord_space
+        else:
+            if not isinstance(region_transform, ScaleTransform):
+                raise NotImplementedError(
+                    f"Support for reading levels with a region tranform of type "
+                    f"{type(region_transform)!r} is not yet supported."
+                )
+            # Create or check output coordinates.
+            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[misc]
+                )
+            elif len(region_coord_space) != len(data_coord_space):
+                raise ValueError(
+                    "The number of output coordinates must match the number of "
+                    "input coordinates."
+                )
+            if region_transform.output_axes != self._coord_space.axis_names:
+                raise ValueError(
+                    f"The output axes of '{region_transform.output_axes}' of the "
+                    f"region transform must match the axes "
+                    f"'{self._coord_space.axis_names}' of the coordinate space of "
+                    f"this multiscale image."
+                )
+            region_transform = group_to_level @ region_transform
+            assert isinstance(region_transform, ScaleTransform)
+
+        # Convert coordinates to new coordinate system.
+        coords, data_region, inv_transform = process_image_region(
+            region,
+            region_transform,
+            channel_coords,
+            self._schema.reference_level_properties.image_type,
+        )
+
+        # Get the array.
+        array_name = level if isinstance(level, str) else self._levels[level].name
+        try:
+            array = self[array_name]
+        except KeyError as ke:
+            raise SOMAError(
+                f"Unable to open the dense array with name '{array_name}'."
+            ) from ke
+        return somacore.SpatialRead(
+            array.read(
+                coords,
+                result_order=result_order,
+                platform_config=platform_config,
+            ),
+            data_coord_space,
+            region_coord_space,
+            inv_transform,
+        )
 
     # Metadata operations
 
@@ -169,7 +477,7 @@ def axis_names(self) -> Tuple[str, ...]:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return self._schema.axis_names
 
     @property
     def coordinate_space(self) -> CoordinateSpace:
@@ -178,7 +486,7 @@ def coordinate_space(self) -> CoordinateSpace:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return self._coord_space
 
     @coordinate_space.setter
     def coordinate_space(self, value: CoordinateSpace) -> None:
@@ -187,7 +495,17 @@ def coordinate_space(self, value: CoordinateSpace) -> None:
         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 multiscale image. Existing axis "
+                    f"names are {self._coord_space.axis_names}. New coordinate space "
+                    f"has axis names {value.axis_names}."
+                )
+        self.metadata[SOMA_COORDINATE_SPACE_METADATA_KEY] = coordinate_space_to_json(
+            value
+        )
+        self._coord_space = value
 
     def get_transform_from_level(self, level: Union[int, str]) -> ScaleTransform:
         """Returns the transformation from user requested level to image reference
@@ -196,7 +514,34 @@ def get_transform_from_level(self, level: Union[int, str]) -> ScaleTransform:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        if isinstance(level, str):
+            for val in self._levels:
+                if val.name == level:
+                    level_props = val
+                else:
+                    raise KeyError("No level with name '{level}'")
+        else:
+            level_props = self._levels[level]
+        ref_level_props = self._schema.reference_level_properties
+        if ref_level_props.depth is None:
+            return ScaleTransform(
+                input_axes=self._coord_space.axis_names,
+                output_axes=self._coord_space.axis_names,
+                scale_factors=[
+                    ref_level_props.width / level_props.width,
+                    ref_level_props.height / level_props.height,
+                ],
+            )
+        assert level_props.depth is not None
+        return ScaleTransform(
+            input_axes=self._coord_space.axis_names,
+            output_axes=self._coord_space.axis_names,
+            scale_factors=[
+                ref_level_props.width / level_props.width,
+                ref_level_props.height / level_props.height,
+                ref_level_props.depth / level_props.depth,
+            ],
+        )
 
     def get_transform_to_level(self, level: Union[int, str]) -> ScaleTransform:
         """Returns the transformation from the image reference level to the user
@@ -205,7 +550,34 @@ def get_transform_to_level(self, level: Union[int, str]) -> ScaleTransform:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        if isinstance(level, str):
+            for val in self._levels:
+                if val.name == level:
+                    level_props = val
+                else:
+                    raise KeyError("No level with name '{level}'")
+        else:
+            level_props = self._levels[level]
+        ref_level_props = self._schema.reference_level_properties
+        if ref_level_props.depth is None:
+            return ScaleTransform(
+                input_axes=self._coord_space.axis_names,
+                output_axes=self._coord_space.axis_names,
+                scale_factors=[
+                    level_props.width / ref_level_props.width,
+                    level_props.height / ref_level_props.height,
+                ],
+            )
+        assert level_props.depth is not None
+        return ScaleTransform(
+            input_axes=self._coord_space.axis_names,
+            output_axes=self._coord_space.axis_names,
+            scale_factors=[
+                level_props.width / ref_level_props.width,
+                level_props.height / ref_level_props.height,
+                level_props.depth / ref_level_props.depth,
+            ],
+        )
 
     @property
     def image_type(self) -> str:
@@ -214,7 +586,7 @@ def image_type(self) -> str:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return self._schema.reference_level_properties.image_type
 
     @property
     def level_count(self) -> int:
@@ -223,7 +595,7 @@ def level_count(self) -> int:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return len(self._levels)
 
     def level_properties(self, level: Union[int, str]) -> somacore.ImageProperties:
         """The properties of an image at the specified level.
@@ -231,7 +603,11 @@ def level_properties(self, level: Union[int, str]) -> somacore.ImageProperties:
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        if isinstance(level, str):
+            raise NotImplementedError(
+                "Support for getting level properties by name is not yet implemented."
+            )  # TODO
+        return self._levels[level]
 
     @property
     def reference_level(self) -> Optional[int]:
@@ -252,4 +628,70 @@ def reference_level_properties(self) -> "ImageProperties":
         Lifecycle:
             Experimental.
         """
-        raise NotImplementedError()
+        return self._schema.reference_level_properties
+
+
+# TODO: Push down to C++ layer
+@attrs.define
+class MultiscaleImageSchema:
+
+    reference_level_properties: ImageProperties
+    axis_names: Tuple[str, ...]
+    datatype: pa.DataType
+
+    def __attrs_post_init__(self):  # type: ignore[no-untyped-def]
+        ndim = len(self.reference_level_properties.shape)
+        if len(self.axis_names) != ndim:
+            raise ValueError(
+                f"Invalid axis names '{self.axis_names}'. {ndim} axis names must be "
+                f"provided for a multiscale image with image type "
+                f"{self.reference_level_properties.image_type}. "
+            )
+
+    # mypy false positive https://github.com/python/mypy/issues/5313
+    def create_coordinate_space(self) -> CoordinateSpace:
+        return CoordinateSpace(
+            tuple(Axis(name) for name in self.get_coordinate_space_axis_names())  # type: ignore[misc]
+        )
+
+    def get_coordinate_space_axis_names(self) -> Tuple[str, ...]:
+        # TODO: Setting axes and the coordinate space is going to be updated
+        # in a future PR.
+        x_name: Optional[str] = None
+        y_name: Optional[str] = None
+        z_name: Optional[str] = None
+        for axis_name, axis_type in zip(
+            self.axis_names, self.reference_level_properties.image_type
+        ):
+            if axis_type == "X":
+                x_name = axis_name
+            elif axis_type == "Y":
+                y_name = axis_name
+            elif axis_type == "Z":
+                z_name = axis_name
+        assert x_name is not None  # For mypy (already validated)
+        assert y_name is not None  # For mypy (already validated)
+        if z_name is None:
+            return (x_name, y_name)
+        else:
+            return (x_name, y_name, z_name)
+
+    def to_json(self) -> str:
+        type_str = pyarrow_to_carrow_type(self.datatype)
+        return json.dumps(
+            {
+                "name": self.reference_level_properties.name,
+                "image_type": self.reference_level_properties.image_type,
+                "shape": self.reference_level_properties.shape,
+                "axis_names": self.axis_names,
+                "datatype": type_str,
+            }
+        )
+
+    @classmethod
+    def from_json(cls, data: str) -> Self:
+        kwargs = json.loads(data)
+        axis_names = kwargs.pop("axis_names")
+        type_str = kwargs.pop("datatype")
+        type = carrow_type_to_pyarrow(type_str)
+        return cls(ImageProperties(**kwargs), axis_names, type)
diff --git a/apis/python/src/tiledbsoma/_soma_object.py b/apis/python/src/tiledbsoma/_soma_object.py
index 3b9132cb16..d6c306743e 100644
--- a/apis/python/src/tiledbsoma/_soma_object.py
+++ b/apis/python/src/tiledbsoma/_soma_object.py
@@ -46,6 +46,7 @@ class SOMAObject(somacore.SOMAObject, Generic[_WrapperType_co]):
         Type[_tdb_handles.CollectionWrapper],
         Type[_tdb_handles.ExperimentWrapper],
         Type[_tdb_handles.MeasurementWrapper],
+        Type[_tdb_handles.MultiscaleImageWrapper],
     ]
     """Class variable of the Wrapper class used to open this object type."""
 
diff --git a/apis/python/src/tiledbsoma/_spatial_util.py b/apis/python/src/tiledbsoma/_spatial_util.py
index d7769e66b6..57a4866cfa 100644
--- a/apis/python/src/tiledbsoma/_spatial_util.py
+++ b/apis/python/src/tiledbsoma/_spatial_util.py
@@ -112,6 +112,50 @@ def transform_region(
     return shapely.affinity.affine_transform(region, affine)
 
 
+def process_image_region(
+    region: options.SpatialRegion,
+    transform: somacore.CoordinateTransform,
+    channel_coords: options.DenseCoord,
+    image_type: str,
+) -> Tuple[options.DenseNDCoords, options.SpatialRegion, somacore.CoordinateTransform]:
+    # Get the transformed region the user is selecting in the data space.
+    # Note: transform region verifies only 2D data. This function is hard-coded to
+    # make the same assumption.
+    data_region = transform_region(region, transform)
+
+    # Convert the region to a bounding box. Round values of bounding box to integer
+    # values. Include any partially intersected pixels.
+    (x_min, y_min, x_max, y_max) = shapely.bounds(data_region)
+    x_min = max(0, int(np.floor(x_min)))
+    y_min = max(0, int(np.floor(y_min)))
+    x_max = int(np.ceil(x_max))
+    y_max = int(np.ceil(y_max))
+
+    # Get the inverse translation from the data space to the original requested region.
+    if x_min != 0 or y_min != 0:
+        translate = somacore.AffineTransform(
+            transform.output_axes,
+            transform.output_axes,
+            np.array([[1, 0, -x_min], [0, 1, -y_min], [0, 0, 1]]),
+        )
+        transform = translate @ transform
+    inv_transform = transform.inverse_transform()
+
+    coords: options.DenseNDCoords = []
+    for axis in image_type:
+        if axis == "C":
+            coords.append(channel_coords)  # type: ignore[attr-defined]
+        if axis == "X":
+            coords.append(slice(x_min, x_max))  # type: ignore[attr-defined]
+        if axis == "Y":
+            coords.append(slice(y_min, y_max))  # type: ignore[attr-defined]
+        if axis == "Z":
+            raise NotImplementedError(
+                "Spatial queries are currently only supported for 2D coordinates."
+            )
+    return (coords, data_region, inv_transform)
+
+
 def process_spatial_df_region(
     region: Optional[options.SpatialRegion],
     transform: somacore.CoordinateTransform,
diff --git a/apis/python/src/tiledbsoma/_tdb_handles.py b/apis/python/src/tiledbsoma/_tdb_handles.py
index 2284636c82..f910502801 100644
--- a/apis/python/src/tiledbsoma/_tdb_handles.py
+++ b/apis/python/src/tiledbsoma/_tdb_handles.py
@@ -48,6 +48,7 @@
     clib.SOMACollection,
     clib.SOMAMeasurement,
     clib.SOMAExperiment,
+    clib.SOMAMultiscaleImage,
 ]
 _RawHdl_co = TypeVar("_RawHdl_co", bound=RawHandle, covariant=True)
 """A raw TileDB object. Covariant because Handles are immutable enough."""
@@ -84,6 +85,7 @@ def open(
         "somacollection": CollectionWrapper,
         "somaexperiment": ExperimentWrapper,
         "somameasurement": MeasurementWrapper,
+        "somamultiscaleimage": MultiscaleImageWrapper,
     }
 
     try:
@@ -91,12 +93,7 @@ def open(
             soma_object, context
         )
     except KeyError:
-        if soma_object.type.lower() in {
-            "somascene",
-            "somapointclouddataframe",
-            "somageometrydataframe",
-            "somamultiscaleimage",
-        }:
+        if soma_object.type.lower() in {"somascene", "somageometrydataframe"}:
             raise NotImplementedError(
                 f"Support for {soma_object.type!r} is not yet implemented."
             )
@@ -323,6 +320,12 @@ class MeasurementWrapper(SOMAGroupWrapper[clib.SOMAMeasurement]):
     _GROUP_WRAPPED_TYPE = clib.SOMAMeasurement
 
 
+class MultiscaleImageWrapper(SOMAGroupWrapper[clib.SOMAMultiscaleImage]):
+    """Wrapper around a Pybind11 MultiscaleImage handle."""
+
+    _GROUP_WRAPPED_TYPE = clib.SOMAMultiscaleImage
+
+
 _ArrType = TypeVar("_ArrType", bound=clib.SOMAArray)
 
 
diff --git a/apis/python/src/tiledbsoma/soma_collection.cc b/apis/python/src/tiledbsoma/soma_collection.cc
index 27a53a152f..38e6522e32 100644
--- a/apis/python/src/tiledbsoma/soma_collection.cc
+++ b/apis/python/src/tiledbsoma/soma_collection.cc
@@ -74,5 +74,8 @@ void load_soma_collection(py::module& m) {
 
     py::class_<SOMAMeasurement, SOMACollection, SOMAGroup, SOMAObject>(
         m, "SOMAMeasurement");
+
+    py::class_<SOMAMultiscaleImage, SOMACollection, SOMAGroup, SOMAObject>(
+        m, "SOMAMultiscaleImage");
 }
 }  // namespace libtiledbsomacpp
diff --git a/apis/python/src/tiledbsoma/soma_object.cc b/apis/python/src/tiledbsoma/soma_object.cc
index 61561f0756..09b6a54048 100644
--- a/apis/python/src/tiledbsoma/soma_object.cc
+++ b/apis/python/src/tiledbsoma/soma_object.cc
@@ -91,6 +91,9 @@ void load_soma_object(py::module& m) {
                     else if (soma_obj_type == "somameasurement")
                         return py::cast(
                             dynamic_cast<SOMAMeasurement&>(*soma_obj));
+                    else if (soma_obj_type == "somamultiscaleimage")
+                        return py::cast(
+                            dynamic_cast<SOMAMultiscaleImage&>(*soma_obj));
                     return py::none();
                 } catch (...) {
                     return py::none();
diff --git a/apis/python/tests/test_multiscale_image.py b/apis/python/tests/test_multiscale_image.py
new file mode 100644
index 0000000000..7ac8331eb3
--- /dev/null
+++ b/apis/python/tests/test_multiscale_image.py
@@ -0,0 +1,122 @@
+from urllib.parse import urljoin
+
+import numpy as np
+import pyarrow as pa
+import pytest
+
+import tiledbsoma as soma
+
+
+def test_multiscale_image_bad_create(tmp_path):
+    baseuri = urljoin(f"{tmp_path.as_uri()}/", "bad_create")
+
+    # Invalid datetype.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.string(),
+            axis_names=("x", "y", "y"),
+            reference_level_shape=(3, 64, 64),
+        )
+
+    # Repeated axis names.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.uint8(),
+            axis_names=("x", "y", "y"),
+            reference_level_shape=(3, 64, 64),
+        )
+
+    # Repeated axis types.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.uint8(),
+            axis_types=("height", "height", "width"),
+            reference_level_shape=(3, 64, 64),
+        )
+
+    # Invalid axis type.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.uint8(),
+            axis_types=("c", "height", "width"),
+            reference_level_shape=(3, 64, 64),
+        )
+
+    # Mismatch in number of axis names and reference shape.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.uint8(),
+            axis_types=("x", "y", "y"),
+            reference_level_shape=(3, 64, 64),
+        )
+
+    # Mismatch in number of axis names and axis types.
+    with pytest.raises(ValueError):
+        soma.MultiscaleImage.create(
+            baseuri,
+            type=pa.uint8(),
+            reference_level_shape=(64, 64),
+            axis_names=("y", "x"),
+            axis_types=("channels",),
+        )
+
+
+def test_multiscale_basic_no_channels(tmp_path):
+    baseuri = urljoin(f"{tmp_path.as_uri()}/", "basic_read")
+    image_uri = urljoin(baseuri, "default")
+
+    # Create the multiscale image.
+    with soma.MultiscaleImage.create(
+        image_uri,
+        type=pa.uint8(),
+        reference_level_shape=(128, 64),
+        axis_names=("y", "x"),
+        axis_types=("height", "width"),
+    ) as image:
+
+        # Create reference level.
+        image.add_new_level("level0", shape=(128, 64))
+
+        # Create medium sized downsample.
+        image.add_new_level("level1", shape=(64, 32))
+
+        # Create very small downsample and write to it.
+        level2 = image.add_new_level("level2", shape=(8, 4))
+        data = pa.Tensor.from_numpy(np.arange(32, dtype=np.uint8))
+        level2.write((slice(None), slice(None)), data)
+
+    # Open for reading and check metadata.
+    with soma.MultiscaleImage.open(image_uri, mode="r") as image:
+
+        # Check the base properties for the image.
+        # assert image.axis_names == ("y", "x")
+        base_props = image.reference_level_properties
+        assert base_props.name == "reference_level"
+        assert base_props.shape == (128, 64)
+        assert base_props.height == 128
+        assert base_props.width == 64
+        assert base_props.nchannels is None
+        assert base_props.depth is None
+        assert base_props.image_type == "YX"
+
+        # Check coordinate space.
+        coord_space = image.coordinate_space
+        assert len(coord_space) == 2
+        assert coord_space.axis_names == ("x", "y")
+
+        # Check the number of levels and level properties.
+        assert image.level_count == 3
+        for index, shape in enumerate([(128, 64), (64, 32), (8, 4)]):
+            props = image.level_properties(index)
+            assert props.nchannels is None
+            assert props.depth is None
+            assert props.image_type == "YX"
+            assert props.name == f"level{index}"
+            assert props.shape == shape
+            assert props.height == shape[0]
+            assert props.width == shape[1]