Repository: shapely/shapely (language: python)

Issue:
Expose method to get / construct from coordinates and offset arrays
In the GeoArrow initiative, we defined a way how to store geometries in a Arrow compatible memory format (using a more "native" storage of the data (using the nested types available in Arrow) compared to using binary WKB): https://github.com/geopandas/geo-arrow-spec/blob/main/format.md (https://github.com/geopandas/geo-arrow-spec/issues/4)
Essentially, this uses a "nested list" logical type (like GeoJSON), or "ragged array", which physically means a single flat array of coordinates with zero or more additional offset arrays to keep track of where each new ring, new geometry or new geometry part starts in the coordinate array. 

This was discussed before in pygeos as well, see:

* [ENH: get coordinates + offset arrays for "ragged array" representation - pygeos/pygeos#75](https://github.com/pygeos/pygeos/issues/75)
* [Add cython algo to get offset arrays for flat coordinates representation - pygeos/pygeos#93 ](https://github.com/pygeos/pygeos/pull/93)
* [Add “offsets” as an alternative to “indices” - pygeos/pygeos#438](https://github.com/pygeos/pygeos/issues/438)


Discussion / hints:
My current code for this (using pygeos, but the same with shapely):

<details>

```python
def _get_arrays_point(arr):
    # only one array of coordinates
    coords = pygeos.get_coordinates(arr)

    return coords.ravel(), None


def _get_arrays_multipoint(arr):
    # explode/flatten the MultiPoints
    arr_flat, part_indices = pygeos.get_parts(arr, return_index=True)
    # the offsets into the multipoint parts
    offsets = np.insert(np.bincount(part_indices).cumsum(), 0, 0)

    # only one array of coordinates
    coords = pygeos.get_coordinates(arr)

    return coords.ravel(), offsets


def _get_arrays_linestring(arr):
    # the coords and offsets into the coordinates of the linestrings
    coords, indices = pygeos.get_coordinates(arr, return_index=True)
    offsets = np.insert(np.bincount(indices).cumsum(), 0, 0)

    return coords.ravel(), offsets


def _get_arrays_multilinestring(arr):
    # explode/flatten the MultiLineStrings
    arr_flat, part_indices = pygeos.get_parts(arr, return_index=True)
    # the offsets into the multilinestring parts
    offsets2 = np.insert(np.bincount(part_indices).cumsum(), 0, 0)

    # the coords and offsets into the coordinates of the linestrings
    coords, indices = pygeos.get_coordinates(arr_flat, return_index=True)
    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)

    return coords.ravel(), (offsets1, offsets2)


def _get_arrays_polygon(arr):
    # explode/flatten the Polygons into Rings
    arr_flat2, ring_indices = pygeos.geometry.get_rings(arr, return_index=True)
    # the offsets into the exterior/interior rings of the multipolygon parts
    offsets2 = np.insert(np.bincount(ring_indices).cumsum(), 0, 0)

    # the coords and offsets into the coordinates of the rings
    coords, indices = pygeos.get_coordinates(arr_flat2, return_index=True)
    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)

    return coords.ravel(), (offsets1, offsets2)


def _get_arrays_multipolygon(arr):
    # explode/flatten the MultiPolygons
    arr_flat, part_indices = pygeos.get_parts(arr, return_index=True)
    # the offsets into the multipolygon parts
    offsets3 = np.insert(np.bincount(part_indices).cumsum(), 0, 0)

    # explode/flatten the Polygons into Rings
    arr_flat2, ring_indices = pygeos.geometry.get_rings(arr_flat, return_index=True)
    # the offsets into the exterior/interior rings of the multipolygon parts
    offsets2 = np.insert(np.bincount(ring_indices).cumsum(), 0, 0)

    # the coords and offsets into the coordinates of the rings
    coords, indices = pygeos.get_coordinates(arr_flat2, return_index=True)
    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)

    return coords.ravel(), (offsets1, offsets2, offsets3)


def get_flat_coords_offset_arrays(arr):
    geom_types = np.unique(pygeos.get_type_id(arr))
    # ignore missing values (type of -1)
    geom_types = geom_types[geom_types >= 0]

    if len(geom_types) == 1 and geom_types[0] == GeometryType.POINT:
        typ = "point"
        coords, offsets = _get_arrays_point(arr)

    elif len(geom_types) == 1 and geom_types[0] == GeometryType.LINESTRING:
        typ = "linestring"
        coords, offsets = _get_arrays_linestring(arr)

    elif len(geom_types) == 1 and geom_types[0] == GeometryType.POLYGON:
        typ = "polygon"
        coords, offsets = _get_arrays_polygon(arr)

    elif all(t in {GeometryType.POINT, GeometryType.MULTIPOINT} for t in geom_types):
        typ = "multipoint"
        coords, offsets = _get_arrays_multipoint(arr)

    elif all(
        t in {GeometryType.LINESTRING, GeometryType.MULTILINESTRING} for t in geom_types
    ):
        typ = "multilinestring"
        coords, offsets = _get_arrays_multilinestring(arr)

    elif all(
        t in {GeometryType.POLYGON, GeometryType.MULTIPOLYGON} for t in geom_types
    ):
        typ = "multipolygon"
        coords, offsets = _get_arrays_multipolygon(arr)

    else:
        raise ValueError(
            f"Geometry type combination is not supported ({list(geom_types)})"
        )

    return typ, coords, offsets


# # coords/offset arrays -> GEOS


def _point_from_flatcoords(coords):
    result = pygeos.points(coords.reshape(-1, 2))

    return result


def _multipoint_from_flatcoords(coords, offsets):
    # recreate points
    points = pygeos.points(coords.reshape(-1, 2))

    # recreate multipoints
    multipoint_parts = np.diff(offsets)
    multipoint_indices = np.repeat(np.arange(len(multipoint_parts)), multipoint_parts)
    result = pygeos.multipoints(points, indices=multipoint_indices)

    return result


def _linestring_from_flatcoords(coords, offsets):
    # recreate linestrings
    linestring_n = np.diff(offsets)
    linestring_indices = np.repeat(np.arange(len(linestring_n)), linestring_n)
    result = pygeos.linestrings(coords.reshape(-1, 2), indices=linestring_indices)

    return result


def _multilinestrings_from_flatcoords(coords, offsets1, offsets2):
    # recreate linestrings
    linestring_n = np.diff(offsets1)
    linestring_indices = np.repeat(np.arange(len(linestring_n)), linestring_n)
    linestrings = pygeos.linestrings(coords.reshape(-1, 2), indices=linestring_indices)

    # recreate multilinestrings
    multilinestring_parts = np.diff(offsets2)
    multilinestring_indices = np.repeat(
        np.arange(len(multilinestring_parts)), multilinestring_parts
    )
    result = pygeos.multilinestrings(linestrings, indices=multilinestring_indices)

    return result


def _polygon_from_flatcoords(coords, offsets1, offsets2):
    # recreate rings
    ring_lengths = np.diff(offsets1)
    ring_indices = np.repeat(np.arange(len(ring_lengths)), ring_lengths)
    rings = pygeos.linearrings(coords.reshape(-1, 2), indices=ring_indices)

    # recreate polygons
    polygon_rings_n = np.diff(offsets2)
    polygon_indices = np.repeat(np.arange(len(polygon_rings_n)), polygon_rings_n)
    result = pygeos.polygons(rings, indices=polygon_indices)

    return result


def _multipolygons_from_flatcoords(coords, offsets1, offsets2, offsets3):
    # recreate polygons
    polygons = _polygon_from_flatcoords(coords, offsets1, offsets2)

    # recreate multipolygons
    multipolygon_parts = np.diff(offsets3)
    multipolygon_indices = np.repeat(
        np.arange(len(multipolygon_parts)), multipolygon_parts
    )
    result = pygeos.multipolygons(polygons, indices=multipolygon_indices)

    return result


def get_geometries_from_flatcoords(typ, coords, offsets):

    if typ == "point":
        return _point_from_flatcoords(coords)
    if typ == "linestring":
        return _linestring_from_flatcoords(coords, offsets)
    if typ == "polygon":
        return _polygon_from_flatcoords(coords, *offsets)
    elif typ == "multipoint":
        return _multipoint_from_flatcoords(coords, offsets)
    elif typ == "multilinestring":
        return _multilinestrings_from_flatcoords(coords, *offsets)
    elif typ == "multipolygon":
        return _multipolygons_from_flatcoords(coords, *offsets)
    else:
        raise ValueError(typ)
```

</details>


Fix patch:
diff --git a/docs/release/2.x.rst b/docs/release/2.x.rst
index fbd15258d..c2ddfbd6d 100644
--- a/docs/release/2.x.rst
+++ b/docs/release/2.x.rst
@@ -391,12 +391,14 @@ prepared. Helper functions :func:`.destroy_prepared` and :func:`.is_prepared`
 are also available.
 
 
-GeoJSON IO
-~~~~~~~~~~
+New IO methods (GeoJSON, ragged arrays)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 * Added GeoJSON input/output capabilities :func:`.from_geojson` and
-  :func:`.to_geojson` for GEOS >= 3.10 (#413).
-
+  :func:`.to_geojson` for GEOS >= 3.10.
+* Added conversion to/from ragged array representation using a contiguous array
+  of coordinates and offset arrays: :func:`.to_ragged_array` and
+  :func:`.from_ragged_array`.
 
 Other improvements
 ~~~~~~~~~~~~~~~~~~
diff --git a/shapely/_ragged_array.py b/shapely/_ragged_array.py
new file mode 100644
index 000000000..68e616da9
--- /dev/null
+++ b/shapely/_ragged_array.py
@@ -0,0 +1,452 @@
+"""
+This modules provides a conversion to / from a ragged (or "jagged") array
+representation of the geometries.
+
+A ragged array is an irregular array of arrays of which each element can have
+a different length. As a result, such an array cannot be represented as a
+standard, rectangular nD array.
+The coordinates of geometries can be represented as arrays of arrays of
+coordinate pairs (possibly multiple levels of nesting, depending on the
+geometry type).
+
+
+Geometries, as a ragged array of coordinates, can be efficiently represented
+as contiguous arrays of coordinates provided that there is another data
+structure that keeps track of which range of coordinate values corresponds
+to a given geometry. This can be done using offsets, counts, or indices.
+
+This module currently implements offsets into the coordinates array. This
+is the ragged array representation defined by the the Apache Arrow project
+as "variable size list array" (https://arrow.apache.org/docs/format/Columnar.html#variable-size-list-layout).
+See for example https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#representations-features
+for different options.
+
+The exact usage of the Arrow list array with varying degrees of nesting for the
+different geometry types is defined by the GeoArrow project:
+https://github.com/geoarrow/geoarrow
+
+"""
+import numpy as np
+
+from . import creation
+from ._geometry import (
+    GeometryType,
+    get_coordinate_dimension,
+    get_parts,
+    get_rings,
+    get_type_id,
+)
+from .coordinates import get_coordinates
+from .predicates import is_empty
+
+__all__ = ["to_ragged_array", "from_ragged_array"]
+
+
+# # GEOS -> coords/offset arrays (to_ragged_array)
+
+
+def _get_arrays_point(arr, include_z):
+    # only one array of coordinates
+    coords = get_coordinates(arr, include_z=include_z)
+
+    # empty points are represented by NaNs
+    empties = is_empty(arr)
+    if empties.any():
+        indices = np.nonzero(empties)[0]
+        indices = indices - np.arange(len(indices))
+        coords = np.insert(coords, indices, np.nan, axis=0)
+
+    return coords, ()
+
+
+def _indices_to_offsets(indices, n):
+    offsets = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    if len(offsets) != n + 1:
+        # last geometries might be empty or missing
+        offsets = np.pad(
+            offsets,
+            (0, n + 1 - len(offsets)),
+            "constant",
+            constant_values=offsets[-1],
+        )
+    return offsets
+
+
+def _get_arrays_multipoint(arr, include_z):
+    # explode/flatten the MultiPoints
+    _, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multipoint parts
+    offsets = _indices_to_offsets(part_indices, len(arr))
+
+    # only one array of coordinates
+    coords = get_coordinates(arr, include_z=include_z)
+
+    return coords, (offsets,)
+
+
+def _get_arrays_linestring(arr, include_z):
+    # the coords and offsets into the coordinates of the linestrings
+    coords, indices = get_coordinates(arr, return_index=True, include_z=include_z)
+    offsets = _indices_to_offsets(indices, len(arr))
+
+    return coords, (offsets,)
+
+
+def _get_arrays_multilinestring(arr, include_z):
+    # explode/flatten the MultiLineStrings
+    arr_flat, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multilinestring parts
+    offsets2 = _indices_to_offsets(part_indices, len(arr))
+
+    # the coords and offsets into the coordinates of the linestrings
+    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2)
+
+
+def _get_arrays_polygon(arr, include_z):
+    # explode/flatten the Polygons into Rings
+    arr_flat, ring_indices = get_rings(arr, return_index=True)
+    # the offsets into the exterior/interior rings of the multipolygon parts
+    offsets2 = _indices_to_offsets(ring_indices, len(arr))
+
+    # the coords and offsets into the coordinates of the rings
+    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2)
+
+
+def _get_arrays_multipolygon(arr, include_z):
+    # explode/flatten the MultiPolygons
+    arr_flat, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multipolygon parts
+    offsets3 = _indices_to_offsets(part_indices, len(arr))
+
+    # explode/flatten the Polygons into Rings
+    arr_flat2, ring_indices = get_rings(arr_flat, return_index=True)
+    # the offsets into the exterior/interior rings of the multipolygon parts
+    offsets2 = np.insert(np.bincount(ring_indices).cumsum(), 0, 0)
+
+    # the coords and offsets into the coordinates of the rings
+    coords, indices = get_coordinates(arr_flat2, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2, offsets3)
+
+
+def to_ragged_array(geometries, include_z=None):
+    """
+    Converts geometries to a ragged array representation using a contiguous
+    array of coordinates and offset arrays.
+
+    This function converts an array of geometries to a ragged array
+    (i.e. irregular array of arrays) of coordinates, represented in memory
+    using a single contiguous array of the coordinates, and
+    up to 3 offset arrays that keep track where each sub-array
+    starts and ends.
+
+    This follows the in-memory layout of the variable size list arrays defined
+    by Apache Arrow, as specified for geometries by the GeoArrow project:
+    https://github.com/geoarrow/geoarrow.
+
+    Parameters
+    ----------
+    geometries : array_like
+        Array of geometries (1-dimensional).
+    include_z : bool, default None
+        If False, return 2D geometries. If True, include the third dimension
+        in the output (if a geometry has no third dimension, the z-coordinates
+        will be NaN). By default, will infer the dimensionality from the
+        input geometries. Note that this inference can be unreliable with
+        empty geometries (for a guaranteed result, it is recommended to
+        specify the keyword).
+
+    Returns
+    -------
+    tuple of (geometry_type, coords, offsets)
+        geometry_type : GeometryType
+            The type of the input geometries (required information for
+            roundtrip).
+        coords : np.ndarray
+            Contiguous array of shape (n, 2) or (n, 3) of all coordinates
+            of all input geometries.
+        offsets: tuple of np.ndarray
+            Offset arrays that make it possible to reconstruct the
+            geometries from the flat coordinates array. The number of
+            offset arrays depends on the geometry type. See
+            https://github.com/geoarrow/geoarrow/blob/main/format.md
+            for details.
+
+    Notes
+    -----
+    Mixed singular and multi geometry types of the same basic type are
+    allowed (e.g., Point and MultiPoint) and all singular types will be
+    treated as multi types.
+    GeometryCollections and other mixed geometry types are not supported.
+
+    See also
+    --------
+    from_ragged_array
+
+    Examples
+    --------
+    Consider a Polygon with one hole (interior ring):
+
+    >>> import shapely
+    >>> polygon = shapely.Polygon(
+    ...     [(0, 0), (10, 0), (10, 10), (0, 10)],
+    ...     holes=[[(2, 2), (3, 2), (2, 3)]]
+    ... )
+    >>> polygon
+    <POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (2 2, 3 2, 2 3, 2 2))>
+
+    This polygon can be thought of as a list of rings (first ring is the
+    exterior ring, subsequent rings are the interior rings), and each ring
+    as a list of coordinate pairs. This is very similar to how GeoJSON
+    represents the coordinates:
+
+    >>> import json
+    >>> json.loads(shapely.to_geojson(polygon))["coordinates"]
+    [[[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+     [[2.0, 2.0], [3.0, 2.0], [2.0, 3.0], [2.0, 2.0]]]
+
+    This function will return a similar list of lists of lists, but
+    using a single contiguous array of coordinates, and multiple arrays of
+    offsets:
+
+    >>> geometry_type, coords, offsets = shapely.to_ragged_array([polygon])
+    >>> geometry_type
+    <GeometryType.POLYGON: 3>
+    >>> coords
+    array([[ 0.,  0.],
+           [10.,  0.],
+           [10., 10.],
+           [ 0., 10.],
+           [ 0.,  0.],
+           [ 2.,  2.],
+           [ 3.,  2.],
+           [ 2.,  3.],
+           [ 2.,  2.]])
+
+    >>> offsets
+    (array([0, 5, 9]), array([0, 2]))
+
+    As an example how to interpret the offsets: the i-th ring in the
+    coordinates is represented by ``offsets[0][i]`` to ``offsets[0][i+1]``:
+
+    >>> exterior_ring_start, exterior_ring_end = offsets[0][0], offsets[0][1]
+    >>> coords[exterior_ring_start:exterior_ring_end]
+    array([[ 0.,  0.],
+           [10.,  0.],
+           [10., 10.],
+           [ 0., 10.],
+           [ 0.,  0.]])
+
+    """
+    geometries = np.asarray(geometries)
+    if include_z is None:
+        include_z = np.any(
+            get_coordinate_dimension(geometries[~is_empty(geometries)]) == 3
+        )
+
+    geom_types = np.unique(get_type_id(geometries))
+    # ignore missing values (type of -1)
+    geom_types = geom_types[geom_types >= 0]
+
+    if len(geom_types) == 1:
+        typ = GeometryType(geom_types[0])
+        if typ == GeometryType.POINT:
+            coords, offsets = _get_arrays_point(geometries, include_z)
+        elif typ == GeometryType.LINESTRING:
+            coords, offsets = _get_arrays_linestring(geometries, include_z)
+        elif typ == GeometryType.POLYGON:
+            coords, offsets = _get_arrays_polygon(geometries, include_z)
+        elif typ == GeometryType.MULTIPOINT:
+            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+        elif typ == GeometryType.MULTILINESTRING:
+            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+        elif typ == GeometryType.MULTIPOLYGON:
+            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+        else:
+            raise ValueError(f"Geometry type {typ.name} is not supported")
+
+    elif len(geom_types) == 2:
+        if set(geom_types) == {GeometryType.POINT, GeometryType.MULTIPOINT}:
+            typ = GeometryType.MULTIPOINT
+            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+        elif set(geom_types) == {GeometryType.LINESTRING, GeometryType.MULTILINESTRING}:
+            typ = GeometryType.MULTILINESTRING
+            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+        elif set(geom_types) == {GeometryType.POLYGON, GeometryType.MULTIPOLYGON}:
+            typ = GeometryType.MULTIPOLYGON
+            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+        else:
+            raise ValueError(
+                "Geometry type combination is not supported "
+                f"({[GeometryType(t).name for t in geom_types]})"
+            )
+    else:
+        raise ValueError(
+            "Geometry type combination is not supported "
+            f"({[GeometryType(t).name for t in geom_types]})"
+        )
+
+    return typ, coords, offsets
+
+
+# # coords/offset arrays -> GEOS (from_ragged_array)
+
+
+def _point_from_flatcoords(coords):
+    result = creation.points(coords)
+
+    # Older versions of GEOS (<= 3.9) don't automatically convert NaNs
+    # to empty points -> do manually
+    empties = np.isnan(coords).all(axis=1)
+    if empties.any():
+        result[empties] = creation.empty(1, geom_type=GeometryType.POINT).item()
+
+    return result
+
+
+def _multipoint_from_flatcoords(coords, offsets):
+    # recreate points
+    points = creation.points(coords)
+
+    # recreate multipoints
+    multipoint_parts = np.diff(offsets)
+    multipoint_indices = np.repeat(np.arange(len(multipoint_parts)), multipoint_parts)
+
+    result = np.empty(len(offsets) - 1, dtype=object)
+    result = creation.multipoints(points, indices=multipoint_indices, out=result)
+    result[multipoint_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTIPOINT
+    ).item()
+
+    return result
+
+
+def _linestring_from_flatcoords(coords, offsets):
+    # recreate linestrings
+    linestring_n = np.diff(offsets)
+    linestring_indices = np.repeat(np.arange(len(linestring_n)), linestring_n)
+
+    result = np.empty(len(offsets) - 1, dtype=object)
+    result = creation.linestrings(coords, indices=linestring_indices, out=result)
+    result[linestring_n == 0] = creation.empty(
+        1, geom_type=GeometryType.LINESTRING
+    ).item()
+    return result
+
+
+def _multilinestrings_from_flatcoords(coords, offsets1, offsets2):
+    # recreate linestrings
+    linestrings = _linestring_from_flatcoords(coords, offsets1)
+
+    # recreate multilinestrings
+    multilinestring_parts = np.diff(offsets2)
+    multilinestring_indices = np.repeat(
+        np.arange(len(multilinestring_parts)), multilinestring_parts
+    )
+
+    result = np.empty(len(offsets2) - 1, dtype=object)
+    result = creation.multilinestrings(
+        linestrings, indices=multilinestring_indices, out=result
+    )
+    result[multilinestring_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTILINESTRING
+    ).item()
+
+    return result
+
+
+def _polygon_from_flatcoords(coords, offsets1, offsets2):
+    # recreate rings
+    ring_lengths = np.diff(offsets1)
+    ring_indices = np.repeat(np.arange(len(ring_lengths)), ring_lengths)
+    rings = creation.linearrings(coords, indices=ring_indices)
+
+    # recreate polygons
+    polygon_rings_n = np.diff(offsets2)
+    polygon_indices = np.repeat(np.arange(len(polygon_rings_n)), polygon_rings_n)
+    result = np.empty(len(offsets2) - 1, dtype=object)
+    result = creation.polygons(rings, indices=polygon_indices, out=result)
+    result[polygon_rings_n == 0] = creation.empty(
+        1, geom_type=GeometryType.POLYGON
+    ).item()
+
+    return result
+
+
+def _multipolygons_from_flatcoords(coords, offsets1, offsets2, offsets3):
+    # recreate polygons
+    polygons = _polygon_from_flatcoords(coords, offsets1, offsets2)
+
+    # recreate multipolygons
+    multipolygon_parts = np.diff(offsets3)
+    multipolygon_indices = np.repeat(
+        np.arange(len(multipolygon_parts)), multipolygon_parts
+    )
+    result = np.empty(len(offsets3) - 1, dtype=object)
+    result = creation.multipolygons(polygons, indices=multipolygon_indices, out=result)
+    result[multipolygon_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTIPOLYGON
+    ).item()
+
+    return result
+
+
+def from_ragged_array(geometry_type, coords, offsets=None):
+    """
+    Creates geometries from a contiguous array of coordinates
+    and offset arrays.
+
+    This function creates geometries from the ragged array representation
+    as returned by ``to_ragged_array``.
+
+    This follows the in-memory layout of the variable size list arrays defined
+    by Apache Arrow, as specified for geometries by the GeoArrow project:
+    https://github.com/geoarrow/geoarrow.
+
+    See :func:`to_ragged_array` for more details.
+
+    Parameters
+    ----------
+    geometry_type : GeometryType
+        The type of geometry to create.
+    coords : np.ndarray
+        Contiguous array of shape (n, 2) ro (n, 3) of all coordinates
+        for the geometries.
+    offsets: tuple of np.ndarray
+        Offset arrays that allow to reconstruct the geometries based on the
+        flat coordinates array. The number of offset arrays depends on the
+        geometry type. See
+        https://github.com/geoarrow/geoarrow/blob/main/format.md for details.
+
+    Returns
+    -------
+    np.ndarray
+        Array of geometries (1-dimensional).
+
+    See Also
+    --------
+    to_ragged_array
+
+    """
+    if geometry_type == GeometryType.POINT:
+        assert offsets is None or len(offsets) == 0
+        return _point_from_flatcoords(coords)
+    if geometry_type == GeometryType.LINESTRING:
+        return _linestring_from_flatcoords(coords, *offsets)
+    if geometry_type == GeometryType.POLYGON:
+        return _polygon_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTIPOINT:
+        return _multipoint_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTILINESTRING:
+        return _multilinestrings_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTIPOLYGON:
+        return _multipolygons_from_flatcoords(coords, *offsets)
+    else:
+        raise ValueError(f"Geometry type {geometry_type.name} is not supported")
diff --git a/shapely/io.py b/shapely/io.py
index a6159cdc7..38631cfe5 100644
--- a/shapely/io.py
+++ b/shapely/io.py
@@ -1,14 +1,19 @@
 import numpy as np
 
 from . import lib
+
+# include ragged array functions here for reference documentation purpose
+from ._ragged_array import from_ragged_array, to_ragged_array
 from .decorators import requires_geos
 from .enum import ParamEnum
 
 __all__ = [
     "from_geojson",
+    "from_ragged_array",
     "from_wkb",
     "from_wkt",
     "to_geojson",
+    "to_ragged_array",
     "to_wkb",
     "to_wkt",
 ]


Test patch:
diff --git a/shapely/tests/test_ragged_array.py b/shapely/tests/test_ragged_array.py
new file mode 100644
index 000000000..2fbc2dc3f
--- /dev/null
+++ b/shapely/tests/test_ragged_array.py
@@ -0,0 +1,362 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+import shapely
+from shapely import MultiLineString, MultiPoint, MultiPolygon
+from shapely.testing import assert_geometries_equal
+
+from .common import (
+    empty_line_string,
+    empty_line_string_z,
+    geometry_collection,
+    line_string,
+    line_string_z,
+    linear_ring,
+    multi_line_string,
+    multi_line_string_z,
+    multi_point,
+    multi_point_z,
+    multi_polygon,
+    multi_polygon_z,
+    point,
+    point_z,
+    polygon,
+    polygon_z,
+)
+
+all_types = (
+    point,
+    line_string,
+    polygon,
+    multi_point,
+    multi_line_string,
+    multi_polygon,
+)
+
+all_types_3d = (
+    point_z,
+    line_string_z,
+    polygon_z,
+    multi_point_z,
+    multi_line_string_z,
+    multi_polygon_z,
+)
+
+all_types_not_supported = (
+    linear_ring,
+    geometry_collection,
+)
+
+
+@pytest.mark.parametrize("geom", all_types + all_types_3d)
+def test_roundtrip(geom):
+    actual = shapely.from_ragged_array(*shapely.to_ragged_array([geom, geom]))
+    assert_geometries_equal(actual, [geom, geom])
+
+
+@pytest.mark.parametrize("geom", all_types)
+def test_include_z(geom):
+    _, coords, _ = shapely.to_ragged_array([geom, geom], include_z=True)
+    # For 2D geoms, z coords are filled in with NaN
+    assert np.isnan(coords[:, 2]).all()
+
+
+@pytest.mark.parametrize("geom", all_types_3d)
+def test_include_z_false(geom):
+    _, coords, _ = shapely.to_ragged_array([geom, geom], include_z=False)
+    # For 3D geoms, z coords are dropped
+    assert coords.shape[1] == 2
+
+
+def test_include_z_default():
+    # corner cases for inferring dimensionality
+
+    # mixed 2D and 3D -> 3D
+    _, coords, _ = shapely.to_ragged_array([line_string, line_string_z])
+    assert coords.shape[1] == 3
+
+    # only empties -> always 2D
+    _, coords, _ = shapely.to_ragged_array([empty_line_string])
+    assert coords.shape[1] == 2
+    _, coords, _ = shapely.to_ragged_array([empty_line_string_z])
+    assert coords.shape[1] == 2
+    # empty collection -> GEOS indicates 2D
+    _, coords, _ = shapely.to_ragged_array(shapely.from_wkt(["MULTIPOLYGON Z EMPTY"]))
+    assert coords.shape[1] == 2
+
+
+@pytest.mark.parametrize("geom", all_types_not_supported)
+def test_raise_geometry_type(geom):
+    with pytest.raises(ValueError):
+        shapely.to_ragged_array([geom, geom])
+
+
+def test_points():
+    arr = shapely.from_wkt(
+        [
+            "POINT (0 0)",
+            "POINT (1 1)",
+            "POINT EMPTY",
+            "POINT EMPTY",
+            "POINT (4 4)",
+            "POINT EMPTY",
+        ]
+    )
+    typ, result, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [[0, 0], [1, 1], [np.nan, np.nan], [np.nan, np.nan], [4, 4], [np.nan, np.nan]]
+    )
+    assert typ == shapely.GeometryType.POINT
+    assert_allclose(result, expected)
+    assert len(offsets) == 0
+
+    geoms = shapely.from_ragged_array(typ, result)
+    assert_geometries_equal(geoms, arr)
+
+
+def test_linestrings():
+    arr = shapely.from_wkt(
+        [
+            "LINESTRING (30 10, 10 30, 40 40)",
+            "LINESTRING (40 40, 30 30, 40 20, 30 10)",
+            "LINESTRING EMPTY",
+            "LINESTRING EMPTY",
+            "LINESTRING (10 10, 20 20, 10 40)",
+            "LINESTRING EMPTY",
+        ]
+    )
+    typ, coords, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [
+            [30.0, 10.0],
+            [10.0, 30.0],
+            [40.0, 40.0],
+            [40.0, 40.0],
+            [30.0, 30.0],
+            [40.0, 20.0],
+            [30.0, 10.0],
+            [10.0, 10.0],
+            [20.0, 20.0],
+            [10.0, 40.0],
+        ]
+    )
+    expected_offsets = np.array([0, 3, 7, 7, 7, 10, 10])
+    assert typ == shapely.GeometryType.LINESTRING
+    assert_allclose(coords, expected)
+    assert len(offsets) == 1
+    assert_allclose(offsets[0], expected_offsets)
+
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    assert_geometries_equal(result, arr)
+
+
+def test_polygons():
+    arr = shapely.from_wkt(
+        [
+            "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))",
+            "POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))",
+            "POLYGON EMPTY",
+            "POLYGON EMPTY",
+            "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))",
+            "POLYGON EMPTY",
+        ]
+    )
+    typ, coords, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [
+            [30.0, 10.0],
+            [40.0, 40.0],
+            [20.0, 40.0],
+            [10.0, 20.0],
+            [30.0, 10.0],
+            [35.0, 10.0],
+            [45.0, 45.0],
+            [15.0, 40.0],
+            [10.0, 20.0],
+            [35.0, 10.0],
+            [20.0, 30.0],
+            [35.0, 35.0],
+            [30.0, 20.0],
+            [20.0, 30.0],
+            [30.0, 10.0],
+            [40.0, 40.0],
+            [20.0, 40.0],
+            [10.0, 20.0],
+            [30.0, 10.0],
+        ]
+    )
+    expected_offsets1 = np.array([0, 5, 10, 14, 19])
+    expected_offsets2 = np.array([0, 1, 3, 3, 3, 4, 4])
+
+    assert typ == shapely.GeometryType.POLYGON
+    assert_allclose(coords, expected)
+    assert len(offsets) == 2
+    assert_allclose(offsets[0], expected_offsets1)
+    assert_allclose(offsets[1], expected_offsets2)
+
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    assert_geometries_equal(result, arr)
+
+
+def test_multipoints():
+    arr = shapely.from_wkt(
+        [
+            "MULTIPOINT (10 40, 40 30, 20 20, 30 10)",
+            "MULTIPOINT (30 10)",
+            "MULTIPOINT EMPTY",
+            "MULTIPOINT EMPTY",
+            "MULTIPOINT (30 10, 10 30, 40 40)",
+            "MULTIPOINT EMPTY",
+        ]
+    )
+    typ, coords, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [
+            [10.0, 40.0],
+            [40.0, 30.0],
+            [20.0, 20.0],
+            [30.0, 10.0],
+            [30.0, 10.0],
+            [30.0, 10.0],
+            [10.0, 30.0],
+            [40.0, 40.0],
+        ]
+    )
+    expected_offsets = np.array([0, 4, 5, 5, 5, 8, 8])
+
+    assert typ == shapely.GeometryType.MULTIPOINT
+    assert_allclose(coords, expected)
+    assert len(offsets) == 1
+    assert_allclose(offsets[0], expected_offsets)
+
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    assert_geometries_equal(result, arr)
+
+
+def test_multilinestrings():
+    arr = shapely.from_wkt(
+        [
+            "MULTILINESTRING ((30 10, 10 30, 40 40))",
+            "MULTILINESTRING ((10 10, 20 20, 10 40),(40 40, 30 30, 40 20, 30 10))",
+            "MULTILINESTRING EMPTY",
+            "MULTILINESTRING EMPTY",
+            "MULTILINESTRING ((35 10, 45 45), (15 40, 10 20), (30 10, 10 30, 40 40))",
+            "MULTILINESTRING EMPTY",
+        ]
+    )
+    typ, coords, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [
+            [30.0, 10.0],
+            [10.0, 30.0],
+            [40.0, 40.0],
+            [10.0, 10.0],
+            [20.0, 20.0],
+            [10.0, 40.0],
+            [40.0, 40.0],
+            [30.0, 30.0],
+            [40.0, 20.0],
+            [30.0, 10.0],
+            [35.0, 10.0],
+            [45.0, 45.0],
+            [15.0, 40.0],
+            [10.0, 20.0],
+            [30.0, 10.0],
+            [10.0, 30.0],
+            [40.0, 40.0],
+        ]
+    )
+    expected_offsets1 = np.array([0, 3, 6, 10, 12, 14, 17])
+    expected_offsets2 = np.array([0, 1, 3, 3, 3, 6, 6])
+
+    assert typ == shapely.GeometryType.MULTILINESTRING
+    assert_allclose(coords, expected)
+    assert len(offsets) == 2
+    assert_allclose(offsets[0], expected_offsets1)
+    assert_allclose(offsets[1], expected_offsets2)
+
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    assert_geometries_equal(result, arr)
+
+
+def test_multipolygons():
+    arr = shapely.from_wkt(
+        [
+            "MULTIPOLYGON (((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30)))",
+            "MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),((20 35, 10 30, 10 10, 30 5, 45 20, 20 35),(30 20, 20 15, 20 25, 30 20)))",
+            "MULTIPOLYGON EMPTY",
+            "MULTIPOLYGON EMPTY",
+            "MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)))",
+            "MULTIPOLYGON EMPTY",
+        ]
+    )
+    typ, coords, offsets = shapely.to_ragged_array(arr)
+    expected = np.array(
+        [
+            [35.0, 10.0],
+            [45.0, 45.0],
+            [15.0, 40.0],
+            [10.0, 20.0],
+            [35.0, 10.0],
+            [20.0, 30.0],
+            [35.0, 35.0],
+            [30.0, 20.0],
+            [20.0, 30.0],
+            [40.0, 40.0],
+            [20.0, 45.0],
+            [45.0, 30.0],
+            [40.0, 40.0],
+            [20.0, 35.0],
+            [10.0, 30.0],
+            [10.0, 10.0],
+            [30.0, 5.0],
+            [45.0, 20.0],
+            [20.0, 35.0],
+            [30.0, 20.0],
+            [20.0, 15.0],
+            [20.0, 25.0],
+            [30.0, 20.0],
+            [40.0, 40.0],
+            [20.0, 45.0],
+            [45.0, 30.0],
+            [40.0, 40.0],
+        ]
+    )
+    expected_offsets1 = np.array([0, 5, 9, 13, 19, 23, 27])
+    expected_offsets2 = np.array([0, 2, 3, 5, 6])
+    expected_offsets3 = np.array([0, 1, 3, 3, 3, 4, 4])
+
+    assert typ == shapely.GeometryType.MULTIPOLYGON
+    assert_allclose(coords, expected)
+    assert len(offsets) == 3
+    assert_allclose(offsets[0], expected_offsets1)
+    assert_allclose(offsets[1], expected_offsets2)
+    assert_allclose(offsets[2], expected_offsets3)
+
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    assert_geometries_equal(result, arr)
+
+
+def test_mixture_point_multipoint():
+    typ, coords, offsets = shapely.to_ragged_array([point, multi_point])
+    assert typ == shapely.GeometryType.MULTIPOINT
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    expected = np.array([MultiPoint([point]), multi_point])
+    assert_geometries_equal(result, expected)
+
+
+def test_mixture_linestring_multilinestring():
+    typ, coords, offsets = shapely.to_ragged_array([line_string, multi_line_string])
+    assert typ == shapely.GeometryType.MULTILINESTRING
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    expected = np.array([MultiLineString([line_string]), multi_line_string])
+    assert_geometries_equal(result, expected)
+
+
+def test_mixture_polygon_multipolygon():
+    typ, coords, offsets = shapely.to_ragged_array([polygon, multi_polygon])
+    assert typ == shapely.GeometryType.MULTIPOLYGON
+    result = shapely.from_ragged_array(typ, coords, offsets)
+    expected = np.array([MultiPolygon([polygon]), multi_polygon])
+    assert_geometries_equal(result, expected)