Repository: geopandas/dask-geopandas

Issue:
Handling of empty geometries in Hilbert/Morton distance methods
Currently the functionality "works" (i.e. doesn't break down) when empty geometries are present:

```
>>>s = geopandas.GeoSeries.from_wkt(["POINT (1 2)", "POINT (2 1)", "POLYGON EMPTY", "POINT (10 10)"])
>>> ds = dask_geopandas.from_geopandas(s, npartitions=1)
>>> ds.hilbert_distance().compute()
0     16691136
1      5509098
2            0
3    715827882
Name: hilbert_distance, dtype: uint32

```

but so empty geometries get filled in with 0. 

Is this OK? Or do we want to raise an error instead if there are empty geometries?

When using this in a context of spatial shuffling, this would mean that all empty ones end up in the same partition (which seems OK if for some reason you want to keep those empty geometries)

Discussion:

Patch:
diff --git a/dask_geopandas/__init__.py b/dask_geopandas/__init__.py
index 8cd6ebe..a776641 100644
--- a/dask_geopandas/__init__.py
+++ b/dask_geopandas/__init__.py
@@ -11,6 +11,7 @@ from .core import (
 )
 from .io.file import read_file
 from .io.parquet import read_parquet, to_parquet
+from .io.arrow import read_feather, to_feather
 from .sjoin import sjoin
 
 
@@ -23,7 +24,9 @@ __all__ = [
     "GeoSeries",
     "from_geopandas",
     "from_dask_dataframe",
+    "read_feather",
     "read_parquet",
+    "to_feather",
     "to_parquet",
     "clip",
     "sjoin",
diff --git a/dask_geopandas/core.py b/dask_geopandas/core.py
index 9843f31..1e261b5 100644
--- a/dask_geopandas/core.py
+++ b/dask_geopandas/core.py
@@ -519,6 +519,11 @@ class GeoDataFrame(_Frame, dd.core.DataFrame):
             meta = self._meta.set_geometry(col)
         return self.map_partitions(M.set_geometry, col, meta=meta)
 
+    @derived_from(geopandas.GeoDataFrame)
+    def rename_geometry(self, col):
+        meta = self._meta.rename_geometry(col)
+        return self.map_partitions(M.rename_geometry, col, meta=meta)
+
     def __getitem__(self, key):
         """
         If the result is a new dask_geopandas.GeoDataFrame/GeoSeries (automatically
@@ -543,6 +548,12 @@ class GeoDataFrame(_Frame, dd.core.DataFrame):
 
         return to_parquet(self, path, *args, **kwargs)
 
+    def to_feather(self, path, *args, **kwargs):
+        """See dask_geopadandas.to_feather docstring for more information"""
+        from .io.arrow import to_feather
+
+        return to_feather(self, path, *args, **kwargs)
+
     def dissolve(self, by=None, aggfunc="first", split_out=1, **kwargs):
         """Dissolve geometries within ``groupby`` into a single geometry.
 
diff --git a/dask_geopandas/geohash.py b/dask_geopandas/geohash.py
index 407e157..1ec95e6 100644
--- a/dask_geopandas/geohash.py
+++ b/dask_geopandas/geohash.py
@@ -8,6 +8,8 @@ The vectorized implementation for quantization and bit interleaving is in turn b
 "Geohash in Golang Assembly" blog (https://mmcloughlin.com/posts/geohash-assembly).
 
 """
+import warnings
+
 import numpy as np
 import pandas as pd
 
@@ -31,6 +33,15 @@ def _geohash(gdf, as_string, precision):
     type : pandas.Series
         Series containing geohash
     """
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", "GeoSeries.isna() previously returned True", UserWarning
+        )
+        if gdf.is_empty.any() | gdf.geometry.isna().any():
+            raise ValueError(
+                "Geohash cannot be computed on a GeoSeries with empty or "
+                "missing geometries.",
+            )
 
     # Calculate bounds
     bounds = gdf.bounds.to_numpy()
diff --git a/dask_geopandas/hilbert_distance.py b/dask_geopandas/hilbert_distance.py
index 1b7b9e8..a05c60a 100644
--- a/dask_geopandas/hilbert_distance.py
+++ b/dask_geopandas/hilbert_distance.py
@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pandas as pd
 
@@ -24,6 +26,15 @@ def _hilbert_distance(gdf, total_bounds=None, level=16):
     Pandas Series containing distances along the Hilbert curve
 
     """
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", "GeoSeries.isna() previously returned True", UserWarning
+        )
+        if gdf.is_empty.any() | gdf.geometry.isna().any():
+            raise ValueError(
+                "Hilbert distance cannot be computed on a GeoSeries with empty or "
+                "missing geometries.",
+            )
     # Calculate bounds as numpy array
     bounds = gdf.bounds.to_numpy()
 
diff --git a/dask_geopandas/io/arrow.py b/dask_geopandas/io/arrow.py
new file mode 100644
index 0000000..cfdc13c
--- /dev/null
+++ b/dask_geopandas/io/arrow.py
@@ -0,0 +1,414 @@
+import copy
+import json
+import math
+from typing import TYPE_CHECKING
+
+from dask.base import compute_as_if_collection, tokenize
+from dask.dataframe.core import new_dd_object, Scalar
+from dask.highlevelgraph import HighLevelGraph
+from dask.layers import DataFrameIOLayer
+from dask.utils import natural_sort_key, apply
+
+import pandas as pd
+import geopandas
+import shapely.geometry
+
+from fsspec.core import get_fs_token_paths
+
+
+if TYPE_CHECKING:
+    import pyarrow
+
+
+def _update_meta_to_geodataframe(meta, schema_metadata):
+    """
+    Convert meta to a GeoDataFrame and update with potential GEO metadata
+    """
+    if schema_metadata and b"geo" in schema_metadata:
+        geo_meta = json.loads(schema_metadata[b"geo"])
+        geometry_column_name = geo_meta["primary_column"]
+        crs = geo_meta["columns"][geometry_column_name]["crs"]
+        geometry_columns = geo_meta["columns"]
+    else:
+        # TODO we could allow the user to pass those explicitly if not
+        # stored in the metadata
+        raise ValueError(
+            "Missing geo metadata in the Parquet/Feather file. "
+            "Use dask.dataframe.read_parquet/pandas.read_feather() instead."
+        )
+
+    # Update meta to be a GeoDataFrame
+    meta = geopandas.GeoDataFrame(meta, geometry=geometry_column_name, crs=crs)
+    for col, item in geometry_columns.items():
+        if not col == meta._geometry_column_name:
+            meta[col] = geopandas.GeoSeries(meta[col], crs=item["crs"])
+
+    return meta
+
+
+def _get_partition_bounds(schema_metadata):
+    """
+    Get the partition bounds, if available, for the dataset fragment.
+    """
+    if not (schema_metadata and b"geo" in schema_metadata):
+        return None
+
+    metadata = json.loads(schema_metadata[b"geo"].decode("utf-8"))
+
+    # for now only check the primary column (TODO generalize this to follow
+    # the logic of geopandas to fallback to other geometry columns)
+    geometry = metadata["primary_column"]
+    bbox = metadata["columns"][geometry].get("bbox", None)
+    if bbox is None or all(math.isnan(val) for val in bbox):
+        return None
+    return shapely.geometry.box(*bbox)
+
+
+class ArrowDatasetEngine:
+    """
+    Custom IO engine based on pyarrow.dataset.
+
+    This is designed after dask's ArrowDatasetEngine for Parquet IO (but simpler
+    with less options, and not dealing with a legacy engine) and ArrowORCEngine
+    for ORC IO (but using pyarrow.dataset for the read_metadata discovery).
+    """
+
+    file_format: str
+
+    @classmethod
+    def read_metadata(cls, fs, paths, columns, filters, index):
+
+        import pyarrow.dataset as ds
+        from pyarrow.parquet import _filters_to_expression
+
+        # dataset discovery
+        if len(paths) == 1:
+            # list of 1 directory path is not supported
+            paths = paths[0]
+        dataset = ds.dataset(
+            paths, partitioning="hive", filesystem=fs, format=cls.file_format
+        )
+
+        # Get all (filtered) fragments
+        if filters is not None:
+            filter = _filters_to_expression(filters)
+        else:
+            filter = None
+
+        fragments = list(dataset.get_fragments(filter=filter))
+
+        # numeric rather than glob ordering
+        # TODO how does this handle different partitioned directories?
+        fragments = sorted(fragments, key=lambda f: natural_sort_key(f.path))
+
+        # TODO potential splitting / aggregating of fragments
+
+        # Create dask meta
+        schema = dataset.schema
+        # TODO add support for `categories`keyword
+        meta = schema.empty_table().to_pandas()
+
+        if index:
+            meta = meta.set_index(index)
+
+        if columns is not None:
+            ex = set(columns) - set(meta.columns)
+            if ex:
+                raise ValueError(
+                    f"Requested columns {ex} not in schema {set(meta.columns)}"
+                )
+            meta = meta[columns]
+
+        return fragments, meta, schema, filter
+
+    @classmethod
+    def _arrow_table_to_pandas(
+        cls, arrow_table: "pyarrow.Table", categories, **kwargs
+    ) -> pd.DataFrame:
+        _kwargs = kwargs.get("arrow_to_pandas", {})
+        _kwargs.update({"use_threads": False, "ignore_metadata": False})
+
+        return arrow_table.to_pandas(categories=categories, **_kwargs)
+
+    @classmethod
+    def read_partition(cls, fs, fragment, schema, columns, filter, **kwargs):
+
+        table = fragment.to_table(
+            schema=schema, columns=columns, filter=filter, use_threads=False
+        )
+        df = cls._arrow_table_to_pandas(table, None)
+        return df
+
+    @classmethod
+    def write_partition(cls, df, path, fs, filename, **kwargs):
+        from pyarrow import feather
+
+        table = cls._pandas_to_arrow_table(df, preserve_index=None)
+        # TODO using the datasets API could automatically support partitioning
+        # on columns
+        with fs.open(fs.sep.join([path, filename]), "wb") as f:
+            feather.write_feather(table, f)
+
+
+class GeoDatasetEngine:
+    """
+    Mixin to combine with an IO Engine (the custom engine defined above for
+    Feather IO, or dask's engine for Parquet IO) that holds the custom logic
+    for geospatial data: overriding the arrow <-> pandas conversions to ensure
+    we read/write GeoDataFrames.
+
+    """
+
+    @classmethod
+    def _arrow_table_to_pandas(
+        cls, arrow_table: "pyarrow.Table", categories, **kwargs
+    ) -> pd.DataFrame:
+        from geopandas.io.arrow import _arrow_to_geopandas
+
+        _kwargs = kwargs.get("arrow_to_pandas", {})
+        _kwargs.update({"use_threads": False, "ignore_metadata": False})
+
+        # TODO support additional keywords
+        try:
+            return _arrow_to_geopandas(arrow_table)
+        except ValueError as err:
+            # when no geometry column is selected, the above will error.
+            # We want to fallback to reading it as a plain dask object, because
+            # the column selection can be an automatic pushdown (eg `ddf['col']`)
+            # TODO more robust detection of when to fall back?
+            if "No geometry columns are included" in str(err):
+                return super()._arrow_table_to_pandas(
+                    arrow_table, categories=categories, **kwargs
+                )
+            else:
+                raise
+
+    @classmethod
+    def _pandas_to_arrow_table(
+        cls, df: pd.DataFrame, preserve_index=False, schema=None
+    ) -> "pyarrow.Table":
+        from geopandas.io.arrow import _geopandas_to_arrow
+
+        # TODO add support for schema
+        if schema is not None:
+            raise NotImplementedError("Passing 'schema' is not yet supported")
+
+        table = _geopandas_to_arrow(df, index=preserve_index)
+        return table
+
+
+class FeatherDatasetEngine(GeoDatasetEngine, ArrowDatasetEngine):
+    file_format = "feather"
+
+
+class FeatherFunctionWrapper:
+    """
+    Feather Function-Wrapper Class
+    Reads Feather data from disk to produce a partition.
+    """
+
+    def __init__(self, engine, fs, columns, filter, schema, index):
+        self.engine = engine
+        self.fs = fs
+        self.columns = columns
+        self.filter = filter
+        self.schema = schema
+        self.index = index
+
+    def project_columns(self, columns):
+        """Return a new FeatherFunctionWrapper object with
+        a sub-column projection.
+        """
+        if columns == self.columns:
+            return self
+        func = copy.deepcopy(self)
+        func.columns = columns
+        return func
+
+    def __call__(self, parts):
+        _df = self.engine.read_partition(
+            self.fs, parts, self.schema, self.columns, self.filter
+        )
+        if self.index:
+            _df.set_index(self.index, inplace=True)
+        return _df
+
+
+def read_feather(
+    path,
+    columns=None,
+    filters=None,
+    index=None,
+    storage_options=None,
+):
+    """Read a Feather dataset into a Dask-GeoPandas DataFrame.
+
+    Parameters
+    ----------
+    path: str or list(str)
+        Source directory for data, or path(s) to individual Feather files.
+        Paths can be a full URL with protocol specifier, and may include
+        glob character if a single string.
+    columns: None or list(str)
+        Columns to load. If None, loads all.
+    filters : list (of list) of tuples or pyarrow.dataset.Expression, default None
+        Row-wise filter to apply while reading the dataset. Can be specified
+        as a ``pyarrow.dataset.Expression`` object or using a list of tuples
+        notation, like ``[[('col1', '==', 0), ...], ...]``. The filter will
+        be applied both at the partition level, this is to prevent the loading
+        of some files, as at the file level to filter the actual rows.
+
+        For the list of tuples format, predicates can be expressed in disjunctive
+        normal form (DNF). This means that the innermost tuple describes a single
+        column predicate. These inner predicates are combined with an AND
+        conjunction into a larger predicate. The outer-most list then combines all
+        of the combined filters with an OR disjunction.
+
+        Predicates can also be expressed as a List[Tuple]. These are evaluated
+        as an AND conjunction. To express OR in predictates, one must use the
+        List[List[Tuple]] notation.
+    index : str, list or False, default None
+        Field name(s) to use as the output frame index. By default will be
+        inferred from the pandas metadata (if present in the files). Use False
+        to read all fields as columns.
+    storage_options : dict, default None
+        Key/value pairs to be passed on to the file-system backend, if any
+        (inferred from the path, such as "s3://...").
+        Please see ``fsspec`` for more details.
+
+    Returns
+    -------
+    dask_geopandas.GeoDataFrame (even if there is only one column)
+
+    """
+    if index is False:
+        raise NotImplementedError("Specifying index=False is not yet implemented")
+
+    # Get engine
+    engine = FeatherDatasetEngine
+
+    # Process file path(s)
+    storage_options = storage_options or {}
+    fs, _, paths = get_fs_token_paths(path, mode="rb", storage_options=storage_options)
+    paths = sorted(paths, key=natural_sort_key)  # numeric rather than glob ordering
+
+    # Let backend engine generate a list of parts from the dataset metadata
+    parts, meta, schema, filter = engine.read_metadata(
+        fs,
+        paths,
+        columns,
+        filters,
+        index,
+    )
+
+    # Update meta to be a GeoDataFrame
+    meta = _update_meta_to_geodataframe(meta, schema.metadata)
+
+    # Construct spatial partitioning information, if available
+    spatial_partitions = geopandas.GeoSeries(
+        [_get_partition_bounds(frag.physical_schema.metadata) for frag in parts],
+        crs=meta.crs,
+    )
+    if spatial_partitions.isna().any():
+        spatial_partitions = None
+
+    # Construct and return a Blockwise layer
+    label = "read-feather-"
+    output_name = label + tokenize(path, columns, filters, index)
+    layer = DataFrameIOLayer(
+        output_name,
+        columns,
+        parts,
+        FeatherFunctionWrapper(engine, fs, columns, filter, schema, index),
+        label=label,
+    )
+    graph = HighLevelGraph({output_name: layer}, {output_name: set()})
+    result = new_dd_object(graph, output_name, meta, [None] * (len(parts) + 1))
+    result.spatial_partitions = spatial_partitions
+    return result
+
+
+def to_feather(
+    df,
+    path,
+    write_index=True,
+    storage_options=None,
+    compute=True,
+    compute_kwargs=None,
+):
+    """Store Dask.dataframe to Feather files
+
+    Notes
+    -----
+    Each partition will be written to a separate file.
+
+    Parameters
+    ----------
+    df : dask_geopandas.GeoDataFrame
+    path : string or pathlib.Path
+        Destination directory for data.  Prepend with protocol like ``s3://``
+        or ``hdfs://`` for remote data.
+    write_index : boolean, default True
+        Whether or not to write the index. Defaults to True.
+    storage_options : dict, default None
+        Key/value pairs to be passed on to the file-system backend, if any
+        (inferred from the path, such as "s3://...").
+        Please see ``fsspec`` for more details.
+    compute : bool, default True
+        If True (default) then the result is computed immediately. If False
+        then a ``dask.delayed`` object is returned for future computation.
+    compute_kwargs : dict, default True
+        Options to be passed in to the compute method
+
+    See Also
+    --------
+    dask_geopandas.read_feather: Read Feather data to dask.dataframe
+    """
+    # based on the to_orc function from dask
+
+    # Get engine
+    engine = FeatherDatasetEngine
+
+    # Process file path
+    storage_options = storage_options or {}
+    fs, _, _ = get_fs_token_paths(path, mode="wb", storage_options=storage_options)
+    # Trim any protocol information from the path before forwarding
+    path = fs._strip_protocol(path)
+
+    if not write_index:
+        # Not writing index - might as well drop it
+        df = df.reset_index(drop=True)
+
+    # Use df.npartitions to define file-name list
+    fs.mkdirs(path, exist_ok=True)
+    filenames = [f"part.{i}.feather" for i in range(df.npartitions)]
+
+    # Construct IO graph
+    dsk = {}
+    name = "to-feather-" + tokenize(df, fs, path, write_index, storage_options)
+    part_tasks = []
+    for d, filename in enumerate(filenames):
+        dsk[(name, d)] = (
+            apply,
+            engine.write_partition,
+            [
+                (df._name, d),
+                path,
+                fs,
+                filename,
+            ],
+        )
+        part_tasks.append((name, d))
+    dsk[name] = (lambda x: None, part_tasks)
+    graph = HighLevelGraph.from_collections(name, dsk, dependencies=[df])
+
+    # Compute or return future
+    if compute:
+        if compute_kwargs is None:
+            compute_kwargs = dict()
+        from dask_geopandas import GeoDataFrame
+
+        return compute_as_if_collection(
+            GeoDataFrame, graph, part_tasks, **compute_kwargs
+        )
+    return Scalar(graph, name, "")
diff --git a/dask_geopandas/io/parquet.py b/dask_geopandas/io/parquet.py
index a2315a5..a4eb1d0 100644
--- a/dask_geopandas/io/parquet.py
+++ b/dask_geopandas/io/parquet.py
@@ -1,26 +1,19 @@
 from functools import partial
-import json
-import math
-from typing import TYPE_CHECKING
-
-import pandas as pd
 
 import geopandas
-import shapely.geometry
 
 import dask.dataframe as dd
 
+from .arrow import GeoDatasetEngine, _get_partition_bounds, _update_meta_to_geodataframe
+
 try:
     # pyarrow is imported here, but is an optional dependency
     from dask.dataframe.io.parquet.arrow import ArrowEngine
 except ImportError:
     ArrowEngine = object
 
-if TYPE_CHECKING:
-    import pyarrow
 
-
-def _get_partition_bounds(part, fs):
+def _get_partition_bounds_parquet(part, fs):
     """
     Based on the part information gathered by dask, get the partition bounds
     if available.
@@ -40,29 +33,24 @@ def _get_partition_bounds(part, fs):
     if pq_metadata is None:
         return None
 
-    metadata_str = pq_metadata.metadata.get(b"geo", None)
-    if metadata_str is None:
-        return None
+    return _get_partition_bounds(pq_metadata.metadata)
 
-    metadata = json.loads(metadata_str.decode("utf-8"))
 
-    # for now only check the primary column (TODO generalize this to follow
-    # the logic of geopandas to fallback to other geometry columns)
-    geometry = metadata["primary_column"]
-    bbox = metadata["columns"][geometry].get("bbox", None)
-    if bbox is None or all(math.isnan(val) for val in bbox):
-        return None
-    return shapely.geometry.box(*bbox)
+class GeoArrowEngine(GeoDatasetEngine, ArrowEngine):
+    """
+    Engine for reading geospatial Parquet datasets. Subclasses dask's
+    ArrowEngine for Parquet, but overriding some methods to ensure we
+    correctly read/write GeoDataFrames.
 
+    """
 
-class GeoArrowEngine(ArrowEngine):
     @classmethod
     def read_metadata(cls, fs, paths, **kwargs):
         meta, stats, parts, index = super().read_metadata(fs, paths, **kwargs)
 
         # get spatial partitions if available
         regions = geopandas.GeoSeries(
-            [_get_partition_bounds(part, fs) for part in parts], crs=meta.crs
+            [_get_partition_bounds_parquet(part, fs) for part in parts], crs=meta.crs
         )
         if regions.notna().all():
             # a bit hacky, but this allows us to get this passed through
@@ -75,25 +63,7 @@ class GeoArrowEngine(ArrowEngine):
         """
         Convert meta to a GeoDataFrame and update with potential GEO metadata
         """
-        if schema.metadata and b"geo" in schema.metadata:
-            geo_meta = json.loads(schema.metadata[b"geo"])
-            geometry_column_name = geo_meta["primary_column"]
-            crs = geo_meta["columns"][geometry_column_name]["crs"]
-            geometry_columns = geo_meta["columns"]
-        else:
-            # TODO we could allow the user to pass those explicitly if not
-            # stored in the metadata
-            geometry_column_name = None
-            crs = None
-            geometry_columns = {}
-
-        # Update meta to be a GeoDataFrame
-        meta = geopandas.GeoDataFrame(meta, geometry=geometry_column_name, crs=crs)
-        for col, item in geometry_columns.items():
-            if not col == meta._geometry_column_name:
-                meta[col] = geopandas.GeoSeries(meta[col], crs=item["crs"])
-
-        return meta
+        return _update_meta_to_geodataframe(meta, schema.metadata)
 
     @classmethod
     def _generate_dd_meta(cls, schema, index, categories, partition_info):
@@ -112,43 +82,6 @@ class GeoArrowEngine(ArrowEngine):
         meta = cls._update_meta(meta, schema)
         return meta
 
-    @classmethod
-    def _arrow_table_to_pandas(
-        cls, arrow_table: "pyarrow.Table", categories, **kwargs
-    ) -> pd.DataFrame:
-        from geopandas.io.arrow import _arrow_to_geopandas
-
-        _kwargs = kwargs.get("arrow_to_pandas", {})
-        _kwargs.update({"use_threads": False, "ignore_metadata": False})
-
-        # TODO support additional keywords
-        try:
-            return _arrow_to_geopandas(arrow_table)
-        except ValueError as err:
-            # when no geometry column is selected, the above will error.
-            # We want to fallback to reading it as a plain dask object, because
-            # the column selection can be an automatic pushdown (eg `ddf['col']`)
-            # TODO more robust detection of when to fall back?
-            if "No geometry columns are included" in str(err):
-                return super()._arrow_table_to_pandas(
-                    arrow_table, categories=categories, **kwargs
-                )
-            else:
-                raise
-
-    @classmethod
-    def _pandas_to_arrow_table(
-        cls, df: pd.DataFrame, preserve_index=False, schema=None
-    ) -> "pyarrow.Table":
-        from geopandas.io.arrow import _geopandas_to_arrow
-
-        # TODO add support for schema
-        if schema is not None:
-            raise NotImplementedError("Passing 'schema' is not yet supported")
-
-        table = _geopandas_to_arrow(df, index=preserve_index)
-        return table
-
 
 to_parquet = partial(dd.to_parquet, engine=GeoArrowEngine)
 to_parquet.__doc__ = dd.to_parquet.__doc__
diff --git a/dask_geopandas/morton_distance.py b/dask_geopandas/morton_distance.py
index 1272bd8..e0fbb4e 100644
--- a/dask_geopandas/morton_distance.py
+++ b/dask_geopandas/morton_distance.py
@@ -1,3 +1,5 @@
+import warnings
+
 import pandas as pd
 from dask_geopandas.hilbert_distance import _continuous_to_discrete_coords
 
@@ -22,6 +24,15 @@ def _morton_distance(gdf, total_bounds, level):
         Series containing distances from Morton curve
 
     """
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", "GeoSeries.isna() previously returned True", UserWarning
+        )
+        if gdf.is_empty.any() | gdf.geometry.isna().any():
+            raise ValueError(
+                "Morton distance cannot be computed on a GeoSeries with empty or "
+                "missing geometries.",
+            )
     # Calculate bounds as numpy array
     bounds = gdf.bounds.to_numpy()
     # Calculate discrete coords based on total bounds and bounds
diff --git a/doc/source/docs/reference/geodataframe.rst b/doc/source/docs/reference/geodataframe.rst
index 0021c23..0e69dcc 100644
--- a/doc/source/docs/reference/geodataframe.rst
+++ b/doc/source/docs/reference/geodataframe.rst
@@ -38,6 +38,7 @@ Active geometry handling
    :toctree: api/
 
    GeoDataFrame.set_geometry
+   GeoDataFrame.rename_geometry
 
 Aggregating and exploding
 -------------------------


Test patch:
diff --git a/tests/io/test_arrow.py b/tests/io/test_arrow.py
new file mode 100644
index 0000000..3ae1367
--- /dev/null
+++ b/tests/io/test_arrow.py
@@ -0,0 +1,215 @@
+import geopandas
+import dask_geopandas
+import dask.dataframe as dd
+
+import pytest
+from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
+from pandas.testing import assert_index_equal
+
+
+pa = pytest.importorskip("pyarrow")
+ds = pytest.importorskip("pyarrow.dataset")
+
+
+pytestmark = pytest.mark.filterwarnings(
+    "ignore:this is an initial implementation:UserWarning"
+)
+
+
+def test_read(tmp_path):
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+
+    # writing a partitioned dataset with geopandas (to not rely on roundtrip)
+    basedir = tmp_path / "dataset"
+    basedir.mkdir()
+    df.iloc[:100].to_feather(basedir / "data.0.feather")
+    df.iloc[100:].to_feather(basedir / "data.1.feather")
+
+    result = dask_geopandas.read_feather(basedir)
+    assert isinstance(result, dask_geopandas.GeoDataFrame)
+    assert result.npartitions == 2
+    assert result.crs == df.crs
+    assert result.spatial_partitions is not None
+    # TODO this reset_index should not be necessary
+    result_gpd = result.compute().reset_index(drop=True)
+    assert_geodataframe_equal(result_gpd, df)
+
+
+def test_write(tmp_path):
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+
+    basedir = tmp_path / "dataset"
+    ddf.to_feather(basedir)
+
+    # each partition (4) is written as a feather file
+    paths = list(basedir.glob("*.feather"))
+    assert len(paths) == 4
+
+    # each individual file is a valid feather file
+    result_part0 = geopandas.read_feather(basedir / "part.0.feather")
+    result_part0.index.name = None
+    assert_geodataframe_equal(result_part0, df.iloc[:45])
+
+    # TODO geopandas doesn't actually support this for "feather" format
+    # # the written dataset is also readable by plain geopandas
+    # result_gpd = geopandas.read_feather(basedir)
+    # # the dataset written by dask has "__null_dask_index__" index column name
+    # result_gpd.index.name = None
+    # assert_geodataframe_equal(result_gpd, df)
+
+
+@pytest.mark.xfail  # https://github.com/dask/dask/issues/8022
+def test_write_delayed(tmp_path):
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+
+    basedir = tmp_path / "dataset"
+    dataset = ddf.to_feather(basedir, compute=False)
+    dataset.compute()
+    result = dask_geopandas.read_feather(basedir)
+    assert result.npartitions == 4
+    # TODO this reset_index should not be necessary
+    result_gpd = result.compute().reset_index(drop=True)
+    assert_geodataframe_equal(result_gpd, df)
+
+
+def test_roundtrip(tmp_path):
+    # basic roundtrip
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+
+    basedir = tmp_path / "dataset"
+    ddf.to_feather(basedir)
+
+    # reading back gives identical GeoDataFrame
+    result = dask_geopandas.read_feather(basedir)
+    assert result.npartitions == 4
+    assert result.crs == df.crs
+    # TODO this reset_index should not be necessary
+    result_gpd = result.compute().reset_index(drop=True)
+    assert_geodataframe_equal(result_gpd, df)
+    # reading back also populates the spatial partitioning property
+    ddf.calculate_spatial_partitions()
+    assert_geoseries_equal(result.spatial_partitions, ddf.spatial_partitions.envelope)
+
+
+def test_roundtrip_s3(s3_resource, s3_storage_options):
+    fs, endpoint_url = s3_resource
+
+    # basic roundtrip to S3
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+
+    uri = "s3://geopandas-test/dataset.feather"
+    ddf.to_feather(uri, storage_options=s3_storage_options)
+
+    # reading back gives identical GeoDataFrame
+    result = dask_geopandas.read_feather(uri, storage_options=s3_storage_options)
+    assert result.npartitions == 4
+    assert_geodataframe_equal(result.compute().reset_index(drop=True), df)
+    # reading back correctly sets the CRS in meta
+    assert result.crs == df.crs
+    # reading back also populates the spatial partitioning property
+    assert result.spatial_partitions is not None
+
+
+def test_column_selection_push_down(tmp_path):
+    # set up dataset
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+    basedir = tmp_path / "dataset"
+    # TODO awaiting a `to_feather` implementation
+    # ddf.to_feather(basedir)
+    basedir.mkdir()
+    for i, part in enumerate(ddf.partitions):
+        part.compute().to_feather(basedir / f"part.{i}.feather")
+
+    ddf = dask_geopandas.read_feather(basedir)
+
+    # selecting columns including geometry column still gives GeoDataFrame
+    ddf_subset = ddf[["pop_est", "geometry"]]
+    assert type(ddf_subset) is dask_geopandas.GeoDataFrame
+    # and also preserves the spatial partitioning information
+    assert ddf_subset.spatial_partitions is not None
+
+    # selecting a single non-geometry column on the dataframe should work
+    s = ddf["pop_est"]
+    assert type(s) is dd.Series
+    assert s.max().compute() == df["pop_est"].max()
+
+
+def test_missing_metadata(tmp_path):
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    path = tmp_path / "test.feather"
+
+    # convert to DataFrame with wkb -> writing to feather will have only pandas metadata
+    df = df.to_wkb()
+    df.to_feather(path)
+
+    with pytest.raises(ValueError, match="Missing geo metadata"):
+        dask_geopandas.read_feather(path)
+
+    # remove metadata completely
+    from pyarrow import feather
+
+    table = feather.read_table(path)
+    feather.write_feather(table.replace_schema_metadata(), path)
+
+    with pytest.raises(ValueError, match="Missing geo metadata"):
+        dask_geopandas.read_feather(path)
+
+
+@pytest.mark.parametrize(
+    "filter", [[("continent", "=", "Africa")], ds.field("continent") == "Africa"]
+)
+def test_filters(tmp_path, filter):
+    # set up dataset
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+    basedir = tmp_path / "dataset"
+    ddf.to_feather(basedir)
+
+    # specifying filters argument
+    result = dask_geopandas.read_feather(basedir, filters=filter)
+    assert result.npartitions == 4
+
+    result_gpd = result.compute().reset_index(drop=True)
+    expected = df[df["continent"] == "Africa"].reset_index(drop=True)
+    assert_geodataframe_equal(result_gpd, expected)
+
+
+def test_index(tmp_path):
+    # set up dataset
+    df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+    # get meaningful index by shuffling (hilbert distance)
+    df = dask_geopandas.from_geopandas(df, npartitions=1).spatial_shuffle().compute()
+    ddf = dask_geopandas.from_geopandas(df, npartitions=4)
+
+    # roundtrip preserves the index by default
+    basedir = tmp_path / "dataset"
+    ddf.to_feather(basedir)
+    result = dask_geopandas.read_feather(basedir)
+    assert "hilbert_distance" not in result.columns
+    assert result.index.name == "hilbert_distance"
+    assert_index_equal(result.index.compute(), df.index)
+
+    # TODO not setting the index
+    with pytest.raises(NotImplementedError):
+        result = dask_geopandas.read_feather(basedir, index=False)
+    # assert "hilbert_distance" in result.columns
+    # assert result.index.name is None
+
+    # setting specific columns as the index
+    result = dask_geopandas.read_feather(basedir, index="iso_a3")
+    assert "iso_a3" not in result.columns
+    assert result.index.name == "iso_a3"
+    assert_geodataframe_equal(result.compute(), df.set_index("iso_a3"))
+
+    # not writing the index
+    basedir = tmp_path / "dataset"
+    ddf.to_feather(basedir, write_index=False)
+    result = dask_geopandas.read_feather(basedir)
+    assert "hilbert_distance" not in result.columns
+    assert result.index.name is None
+    assert result.index.compute()[0] == 0
diff --git a/tests/test_core.py b/tests/test_core.py
index 8571e7e..29886c6 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -363,6 +363,26 @@ def test_set_geometry_with_dask_geoseries():
     assert_geoseries_equal(dask_obj.geometry.compute(), expected.geometry)
 
 
+def test_rename_geometry(geodf_points):
+    df = geodf_points
+    dask_obj = dask_geopandas.from_geopandas(df, npartitions=2)
+    renamed = dask_obj.rename_geometry("points")
+    assert renamed._meta.geometry.name == "points"
+
+    for part in renamed.partitions:
+        assert part.compute().geometry.name == "points"
+
+    result = renamed.compute()
+    assert_geodataframe_equal(result, df.rename_geometry("points"))
+
+
+def test_rename_geometry_error(geodf_points):
+    df = geodf_points
+    dask_obj = dask_geopandas.from_geopandas(df, npartitions=2)
+    with pytest.raises(ValueError, match="Column named value1 already exists"):
+        dask_obj.rename_geometry("value1")
+
+
 def test_from_dask_dataframe_with_dask_geoseries():
     df = pd.DataFrame({"x": [0, 1, 2, 3], "y": [1, 2, 3, 4]})
     dask_obj = dd.from_pandas(df, npartitions=2)
diff --git a/tests/test_geohash.py b/tests/test_geohash.py
index f43f585..2533e77 100644
--- a/tests/test_geohash.py
+++ b/tests/test_geohash.py
@@ -8,6 +8,7 @@ from dask_geopandas.geohash import _calculate_mid_points
 from dask_geopandas import from_geopandas
 import geopandas
 from shapely.geometry import Point, LineString, Polygon
+from shapely.wkt import loads
 
 
 @pytest.fixture
@@ -80,3 +81,20 @@ def test_world():
     geohash_dask(
         geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres")).iloc[1:]
     )
+
+
+@pytest.mark.parametrize(
+    "empty",
+    [
+        None,
+        loads("POLYGON EMPTY"),
+    ],
+)
+def test_empty(geoseries_polygons, empty):
+    s = geoseries_polygons
+    s.iloc[-1] = empty
+    dask_obj = from_geopandas(s, npartitions=2)
+    with pytest.raises(
+        ValueError, match="cannot be computed on a GeoSeries with empty"
+    ):
+        dask_obj.geohash().compute()
diff --git a/tests/test_hilbert_distance.py b/tests/test_hilbert_distance.py
index 4240781..3e8750f 100644
--- a/tests/test_hilbert_distance.py
+++ b/tests/test_hilbert_distance.py
@@ -11,6 +11,7 @@ from dask_geopandas.hilbert_distance import (
 from dask_geopandas import from_geopandas
 import geopandas
 from shapely.geometry import Point, LineString, Polygon
+from shapely.wkt import loads
 
 
 def test_hilbert_distance():
@@ -112,3 +113,20 @@ def test_world():
     hilbert_distance_dask(
         geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres")).iloc[1:]
     )
+
+
+@pytest.mark.parametrize(
+    "empty",
+    [
+        None,
+        loads("POLYGON EMPTY"),
+    ],
+)
+def test_empty(geoseries_polygons, empty):
+    s = geoseries_polygons
+    s.iloc[-1] = empty
+    dask_obj = from_geopandas(s, npartitions=2)
+    with pytest.raises(
+        ValueError, match="cannot be computed on a GeoSeries with empty"
+    ):
+        dask_obj.hilbert_distance().compute()
diff --git a/tests/test_morton_distance.py b/tests/test_morton_distance.py
index 1df7e93..a6ebf45 100644
--- a/tests/test_morton_distance.py
+++ b/tests/test_morton_distance.py
@@ -6,6 +6,7 @@ from dask_geopandas.hilbert_distance import _continuous_to_discrete_coords
 from dask_geopandas import from_geopandas
 import geopandas
 from shapely.geometry import Point, LineString, Polygon
+from shapely.wkt import loads
 
 
 @pytest.fixture
@@ -81,3 +82,20 @@ def test_world():
     morton_distance_dask(
         geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres")).iloc[1:]
     )
+
+
+@pytest.mark.parametrize(
+    "empty",
+    [
+        None,
+        loads("POLYGON EMPTY"),
+    ],
+)
+def test_empty(geoseries_polygons, empty):
+    s = geoseries_polygons
+    s.iloc[-1] = empty
+    dask_obj = from_geopandas(s, npartitions=2)
+    with pytest.raises(
+        ValueError, match="cannot be computed on a GeoSeries with empty"
+    ):
+        dask_obj.morton_distance().compute()