Repository: nmandery/h3ronpy (language: python)

Issue:
Rebuild ontop h3o and arrow2
Use arrow2 to exchange arrays with python. Migrate to h3o. Make use of [h3arrow](https://github.com/nmandery/h3arrow) and [rasterh3](https://github.com/nmandery/rasterh3)

By building ontop arrow, we can also provide direct support for `polars` without requiring `pandas` for most functionalities


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0a1ad66..1448705 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -91,6 +91,8 @@ jobs:
 
       - name: Build Python package
         shell: bash
+        env:
+          RUSTFLAGS: "-C target-feature=+fxsr,+sse,+sse2,+sse3,+ssse3,+sse4.1,+sse4.2,+popcnt,+avx,+fma"
         run: |
           source $HOME/.cargo/env
           export PATH="/opt/python/cp38-cp38/bin/:$PATH"
@@ -128,6 +130,8 @@ jobs:
           python-version: "3.8"
           architecture: x64
       - name: Install Python dependencies and build
+        env:
+          RUSTFLAGS: "-C target-feature=+fxsr,+sse,+sse2,+sse3,+ssse3,+sse4.1,+sse4.2,+popcnt,+avx,+fma"
         run: |
           python -V
           pip install wheel
@@ -160,6 +164,8 @@ jobs:
       - name: Set LIBCLANG_PATH
         run: echo "LIBCLANG_PATH=$((gcm clang).source -replace "clang.exe")" >> $env:GITHUB_ENV
       - name: Install Python dependencies and build
+        env:
+          RUSTFLAGS: "-C target-feature=+fxsr,+sse,+sse2,+sse3,+ssse3,+sse4.1,+sse4.2,+popcnt,+avx,+fma"
         # binary dependencies without wheels on pypi are installed via conda
         run: |
           python -V
diff --git a/.gitignore b/.gitignore
index b6d1de6..a587eef 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,6 @@ flamegraph*.svg
 /venv
 __pycache__
 .pytest_cache
+.ruff_cache
+/docs/build/
+/docs/source/generated/
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 0000000..bad159b
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,16 @@
+version: 2
+
+sphinx:
+  builder: html
+
+build:
+  os: "ubuntu-22.04"
+  tools:
+    python: "3.9"
+    rust: "1.64"
+
+python:
+  install:
+    - requirements: docs/requirements.txt
+    - method: pip
+      path: .
diff --git a/CHANGES.md b/CHANGES.md
deleted file mode 100644
index 21fb73b..0000000
--- a/CHANGES.md
+++ /dev/null
@@ -1,98 +0,0 @@
-# Changelog
-
-All notable changes to this project will be documented in this file.
-
-The format is loosely based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres
-to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
-
-## h3ronpy [Unreleased]
-
-## h3ronpy [0.16.1] - 2023-02-16
-* Upgrade from pyo3 0.17 to 0.18.
-* Support minimum supported python version to 3.7 again
-
-## h3ronpy [0.16.0] - 2022-12-28
-
-* Directly support GeoSeries in vector to H3 conversion by automatically exchanging geometries using WKB. [#7](https://github.com/nmandery/h3ronpy/pull/7)
-* Raise minimum supported python version to 3.8.
-* `intersecting` argument for `geodataframe_to_h3` to also include cells which are only intersecting with the geometry, but whose centroid is not contained in the geometry.
-
-## h3ronpy [0.15.1] - 2022-10-28
-
-* Upgrade to h3ron v0.16.0.
-* Reduced the durations the GIL is held.
-
-## h3ronpy [0.15.0] - 2022-09-11
-
-- Dependency upgrades incl. upgrade to H3 v4.0.0 / h3ron v0.15.0.
-  - The migration to H3 v4.0.0 comes with renaming a few functions to follow H3 conventions:
-    - In `h3ronpy.op` module:
-      - `kring_distances` -> `grid_disk_distances`
-      - `kring_distances_agg` -> `grid_disk_distances_agg`
-      - `kring_distances_agg_np` -> `grid_disk_distances_agg_np`
-     
-## h3ronpy [0.14.0] - 2022-01-29
-
-### Added
-- Add `h3ronpy.op.change_resolution` and `h3ronpy.op.change_resolution_paired`
-
-### Changed
-- Replace usage of `wkb` crate with `geozero` because of licensing.
-- Stop supporting Python 3.6 (was EOL on 23 Dec 2021). Minimum supported python version now is Python 3.7.
-- Omit empty geometries when converting vector data.
-- Dependency upgrades.
-
-## h3ronpy [0.13.1] - 2021-12-13
-
-### Added 
-- Building wheels for Windows and Mac using github actions.
-
-## h3ronpy [0.13.0] - 2021-12-10
-### Added
-- `h3ronpy.op.kring_distances` and `h3ronpy.op.kring_distances_agg`.
-
-### Changed
-- The python extension has been removed from this repository and moved to its own repository at
-  [github.com/nmandery/h3ronpy](https://github.com/nmandery/h3ronpy).
-- Upgrade h3ron dependency to 0.13, h3ron-ndarray to 0.13.
-- Upgrade `pyo3` and `rust-numpy` dependencies to 0.15.
-- Raise `geopandas` version requirement from 0.8 to 0.10.
-
-## h3ronpy [0.12.0] - 2021-08-10
-### Changed
-- dependency updates
-
-## h3ronpy [0.11.0] - 2021-06-12
-### Added
-- Support for transforming `numpy.float32` and `numpy.float64` raster arrays to H3 dataframes by warping the array values in `OrderedFloat<T>`.
-
-### Changed
-- Fix `ValueError` when converting empty dataframes. [#17](https://github.com/nmandery/h3ron/issues/17)
-- Deprecate `h3ronpy.util.h3index_column_to_geodataframe` in favor of `h3ronpy.util.dataframe_to_geodataframe`.
-- Update dependencies: `geo-types` 0.6->0.7, `ndarray` 0.14->0.15
-
-## h3ronpy [0.10.0] - 2021-04-24
-### Added
-- Unittests for `raster_to_dataframe` and `geodataframe_to_h3` using `pytest`
-
-### Changed
-### Removed
-
-## h3ronpy [0.9.0] - 2021-04-11
-### Added
-- Integration with geopandas `GeoDataFrame` to convert the contained geometries to H3.
-- Update of `maturin` to 0.10.2
-
-### Changed
-- Simplified API of raster integration.
-
-## Earlier versions
-
-The changes done in earlier versions where not documented in this changelog and can only be reconstructed from the
-commits in git.
-
-[Unreleased]: https://github.com/nmandery/h3ronpy/compare/v0.12.0...HEAD
-[0.12.0]: https://github.com/nmandery/h3ronpy/compare/v0.11.0...v0.12.0
-[0.11.0]: https://github.com/nmandery/h3ronpy/compare/v0.10.0...v0.11.0
-[0.10.0]: https://github.com/nmandery/h3ronpy/compare/v0.9.0...v0.10.0
-[0.9.0]: https://github.com/nmandery/h3ronpy/compare/v0.8.1...v0.9.0
diff --git a/CHANGES.rst b/CHANGES.rst
new file mode 100644
index 0000000..94c6be9
--- /dev/null
+++ b/CHANGES.rst
@@ -0,0 +1,191 @@
+Changelog
+=========
+
+All notable changes to this project will be documented in this file.
+
+The format is loosely based on `Keep a
+Changelog <https://keepachangelog.com/en/1.0.0/>`__, and this project
+adheres to `Semantic
+Versioning <https://semver.org/spec/v2.0.0.html>`__.
+
+0.17.0 - Unreleased
+-------------------
+
+In this release the project migrated to the `arrow memory model <https://arrow.apache.org/>`_ and switched
+from the `h3ron library <https://github.com/nmandery/h3ron>`_ to the Rust implementation of H3 named `h3o <https://github.com/HydroniumLabs/h3o>`_.
+This comes along with safer code and `performance improvements <https://github.com/nmandery/rasterh3/issues/1>`_. Although `h3ronpy` is no longer
+build on the `h3ron` rust crate, the name remains.
+
+As a result of these migrations, the Python API has completely changed. Due to this it is easier to refer to new
+new documentation than listing the changes here. That is another aspect of these changes - there now exists a sphinx generated
+documentation.
+
+Other changes:
+
+-  Upgrade from pyo3 0.18 to 0.19.
+
+0.16.1 - 2023-02-16
+-------------------
+
+-  Upgrade from pyo3 0.17 to 0.18.
+-  Support minimum supported python version to 3.7 again
+
+0.16.0 - 2022-12-28
+--------------------
+
+-  Directly support GeoSeries in vector to H3 conversion by
+   automatically exchanging geometries using WKB.
+   `#7 <https://github.com/nmandery/h3ronpy/pull/7>`__
+-  Raise minimum supported python version to 3.8.
+-  ``intersecting`` argument for ``geodataframe_to_h3`` to also include
+   cells which are only intersecting with the geometry, but whose
+   centroid is not contained in the geometry.
+
+0.15.1 - 2022-10-28
+-------------------
+
+-  Upgrade to h3ron v0.16.0.
+-  Reduced the durations the GIL is held.
+
+0.15.0 - 2022-09-11
+-------------------
+
+-  Dependency upgrades incl. upgrade to H3 v4.0.0 / h3ron v0.15.0.
+
+   -  The migration to H3 v4.0.0 comes with renaming a few functions to
+      follow H3 conventions:
+
+      -  In ``h3ronpy.op`` module:
+
+         -  ``kring_distances`` -> ``grid_disk_distances``
+         -  ``kring_distances_agg`` -> ``grid_disk_distances_agg``
+         -  ``kring_distances_agg_np`` -> ``grid_disk_distances_agg_np``
+
+0.14.0 - 2022-01-29
+-------------------
+
+Added
+~~~~~
+
+-  Add ``h3ronpy.op.change_resolution`` and
+   ``h3ronpy.op.change_resolution_paired``
+
+Changed
+~~~~~~~
+
+-  Replace usage of ``wkb`` crate with ``geozero`` because of licensing.
+-  Stop supporting Python 3.6 (was EOL on 23 Dec 2021). Minimum
+   supported python version now is Python 3.7.
+-  Omit empty geometries when converting vector data.
+-  Dependency upgrades.
+
+0.13.1 - 2021-12-13
+-------------------
+
+.. _added-1:
+
+Added
+~~~~~
+
+-  Building wheels for Windows and Mac using github actions.
+
+0.13.0 - 2021-12-10
+-------------------
+
+.. _added-2:
+
+Added
+~~~~~
+
+-  ``h3ronpy.op.kring_distances`` and
+   ``h3ronpy.op.kring_distances_agg``.
+
+.. _changed-1:
+
+Changed
+~~~~~~~
+
+-  The python extension has been removed from this repository and moved
+   to its own repository at
+   `github.com/nmandery/h3ronpy <https://github.com/nmandery/h3ronpy>`__.
+-  Upgrade h3ron dependency to 0.13, h3ron-ndarray to 0.13.
+-  Upgrade ``pyo3`` and ``rust-numpy`` dependencies to 0.15.
+-  Raise ``geopandas`` version requirement from 0.8 to 0.10.
+
+0.12.0 - 2021-08-10
+-------------------
+
+.. _changed-2:
+
+Changed
+~~~~~~~
+
+-  dependency updates
+
+0.11.0 - 2021-06-12
+-------------------
+
+.. _added-3:
+
+Added
+~~~~~
+
+-  Support for transforming ``numpy.float32`` and ``numpy.float64``
+   raster arrays to H3 dataframes by warping the array values in
+   ``OrderedFloat<T>``.
+
+.. _changed-3:
+
+Changed
+~~~~~~~
+
+-  Fix ``ValueError`` when converting empty dataframes.
+   `#17 <https://github.com/nmandery/h3ron/issues/17>`__
+-  Deprecate ``h3ronpy.util.h3index_column_to_geodataframe`` in favor of
+   ``h3ronpy.util.dataframe_to_geodataframe``.
+-  Update dependencies: ``geo-types`` 0.6->0.7, ``ndarray`` 0.14->0.15
+
+
+0.10.0 - 2021-04-24
+-------------------
+
+.. _added-4:
+
+Added
+~~~~~
+
+-  Unittests for ``raster_to_dataframe`` and ``geodataframe_to_h3``
+   using ``pytest``
+
+.. _changed-4:
+
+Changed
+~~~~~~~
+
+Removed
+~~~~~~~
+
+0.9.0 - 2021-04-11
+------------------
+
+.. _added-5:
+
+Added
+~~~~~
+
+-  Integration with geopandas ``GeoDataFrame`` to convert the contained
+   geometries to H3.
+-  Update of ``maturin`` to 0.10.2
+
+.. _changed-5:
+
+Changed
+~~~~~~~
+
+-  Simplified API of raster integration.
+
+Earlier versions
+----------------
+
+The changes done in earlier versions where not documented in this
+changelog and can only be reconstructed from the commits in git.
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
deleted file mode 100644
index 01353fb..0000000
--- a/CONTRIBUTING.md
+++ /dev/null
@@ -1,36 +0,0 @@
-# Contributing
-
-Please take a moment to review this document to make the contribution process easy and effective
-for everyone involved.
-
-## License
-
-You agree to license your contribution under the
-[MIT License](LICENSE-MIT).
-
-## Submission guidelines
-
-1. Format your contribution with [rustfmt](https://github.com/rust-lang/rustfmt#quick-start).
-1. Resolve all warnings from [clippy](https://github.com/rust-lang/rust-clippy#usage).
-1. Document your contribution per the
-   [Rust documentation guidelines](https://doc.rust-lang.org/1.30.0/book/first-edition/documentation.html).
-
-Provide tests that prove the functionality you're contributing is correct.
-
-Ensure your contribution is free of all warnings and errors during compilation
-and testing.
-maintain quality.
-
-## Pull requests
-
-This project appreciates pull requests that fix bugs, make improvements or add new features.
-These are one of the main points of collaboration for the project.
-
-Pull requests should remain focused in scope and not contain unrelated commits.
-
-The project may use automated tools to verify submitted pull requests.
-Ensure you address any errors these tools may report.
-
-The project may choose to **not** accept your contribution.  When this happens we'll explain why
-and where possible make suggestions for how you might be able to achieve what you set out to
-with your pull request.
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
new file mode 100644
index 0000000..302436c
--- /dev/null
+++ b/CONTRIBUTING.rst
@@ -0,0 +1,44 @@
+Contributing
+============
+
+Please take a moment to review this document to make the contribution
+process easy and effective for everyone involved.
+
+License
+-------
+
+You agree to license your contribution under the `MIT
+License <LICENSE-MIT>`__.
+
+Submission guidelines
+---------------------
+
+1. Format your contribution with
+   `rustfmt <https://github.com/rust-lang/rustfmt#quick-start>`__.
+2. Resolve all warnings from
+   `clippy <https://github.com/rust-lang/rust-clippy#usage>`__.
+3. Document your contribution per the `Rust documentation
+   guidelines <https://doc.rust-lang.org/1.30.0/book/first-edition/documentation.html>`__.
+
+Provide tests that prove the functionality you’re contributing is
+correct.
+
+Ensure your contribution is free of all warnings and errors during
+compilation and testing. maintain quality.
+
+Pull requests
+-------------
+
+This project appreciates pull requests that fix bugs, make improvements
+or add new features. These are one of the main points of collaboration
+for the project.
+
+Pull requests should remain focused in scope and not contain unrelated
+commits.
+
+The project may use automated tools to verify submitted pull requests.
+Ensure you address any errors these tools may report.
+
+The project may choose to **not** accept your contribution. When this
+happens we’ll explain why and where possible make suggestions for how
+you might be able to achieve what you set out to with your pull request.
diff --git a/Cargo.toml b/Cargo.toml
index 8fe2318..66327ba 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,12 +1,11 @@
 [package]
 name = "h3ronpy"
-version = "0.16.1"
+version = "0.17.0"
 authors = ["Nico Mandery <nico@nmandery.net>"]
-description = "python bindings for h3ron"
+description = "Data science toolkit for the H3 geospatial grid"
 edition = "2021"
 license = "MIT"
-keywords = ["geo", "spatial", "h3", "ndarray", "numpy", "python"]
-readme = "README.md"
+keywords = ["geo", "spatial", "h3", "pandas", "geopandas", "arrow", "polars", "numpy", "python"]
 homepage = "https://github.com/nmandery/h3ronpy"
 repository = "https://github.com/nmandery/h3ronpy"
 
@@ -17,12 +16,21 @@ crate-type = ["cdylib"]
 [dependencies]
 env_logger = "^0.10"
 geo-types = "^0.7"
-geozero = { version = "0.9", default-features = false, features = ["with-geo", "with-wkb"] }
-h3ron = { git = "https://github.com/nmandery/h3ron.git", rev = "b72c4220f7fd0e787f2ea6e99b9e7fafdbfec6cd" }
-h3ron-ndarray = { git = "https://github.com/nmandery/h3ron.git", rev = "b72c4220f7fd0e787f2ea6e99b9e7fafdbfec6cd" }
+rasterh3 = { git = "https://github.com/nmandery/rasterh3.git", rev = "a8f46d56ae8197e1978385c13c30b3b8f5a932a0" }
+h3arrow = { git = "https://github.com/nmandery/h3arrow.git", rev = "2ed9cb76b73658b1209440edd810a1c198032299", features = ["geoarrow", "rayon"] }
+#h3arrow = { path = "/home/nicodev/CLionProjects/h3arrow", features = ["geoarrow", "rayon"] }
 ndarray = { version = "0.15", features = ["rayon"] }
-numpy = "0.18"
+numpy = "0.19"
 ordered-float = ">=2.0.1"
-py_geo_interface = { version = "0.5", features = ["f64", "wkb"] }
-pyo3 = { version = "0.18", features = ["extension-module", "abi3", "abi3-py37"] }
+py_geo_interface = { version = "0.6", features = ["f64", "wkb"] }
+pyo3 = { version = "0.19", features = ["extension-module", "abi3", "abi3-py38"] }
 rayon = "^1.5"
+hashbrown = "0.13"
+geo = "0.25"
+itertools = "0.10"
+
+
+[profile.release]
+lto = "thin"
+strip = true
+debug = false
\ No newline at end of file
diff --git a/MANIFEST.in b/MANIFEST.in
index 127ea3f..6190cd0 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,1 +1,1 @@
-exclude venv install-dev-dependencies.py documentation.ipynb requirements.documentation.txt data
+exclude venv install-dev-dependencies.py README.ipynb requirements.documentation.txt data
diff --git a/README.ipynb b/README.ipynb
deleted file mode 100644
index c10168a..0000000
--- a/README.ipynb
+++ /dev/null
@@ -1,421 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# `h3ronpy`\n",
-    "\n",
-    "[![PyPI version](https://img.shields.io/pypi/v/h3ronpy)](https://pypi.python.org/pypi/h3ronpy/) [![dependency status](https://deps.rs/repo/github/nmandery/h3ronpy/status.svg)](https://deps.rs/repo/github/nmandery/h3ronpy) [![DOI](https://zenodo.org/badge/402118389.svg)](https://zenodo.org/badge/latestdoi/402118389)\n",
-    "\n",
-    "Python extension for the [H3](https://h3geo.org) geospatial indexing system exposing some functionalities of the [h3ron](https://github.com/nmandery/h3ron) rust crates to the python language and integrating with the [numpy](https://numpy.org/), [pandas](https://pandas.pydata.org/),\n",
-    "[geopandas](https://geopandas.org/), [rasterio](https://rasterio.readthedocs.io/en/latest/) and [gdal](https://gdal.org/) libraries. \n",
-    "\n",
-    "One goal is to not duplicate any functions already implemented by the [official H3 python bindings](https://github.com/uber/h3-py).\n",
-    "\n",
-    "This library is in parts parallelized using [rayon](https://github.com/rayon-rs/rayon). The number of threads can be controlled as \n",
-    "described in [the rayon FAQ](https://github.com/rayon-rs/rayon/blob/master/FAQ.md#how-many-threads-will-rayon-spawn)\n",
-    "\n",
-    "\n",
-    "## Installation\n",
-    "\n",
-    "### via PyPI\n",
-    "\n",
-    "Regarding linux: to install the `manylinux_2_24` wheels, a fairly recent version of `pip` is required. So in case `pip` only finds the source distribution, upgrade `pip` using\n",
-    "\n",
-    "```\n",
-    "pip install --upgrade pip\n",
-    "```\n",
-    "\n",
-    "### Building manually\n",
-    "\n",
-    "To build this extension from source, you will need:\n",
-    "\n",
-    "* Rust. Install the latest version using [rustup](https://rustup.rs/)\n",
-    "* A C compiler for the libh3 sources, for example `clang`\n",
-    "* Python 3.x and the corresponding C headers\n",
-    "* The dependencies from the [requirements.dev.txt](./requirements.dev.txt) file.\n",
-    "\n",
-    "On Ubuntu most system-level dependencies should be available after running rustup and \n",
-    "\n",
-    "```shell\n",
-    "sudo apt-get install python3-dev build-essential clang\n",
-    "```\n",
-    "\n",
-    "The build is done using [maturin](https://github.com/PyO3/maturin)\n",
-    "\n",
-    "There are three main commands:\n",
-    "\n",
-    "* `maturin publish` builds the crate into python packages and publishes them to pypi.\n",
-    "* `maturin build` builds the wheels and stores them in a folder (`../target/wheels` by default), but doesn't upload them. It's possible to upload those with [twine](https://github.com/pypa/twine).\n",
-    "* `maturin develop` builds the crate and installs it as a python module directly in the current virtualenv.\n",
-    "\n",
-    "To build the extension just use the `maturin build --release` command for an optimized build.\n",
-    "\n",
-    "\n",
-    "## Logging\n",
-    "\n",
-    "This library uses rusts [`log` crate](https://docs.rs/log) together with the [`env_logger` crate](https://docs.rs/env_logger). \n",
-    "This means that logging to `stdout` can be  controlled via environment variables. Set `RUST_LOG` to `debug`, `error`, `info`, \n",
-    "`warn`, or `trace` for the corresponding log output.\n",
-    "\n",
-    "For more fine-grained logging settings, see the documentation of `env_logger`.\n",
-    "\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Usage\n",
-    "\n",
-    "The following sections depend on the libraries from the `requirements.documentation.txt` file being installed."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# imports\n",
-    "\n",
-    "from matplotlib import pyplot\n",
-    "import rasterio\n",
-    "from rasterio.plot import show\n",
-    "import numpy as np\n",
-    "import h3.api.numpy_int as h3\n",
-    "from scipy import ndimage\n",
-    "import geopandas as gpd\n",
-    "\n",
-    "# increase the plot size\n",
-    "pyplot.rcParams['figure.dpi'] = 100"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Converting raster data to H3\n",
-    "\n",
-    "#### Prepare a dataset using rasterio first"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "src = rasterio.open(\"data/europe-and-north-africa.tif\")\n",
-    "print(src.colorinterp)\n",
-    "\n",
-    "green = src.read(2)\n",
-    "blue = src.read(3)\n",
-    "print(green.shape)\n",
-    "\n",
-    "show(src)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Do image processing, like this messy extraction and smoothing of the greenish pixels done here:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "vegetation_mask = (green < 250) & (blue < 50)\n",
-    "ocean_mask = (green >= 6) & (green <= 14) & (blue >= 47) & (blue <= 54)\n",
-    "vegetation_nodata_value = 0\n",
-    "\n",
-    "vegetation = np.full(green.shape, 10, dtype=\"int8\")\n",
-    "vegetation[ocean_mask] = vegetation_nodata_value\n",
-    "vegetation[vegetation_mask] = 20\n",
-    "#pyplot.imshow(vegetation, cmap='Greens')\n",
-    "\n",
-    "# smooth a bit to remove single pixels\n",
-    "vegetation = ndimage.gaussian_filter(vegetation, sigma=.7)\n",
-    "vegetation[vegetation <= 5] = vegetation_nodata_value\n",
-    "vegetation[(vegetation > 0) & (vegetation < 15)] = 1\n",
-    "vegetation[vegetation >= 15] = 2\n",
-    "vegetation[ocean_mask] = vegetation_nodata_value\n",
-    "\n",
-    "vegetation_plot_args = dict(cmap='Greens', vmin=0, vmax=2)\n",
-    "\n",
-    "pyplot.imshow(vegetation, **vegetation_plot_args)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "vegetation"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Convert the raster numpy array to H3\n",
-    "\n",
-    "Find the h3 resolution to use. See also the docstrings of the used functions and of the `h3ronpy.raster` module."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "from h3ronpy.raster import nearest_h3_resolution\n",
-    "\n",
-    "h3_res = nearest_h3_resolution(vegetation.shape, src.transform, search_mode=\"smaller_than_pixel\")\n",
-    "print(f\"Using H3 resolution {h3_res}\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### ... directly into a pandas (geo-)dataframe:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "from h3ronpy.raster import raster_to_dataframe\n",
-    "\n",
-    "print(\"conversion start\")\n",
-    "vegetation_h3_df = raster_to_dataframe(vegetation, src.transform, h3_res, nodata_value=vegetation_nodata_value, compacted=True, geo=True)\n",
-    "print(\"conversion done\")\n",
-    "\n",
-    "print(\"plotting ... this may take a bit\")\n",
-    "vegetation_h3_df.plot(column=\"value\", linewidth=0.2, edgecolor=\"black\", **vegetation_plot_args)\n",
-    "pyplot.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Using GDAL instead of rasterio\n",
-    "\n",
-    "To not repeat any steps from above, we just save the processed data using `rasterio` and load it again using `gdal`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "# ... TBW ..."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Convert vectordata to H3 indexes\n",
-    "\n",
-    "The conversion of GeoDataFrames is parallelized using the available CPUs."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import geopandas as gpd\n",
-    "from h3ronpy import vector, util\n",
-    "\n",
-    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
-    "africa = world[world[\"continent\"] == \"Africa\"]\n",
-    "africa.plot(column=\"name\")\n",
-    "\n",
-    "# convert to h3 cells and build a new geodataframe\n",
-    "gdf = util.h3index_column_to_geodataframe(vector.geodataframe_to_h3(africa, 3))\n",
-    "gdf.plot(column=\"name\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Polygon to H3 conversion is based on the H3 `polyfill` function, which fills based on the containment of the H3 cells centroid in the polygon geometry. Depending on the shape of the geometry the resulting cells may look like below:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def fill_namibia(**kw):\n",
-    "    namibia = world[world[\"name\"] == \"Namibia\"]\n",
-    "\n",
-    "    cell_ax = util.h3index_column_to_geodataframe(vector.geodataframe_to_h3(namibia, 3, **kw)).plot()\n",
-    "\n",
-    "    return namibia.plot(ax=cell_ax, facecolor=(0,0,0,0), edgecolor='black')\n",
-    "\n",
-    "fill_namibia()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The `intersecting` argument extends the polyfill a bit to include all cells which intersect with the polygon geometry:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fill_namibia(intersecting=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Converting H3 indexes to vector geometries\n",
-    "\n",
-    "The main functionality already exists with the `LinkedPolygons` in the H3 library itself. This library here just integrates this and adds some polish - like the smoothing - around it.\n",
-    "\n",
-    "Additionally the polygon objects implement the [python geo interface](https://gist.github.com/sgillies/2217756), so they are compatible to most mapping libraries - see the geopandas plotting below."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from h3ronpy import Polygon\n",
-    "\n",
-    "squares_and_rings = np.concatenate([\n",
-    "    h3.hex_ring(h3.geo_to_h3(12.3, 15.4, 10), 6),\n",
-    "    h3.hex_ring(h3.geo_to_h3(12.31, 15.39, 10), 10),\n",
-    "    h3.hex_ring(h3.geo_to_h3(12.29, 15.395, 10), 8),\n",
-    "    h3.polyfill_polygon([\n",
-    "        (12.25, 15.28),\n",
-    "        (12.25, 15.37),\n",
-    "        (12.35, 15.37),\n",
-    "        (12.35, 15.28), \n",
-    "        (12.25, 15.28)\n",
-    "    ], 10, lnglat_order=False)\n",
-    "])\n",
-    "sqr_geoms = Polygon.from_h3indexes(squares_and_rings) # returns a list of polygons\n",
-    "gpd.GeoDataFrame({}, geometry=sqr_geoms).plot(edgecolor=\"black\", linewidth=0.5)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For a cleaner look, the generated vector geometries can also be smoothed using a variation of [Chaikins smoothing algorithm](http://www.idav.ucdavis.edu/education/CAGDNotes/Chaikins-Algorithm/Chaikins-Algorithm.html). This includes a non-destructive simplification where the number of vertices gets reduced."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "sqr_smoothed = Polygon.from_h3indexes(squares_and_rings, smoothen=True)\n",
-    "gpd.GeoDataFrame({}, geometry=sqr_smoothed).plot(edgecolor=\"black\", linewidth=0.5)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "There is also the `align_to_h3_resolution` variant of this which will generate polygons not exceeding the area of a given parent resolution.\n",
-    "\n",
-    "Corners will be aligned to the corners of the parent resolution when they are less than an\n",
-    "edge length away from them. This is to avoid gaps when `smoothen` is set to true."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "box_indexes = h3.polyfill_polygon([\n",
-    "        (12.25, 15.30),\n",
-    "        (12.25, 15.35),\n",
-    "        (12.30, 15.35),\n",
-    "        (12.30, 15.30), \n",
-    "        (12.25, 15.30)\n",
-    "], 10)\n",
-    "\n",
-    " # aling the geometries to h3 resolution 2 where complete edges touch. This reduced the number of vertices \n",
-    "# a fair amount\n",
-    "box_aligned_geoms = Polygon.from_h3indexes_aligned(box_indexes, 7, smoothen=True)\n",
-    "gpd.GeoDataFrame({}, geometry=box_aligned_geoms).plot(edgecolor=\"black\", linewidth=0.5)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Vectorized operations on h3indexes\n",
-    "\n",
-    "The `h3ronpy.op` module contains vectorized functions for working with h3indexes. Currently mostly for working with k-rings - more will propably be added in the future."
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.10.6"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/README.md b/README.md
deleted file mode 100644
index 986af18..0000000
--- a/README.md
+++ /dev/null
@@ -1,399 +0,0 @@
-# `h3ronpy`
-
-[![PyPI version](https://img.shields.io/pypi/v/h3ronpy)](https://pypi.python.org/pypi/h3ronpy/) [![dependency status](https://deps.rs/repo/github/nmandery/h3ronpy/status.svg)](https://deps.rs/repo/github/nmandery/h3ronpy) [![DOI](https://zenodo.org/badge/402118389.svg)](https://zenodo.org/badge/latestdoi/402118389)
-
-Python extension for the [H3](https://h3geo.org) geospatial indexing system exposing some functionalities of the [h3ron](https://github.com/nmandery/h3ron) rust crates to the python language and integrating with the [numpy](https://numpy.org/), [pandas](https://pandas.pydata.org/),
-[geopandas](https://geopandas.org/), [rasterio](https://rasterio.readthedocs.io/en/latest/) and [gdal](https://gdal.org/) libraries. 
-
-One goal is to not duplicate any functions already implemented by the [official H3 python bindings](https://github.com/uber/h3-py).
-
-This library is in parts parallelized using [rayon](https://github.com/rayon-rs/rayon). The number of threads can be controlled as 
-described in [the rayon FAQ](https://github.com/rayon-rs/rayon/blob/master/FAQ.md#how-many-threads-will-rayon-spawn)
-
-
-## Installation
-
-### via PyPI
-
-Regarding linux: to install the `manylinux_2_24` wheels, a fairly recent version of `pip` is required. So in case `pip` only finds the source distribution, upgrade `pip` using
-
-```
-pip install --upgrade pip
-```
-
-### Building manually
-
-To build this extension from source, you will need:
-
-* Rust. Install the latest version using [rustup](https://rustup.rs/)
-* A C compiler for the libh3 sources, for example `clang`
-* Python 3.x and the corresponding C headers
-* The dependencies from the [requirements.dev.txt](./requirements.dev.txt) file.
-
-On Ubuntu most system-level dependencies should be available after running rustup and 
-
-```shell
-sudo apt-get install python3-dev build-essential clang
-```
-
-The build is done using [maturin](https://github.com/PyO3/maturin)
-
-There are three main commands:
-
-* `maturin publish` builds the crate into python packages and publishes them to pypi.
-* `maturin build` builds the wheels and stores them in a folder (`../target/wheels` by default), but doesn't upload them. It's possible to upload those with [twine](https://github.com/pypa/twine).
-* `maturin develop` builds the crate and installs it as a python module directly in the current virtualenv.
-
-To build the extension just use the `maturin build --release` command for an optimized build.
-
-
-## Logging
-
-This library uses rusts [`log` crate](https://docs.rs/log) together with the [`env_logger` crate](https://docs.rs/env_logger). 
-This means that logging to `stdout` can be  controlled via environment variables. Set `RUST_LOG` to `debug`, `error`, `info`, 
-`warn`, or `trace` for the corresponding log output.
-
-For more fine-grained logging settings, see the documentation of `env_logger`.
-
-
-
-## Usage
-
-The following sections depend on the libraries from the `requirements.documentation.txt` file being installed.
-
-
-```python
-# imports
-
-from matplotlib import pyplot
-import rasterio
-from rasterio.plot import show
-import numpy as np
-import h3.api.numpy_int as h3
-from scipy import ndimage
-import geopandas as gpd
-
-# increase the plot size
-pyplot.rcParams['figure.dpi'] = 100
-```
-
-### Converting raster data to H3
-
-#### Prepare a dataset using rasterio first
-
-
-```python
-src = rasterio.open("data/europe-and-north-africa.tif")
-print(src.colorinterp)
-
-green = src.read(2)
-blue = src.read(3)
-print(green.shape)
-
-show(src)
-```
-
-    (<ColorInterp.red: 3>, <ColorInterp.green: 4>, <ColorInterp.blue: 5>)
-    (284, 327)
-
-
-
-    
-![png](README_files/README_4_1.png)
-    
-
-
-
-
-
-    <AxesSubplot: >
-
-
-
-Do image processing, like this messy extraction and smoothing of the greenish pixels done here:
-
-
-```python
-vegetation_mask = (green < 250) & (blue < 50)
-ocean_mask = (green >= 6) & (green <= 14) & (blue >= 47) & (blue <= 54)
-vegetation_nodata_value = 0
-
-vegetation = np.full(green.shape, 10, dtype="int8")
-vegetation[ocean_mask] = vegetation_nodata_value
-vegetation[vegetation_mask] = 20
-#pyplot.imshow(vegetation, cmap='Greens')
-
-# smooth a bit to remove single pixels
-vegetation = ndimage.gaussian_filter(vegetation, sigma=.7)
-vegetation[vegetation <= 5] = vegetation_nodata_value
-vegetation[(vegetation > 0) & (vegetation < 15)] = 1
-vegetation[vegetation >= 15] = 2
-vegetation[ocean_mask] = vegetation_nodata_value
-
-vegetation_plot_args = dict(cmap='Greens', vmin=0, vmax=2)
-
-pyplot.imshow(vegetation, **vegetation_plot_args)
-```
-
-
-
-
-    <matplotlib.image.AxesImage at 0x7f46942b8fd0>
-
-
-
-
-    
-![png](README_files/README_6_1.png)
-    
-
-
-
-```python
-vegetation
-```
-
-
-
-
-    array([[1, 1, 1, ..., 0, 0, 0],
-           [1, 1, 1, ..., 0, 0, 0],
-           [1, 1, 1, ..., 0, 0, 0],
-           ...,
-           [0, 0, 0, ..., 0, 0, 0],
-           [0, 0, 0, ..., 0, 0, 0],
-           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)
-
-
-
-#### Convert the raster numpy array to H3
-
-Find the h3 resolution to use. See also the docstrings of the used functions and of the `h3ronpy.raster` module.
-
-
-```python
-from h3ronpy.raster import nearest_h3_resolution
-
-h3_res = nearest_h3_resolution(vegetation.shape, src.transform, search_mode="smaller_than_pixel")
-print(f"Using H3 resolution {h3_res}")
-```
-
-    Using H3 resolution 5
-
-
-#### ... directly into a pandas (geo-)dataframe:
-
-
-```python
-from h3ronpy.raster import raster_to_dataframe
-
-print("conversion start")
-vegetation_h3_df = raster_to_dataframe(vegetation, src.transform, h3_res, nodata_value=vegetation_nodata_value, compacted=True, geo=True)
-print("conversion done")
-
-print("plotting ... this may take a bit")
-vegetation_h3_df.plot(column="value", linewidth=0.2, edgecolor="black", **vegetation_plot_args)
-pyplot.show()
-```
-
-    conversion start
-    conversion done
-    plotting ... this may take a bit
-
-
-
-    
-![png](README_files/README_11_1.png)
-    
-
-
-#### Using GDAL instead of rasterio
-
-To not repeat any steps from above, we just save the processed data using `rasterio` and load it again using `gdal`.
-
-
-```python
-# ... TBW ...
-```
-
-### Convert vectordata to H3 indexes
-
-The conversion of GeoDataFrames is parallelized using the available CPUs.
-
-
-```python
-import geopandas as gpd
-from h3ronpy import vector, util
-
-world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
-africa = world[world["continent"] == "Africa"]
-africa.plot(column="name")
-
-# convert to h3 cells and build a new geodataframe
-gdf = util.h3index_column_to_geodataframe(vector.geodataframe_to_h3(africa, 3))
-gdf.plot(column="name")
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_15_1.png)
-    
-
-
-
-    
-![png](README_files/README_15_2.png)
-    
-
-
-Polygon to H3 conversion is based on the H3 `polyfill` function, which fills based on the containment of the H3 cells centroid in the polygon geometry. Depending on the shape of the geometry the resulting cells may look like below:
-
-
-```python
-def fill_namibia(**kw):
-    namibia = world[world["name"] == "Namibia"]
-
-    cell_ax = util.h3index_column_to_geodataframe(vector.geodataframe_to_h3(namibia, 3, **kw)).plot()
-
-    return namibia.plot(ax=cell_ax, facecolor=(0,0,0,0), edgecolor='black')
-
-fill_namibia()
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_17_1.png)
-    
-
-
-The `intersecting` argument extends the polyfill a bit to include all cells which intersect with the polygon geometry:
-
-
-```python
-fill_namibia(intersecting=True)
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_19_1.png)
-    
-
-
-### Converting H3 indexes to vector geometries
-
-The main functionality already exists with the `LinkedPolygons` in the H3 library itself. This library here just integrates this and adds some polish - like the smoothing - around it.
-
-Additionally the polygon objects implement the [python geo interface](https://gist.github.com/sgillies/2217756), so they are compatible to most mapping libraries - see the geopandas plotting below.
-
-
-```python
-from h3ronpy import Polygon
-
-squares_and_rings = np.concatenate([
-    h3.hex_ring(h3.geo_to_h3(12.3, 15.4, 10), 6),
-    h3.hex_ring(h3.geo_to_h3(12.31, 15.39, 10), 10),
-    h3.hex_ring(h3.geo_to_h3(12.29, 15.395, 10), 8),
-    h3.polyfill_polygon([
-        (12.25, 15.28),
-        (12.25, 15.37),
-        (12.35, 15.37),
-        (12.35, 15.28), 
-        (12.25, 15.28)
-    ], 10, lnglat_order=False)
-])
-sqr_geoms = Polygon.from_h3indexes(squares_and_rings) # returns a list of polygons
-gpd.GeoDataFrame({}, geometry=sqr_geoms).plot(edgecolor="black", linewidth=0.5)
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_21_1.png)
-    
-
-
-For a cleaner look, the generated vector geometries can also be smoothed using a variation of [Chaikins smoothing algorithm](http://www.idav.ucdavis.edu/education/CAGDNotes/Chaikins-Algorithm/Chaikins-Algorithm.html). This includes a non-destructive simplification where the number of vertices gets reduced.
-
-
-```python
-sqr_smoothed = Polygon.from_h3indexes(squares_and_rings, smoothen=True)
-gpd.GeoDataFrame({}, geometry=sqr_smoothed).plot(edgecolor="black", linewidth=0.5)
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_23_1.png)
-    
-
-
-There is also the `align_to_h3_resolution` variant of this which will generate polygons not exceeding the area of a given parent resolution.
-
-Corners will be aligned to the corners of the parent resolution when they are less than an
-edge length away from them. This is to avoid gaps when `smoothen` is set to true.
-
-
-```python
-box_indexes = h3.polyfill_polygon([
-        (12.25, 15.30),
-        (12.25, 15.35),
-        (12.30, 15.35),
-        (12.30, 15.30), 
-        (12.25, 15.30)
-], 10)
-
- # aling the geometries to h3 resolution 2 where complete edges touch. This reduced the number of vertices 
-# a fair amount
-box_aligned_geoms = Polygon.from_h3indexes_aligned(box_indexes, 7, smoothen=True)
-gpd.GeoDataFrame({}, geometry=box_aligned_geoms).plot(edgecolor="black", linewidth=0.5)
-```
-
-
-
-
-    <AxesSubplot: >
-
-
-
-
-    
-![png](README_files/README_25_1.png)
-    
-
-
-### Vectorized operations on h3indexes
-
-The `h3ronpy.op` module contains vectorized functions for working with h3indexes. Currently mostly for working with k-rings - more will propably be added in the future.
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000..e311c10
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,37 @@
+h3ronpy
+=======
+
+A data science toolkit for the `H3 geospatial grid <https://h3geo.org/>`_.
+
+.. image:: https://img.shields.io/pypi/v/h3ronpy
+    :alt: PyPI
+    :target: https://pypi.python.org/pypi/h3ronpy/
+
+.. image:: https://readthedocs.org/projects/h3ronpy/badge/?version=latest
+    :alt: ReadTheDocs
+    :target: https://h3ronpy.readthedocs.io/
+
+.. image:: https://zenodo.org/badge/402118389.svg
+    :alt: DOI
+    :target: https://zenodo.org/badge/latestdoi/402118389
+
+
+This library is not a substitute for the official `python h3 library <https://github.com/uber/h3-py>`_ - instead it provides more
+high-level functions on top of H3 and integrations into common dataframe libraries.
+
+Features
+--------
+
+* Build on `Apache Arrow <https://arrow.apache.org>`_ and `pyarrow <https://arrow.apache.org/docs/python/index.html>`_ for efficient data handling.
+* Dedicated APIs for the the `pandas <https://pandas.pydata.org>`_ and `polars <https://www.pola.rs/>`_ dataframe libraries. The `pandas` support includes `geopandas <https://geopandas.org>`_.
+* Multi-threaded conversion of raster data to the H3 grid using `numpy arrays <https://numpy.org/>`_.
+* Multi-threaded conversion of vector data, including `geopandas` `GeoDataFrames` and any object which supports the python `__geo_interface__` protocol (`shapely`, `geojson`, ...).
+
+
+Documentation is available on `<https://h3ronpy.readthedocs.io/>`_.
+
+
+License
+-------
+
+MIT
\ No newline at end of file
diff --git a/README_files/README_11_1.png b/README_files/README_11_1.png
deleted file mode 100644
index 917a243..0000000
Binary files a/README_files/README_11_1.png and /dev/null differ
diff --git a/README_files/README_15_1.png b/README_files/README_15_1.png
deleted file mode 100644
index 6474e4f..0000000
Binary files a/README_files/README_15_1.png and /dev/null differ
diff --git a/README_files/README_15_2.png b/README_files/README_15_2.png
deleted file mode 100644
index 4e343e6..0000000
Binary files a/README_files/README_15_2.png and /dev/null differ
diff --git a/README_files/README_17_1.png b/README_files/README_17_1.png
deleted file mode 100644
index 5290059..0000000
Binary files a/README_files/README_17_1.png and /dev/null differ
diff --git a/README_files/README_19_1.png b/README_files/README_19_1.png
deleted file mode 100644
index 222a056..0000000
Binary files a/README_files/README_19_1.png and /dev/null differ
diff --git a/README_files/README_21_1.png b/README_files/README_21_1.png
deleted file mode 100644
index 8c0f133..0000000
Binary files a/README_files/README_21_1.png and /dev/null differ
diff --git a/README_files/README_23_1.png b/README_files/README_23_1.png
deleted file mode 100644
index 0c81a93..0000000
Binary files a/README_files/README_23_1.png and /dev/null differ
diff --git a/README_files/README_25_1.png b/README_files/README_25_1.png
deleted file mode 100644
index b575016..0000000
Binary files a/README_files/README_25_1.png and /dev/null differ
diff --git a/README_files/README_4_1.png b/README_files/README_4_1.png
deleted file mode 100644
index 0118d92..0000000
Binary files a/README_files/README_4_1.png and /dev/null differ
diff --git a/README_files/README_6_1.png b/README_files/README_6_1.png
deleted file mode 100644
index 1ac211c..0000000
Binary files a/README_files/README_6_1.png and /dev/null differ
diff --git a/data/naturalearth_110m_admin_0_countries.fgb b/data/naturalearth_110m_admin_0_countries.fgb
new file mode 100644
index 0000000..cc85a44
Binary files /dev/null and b/data/naturalearth_110m_admin_0_countries.fgb differ
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000..d0c3cbf
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = source
+BUILDDIR      = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/docs/make.bat b/docs/make.bat
new file mode 100644
index 0000000..9534b01
--- /dev/null
+++ b/docs/make.bat
@@ -0,0 +1,35 @@
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=source
+set BUILDDIR=build
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+	echo.
+	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+	echo.installed, then set the SPHINXBUILD environment variable to point
+	echo.to the full path of the 'sphinx-build' executable. Alternatively you
+	echo.may add the Sphinx directory to PATH.
+	echo.
+	echo.If you don't have Sphinx installed, grab it from
+	echo.http://sphinx-doc.org/
+	exit /b 1
+)
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+
+:end
+popd
diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 0000000..ed56745
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,8 @@
+sphinx
+sphinx_rtd_theme>=1.2
+sphinx-autodoc-typehints
+h3==3.7.6
+matplotlib
+rasterio==1.3.3
+jupyter-sphinx
+scipy
\ No newline at end of file
diff --git a/docs/source/api/arrow.rst b/docs/source/api/arrow.rst
new file mode 100644
index 0000000..9a97674
--- /dev/null
+++ b/docs/source/api/arrow.rst
@@ -0,0 +1,25 @@
+Arrow
+=====
+
+Functions
+---------
+
+.. automodule:: h3ronpy.arrow
+   :members:
+   :undoc-members:
+
+
+Raster module
+-------------
+
+.. automodule:: h3ronpy.arrow.raster
+   :members:
+   :undoc-members:
+
+
+Vector module
+-------------
+
+.. automodule:: h3ronpy.arrow.vector
+   :members:
+   :undoc-members:
diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst
new file mode 100644
index 0000000..f755736
--- /dev/null
+++ b/docs/source/api/index.rst
@@ -0,0 +1,8 @@
+API documentation
+=================
+
+.. toctree::
+
+   arrow
+   pandas
+   polars
diff --git a/docs/source/api/pandas.rst b/docs/source/api/pandas.rst
new file mode 100644
index 0000000..650cba3
--- /dev/null
+++ b/docs/source/api/pandas.rst
@@ -0,0 +1,25 @@
+Pandas/GeoPandas
+================
+
+Functions
+---------
+
+.. automodule:: h3ronpy.pandas
+   :members:
+   :undoc-members:
+
+
+Raster module
+-------------
+
+.. automodule:: h3ronpy.pandas.raster
+   :members:
+   :undoc-members:
+
+
+Vector module
+-------------
+
+.. automodule:: h3ronpy.pandas.vector
+   :members:
+   :undoc-members:
diff --git a/docs/source/api/polars.rst b/docs/source/api/polars.rst
new file mode 100644
index 0000000..fb5b371
--- /dev/null
+++ b/docs/source/api/polars.rst
@@ -0,0 +1,25 @@
+Polars
+======
+
+Functions
+---------
+
+.. automodule:: h3ronpy.polars
+   :members:
+   :undoc-members:
+
+
+Raster module
+-------------
+
+.. automodule:: h3ronpy.polars.raster
+   :members:
+   :undoc-members:
+
+
+Vector module
+-------------
+
+.. automodule:: h3ronpy.polars.vector
+   :members:
+   :undoc-members:
diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
new file mode 100644
index 0000000..d76c92b
--- /dev/null
+++ b/docs/source/changelog.rst
@@ -0,0 +1,1 @@
+.. include:: ../../CHANGES.rst
diff --git a/docs/source/conf.py b/docs/source/conf.py
new file mode 100644
index 0000000..ef83d6b
--- /dev/null
+++ b/docs/source/conf.py
@@ -0,0 +1,42 @@
+from h3ronpy import __version__ as hp_version
+
+# Configuration file for the Sphinx documentation builder.
+
+# -- Project information
+
+project = "h3ronpy"
+copyright = "2023, the h3ronpy authors"
+author = "Nico Mandery"
+
+
+release = hp_version
+version = hp_version
+
+# -- General configuration
+
+extensions = [
+    "sphinx.ext.duration",
+    "sphinx.ext.doctest",
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx.ext.intersphinx",
+    "sphinx_rtd_theme",
+    "jupyter_sphinx",
+]
+
+intersphinx_mapping = {
+    "python": ("https://docs.python.org/3/", None),
+    "sphinx": ("https://www.sphinx-doc.org/en/master/", None),
+}
+intersphinx_disabled_domains = ["std"]
+
+templates_path = ["_templates"]
+
+# -- Options for HTML output
+
+html_theme = "sphinx_rtd_theme"
+
+# -- Options for EPUB output
+epub_show_urls = "footnote"
+
+autodoc_typehints = "both"
diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst
new file mode 100644
index 0000000..ac7b6bc
--- /dev/null
+++ b/docs/source/contributing.rst
@@ -0,0 +1,1 @@
+.. include:: ../../CONTRIBUTING.rst
diff --git a/docs/source/index.rst b/docs/source/index.rst
new file mode 100644
index 0000000..266fcca
--- /dev/null
+++ b/docs/source/index.rst
@@ -0,0 +1,19 @@
+.. include:: ../../README.rst
+
+.. note::
+
+   This project is under active development.
+
+
+Contents
+--------
+
+.. toctree::
+    :maxdepth: 3
+
+    installation
+    usage/index
+    api/index
+    changelog
+    license
+    contributing
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
new file mode 100644
index 0000000..31d62b0
--- /dev/null
+++ b/docs/source/installation.rst
@@ -0,0 +1,42 @@
+Installation
+============
+
+.. note::
+
+    To avoid pulling in unused dependencies, `h3ronpy` does not declare a dependency to `pandas`,
+    `geopandas` and `polars`. These packages need to be installed separately.
+
+
+From PyPi
+---------
+
+.. code-block:: shell
+
+   pip install h3ronpy
+
+
+From source
+-----------
+
+To build from source a recent version of the `rust language <https://www.rust-lang.org/>`_ is required. The easiest
+way to install is by using `rustup <https://rustup.rs/>`_.
+
+An recent version of `pip` is required - version 23.1.2 works. `pip` can be upgraded by running
+
+.. code-block:: shell
+
+    pip install --upgrade pip
+
+
+After that, everything is set up to build and install `h3ronpy`:
+
+.. code-block:: shell
+
+    git clone https://github.com/nmandery/h3ronpy.git
+    cd h3ronpy
+    pip install .
+
+This will build the rust code using `maturin <https://www.maturin.rs/>`_. For more information on this see its website.
+
+When encountering a circular import error after this installation procedure, just change the directory out of the
+h3ronpy source directory.
diff --git a/docs/source/license.rst b/docs/source/license.rst
new file mode 100644
index 0000000..6f497ec
--- /dev/null
+++ b/docs/source/license.rst
@@ -0,0 +1,6 @@
+License
+=======
+
+MIT
+
+.. include:: ../../LICENSE.txt
\ No newline at end of file
diff --git a/docs/source/usage/grid.rst b/docs/source/usage/grid.rst
new file mode 100644
index 0000000..048d1aa
--- /dev/null
+++ b/docs/source/usage/grid.rst
@@ -0,0 +1,57 @@
+Grid traversal
+==============
+
+Create a few test cells to the examples on this page:
+
+.. jupyter-execute::
+
+    import numpy as np
+    import h3.api.numpy_int as h3
+    import pandas as pd
+
+    from h3ronpy.pandas.vector import cells_dataframe_to_geodataframe
+    from h3ronpy import DEFAULT_CELL_COLUMN_NAME
+
+    cells = np.array(
+        [
+            h3.geo_to_h3(5.2, -5.2, 7),
+            h3.geo_to_h3(5.3, -5.1, 7),
+        ],
+        dtype=np.uint64,
+    )
+
+
+Grid-disks with :py:func:`h3ronpy.pandas.grid_disk`
+---------------------------------------------------
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas import grid_disk
+
+    cells_dataframe_to_geodataframe(
+        pd.DataFrame({
+            DEFAULT_CELL_COLUMN_NAME: grid_disk(cells, 9, flatten=True).unique()}
+        )
+    ).plot()
+
+
+Grid-disk aggregates with :py:func:`h3ronpy.pandas.grid_disk_aggregate_k`
+-------------------------------------------------------------------------
+
+This builds ontop of :py:func:`h3ronpy.pandas.grid_disk_distances` while directly
+performing simple aggregations to avoid returning potentially very large dataframes.
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas import grid_disk_aggregate_k
+
+    cells_dataframe_to_geodataframe(
+        grid_disk_aggregate_k(cells, 9, "min")
+    ).plot(column="k", legend=True, legend_kwds={"label": "k", "orientation": "horizontal"},)
+
+
+.. jupyter-execute::
+
+    cells_dataframe_to_geodataframe(
+        grid_disk_aggregate_k(cells, 9, "max")
+    ).plot(column="k", legend=True, legend_kwds={"label": "k", "orientation": "horizontal"},)
diff --git a/docs/source/usage/index.rst b/docs/source/usage/index.rst
new file mode 100644
index 0000000..61e4440
--- /dev/null
+++ b/docs/source/usage/index.rst
@@ -0,0 +1,18 @@
+.. _usage:
+
+Usage examples
+==============
+
+We are mainly using `pandas` and `geopandas` here to be able to plot the geometries, but most non-geometry features are available
+in the other APIs as well.
+
+.. note::
+
+    These examples are supposed to show some of the provided functionalities. They certainly do not cover all aspects.
+
+.. toctree::
+    :maxdepth: 2
+
+    raster
+    vector
+    grid
diff --git a/docs/source/usage/raster.rst b/docs/source/usage/raster.rst
new file mode 100644
index 0000000..41696d9
--- /dev/null
+++ b/docs/source/usage/raster.rst
@@ -0,0 +1,90 @@
+Converting raster data
+======================
+
+.. jupyter-execute::
+
+    from matplotlib import pyplot
+    import rasterio
+    from rasterio.plot import show
+    import numpy as np
+    import h3.api.numpy_int as h3
+    from scipy import ndimage
+    import geopandas as gpd
+
+    # increase the plot size
+    pyplot.rcParams['figure.dpi'] = 120
+
+
+Prepare a dataset using rasterio first
+--------------------------------------
+
+.. jupyter-execute::
+
+    import rasterio
+    from rasterio.plot import show
+
+    src = rasterio.open("../data/europe-and-north-africa.tif")
+    print(src.colorinterp)
+
+    green = src.read(2)
+    blue = src.read(3)
+    print(green.shape)
+
+    show(src)
+
+Do some image processing - like this messy extraction of a vegetation mask here:
+
+.. jupyter-execute::
+
+    vegetation_mask = (green < 250) & (blue < 50)
+    ocean_mask = (green >= 6) & (green <= 14) & (blue >= 47) & (blue <= 54)
+    vegetation_nodata_value = 0
+
+    vegetation = np.full(green.shape, 10, dtype="int8")
+    vegetation[ocean_mask] = vegetation_nodata_value
+    vegetation[vegetation_mask] = 20
+
+    # smooth a bit to remove single pixels
+    vegetation = ndimage.gaussian_filter(vegetation, sigma=.7)
+    vegetation[vegetation <= 5] = vegetation_nodata_value
+    vegetation[(vegetation > 0) & (vegetation < 15)] = 1
+    vegetation[vegetation >= 15] = 2
+    vegetation[ocean_mask] = vegetation_nodata_value
+
+    vegetation_plot_args = dict(cmap='Greens', vmin=0, vmax=2)
+
+    pyplot.imshow(vegetation, **vegetation_plot_args)
+
+.. jupyter-execute::
+
+    vegetation
+
+Convert the raster numpy array to H3
+------------------------------------
+
+Find the closest H3 resolution to use. See also the docstrings of the used functions and of the `h3ronpy.pandas.raster` module.
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas.raster import nearest_h3_resolution
+
+    h3_res = nearest_h3_resolution(vegetation.shape, src.transform, search_mode="smaller_than_pixel")
+    print(f"Using H3 resolution {h3_res}")
+
+Now we convert the raster directly into a geopandas `GeoDataFrame`:
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas.raster import raster_to_dataframe
+
+    vegetation_h3_df = raster_to_dataframe(
+        vegetation,
+        src.transform,
+        h3_res,
+        nodata_value=vegetation_nodata_value,
+        compact=True,
+        geo=True
+    )
+
+    vegetation_h3_df.plot(column="value", linewidth=0.2, edgecolor="black", **vegetation_plot_args)
+    pyplot.show()
diff --git a/docs/source/usage/vector.rst b/docs/source/usage/vector.rst
new file mode 100644
index 0000000..7fb8910
--- /dev/null
+++ b/docs/source/usage/vector.rst
@@ -0,0 +1,84 @@
+Converting vector data
+======================
+
+.. jupyter-execute::
+
+    from matplotlib import pyplot
+
+    # increase the plot size
+    pyplot.rcParams['figure.dpi'] = 120
+
+
+.. jupyter-execute::
+
+    import geopandas as gpd
+
+    world = gpd.read_file("../data/naturalearth_110m_admin_0_countries.fgb")
+    africa = world[world["CONTINENT"] == "Africa"]
+    africa.plot(column="NAME_EN")
+
+
+Converting a complete `GeoDataFrame` to cells
+---------------------------------------------
+
+Includes building a new `GeoDataFrame` with the cell geometries using :py:func:`h3ronpy.pandas.geodataframe_to_cells`:
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas.vector import geodataframe_to_cells, cells_dataframe_to_geodataframe
+
+    df = geodataframe_to_cells(africa, 3)
+    gdf = cells_dataframe_to_geodataframe(df)
+    gdf.plot(column="NAME_EN")
+
+
+Polygon fill modes
+------------------
+
+Polygon to H3 conversion is based on centroid containment.
+Depending on the shape of the geometry the resulting cells may look like below:
+
+
+.. jupyter-execute::
+
+    namibia = africa[africa["NAME_EN"] == "Namibia"]
+
+    def fill_namibia(**kw):
+        cell_ax = cells_dataframe_to_geodataframe(geodataframe_to_cells(namibia, 3, **kw)).plot()
+        return namibia.plot(ax=cell_ax, facecolor=(0,0,0,0), edgecolor='black')
+
+    fill_namibia()
+
+The `all_intersecting` argument extends this to include all cells which intersect with the polygon geometry:
+
+.. jupyter-execute::
+
+    fill_namibia(all_intersecting=True)
+
+
+Merging cells into larger polygons
+----------------------------------
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas.vector import cells_to_polygons, geoseries_to_cells
+
+    gpd.GeoDataFrame({
+        "geometry": cells_to_polygons(geoseries_to_cells(namibia.geometry, 3, flatten=True), link_cells=True)
+    }).plot()
+
+
+Single geometries
+-----------------
+
+It is also possible to convert single `shapely` geometries or any other type providing the python `__geo_interface__`:
+
+.. jupyter-execute::
+
+    from h3ronpy.pandas.vector import geometry_to_cells
+
+    namibia_geom = namibia["geometry"].iloc[0]
+    print(namibia_geom)
+    geometry_to_cells(namibia_geom, 3)
+
+
diff --git a/h3ronpy/__init__.py b/h3ronpy/__init__.py
index f66ed15..d18f690 100644
--- a/h3ronpy/__init__.py
+++ b/h3ronpy/__init__.py
@@ -1,10 +1,11 @@
-"""
-h3ron python bindings
-"""
-
-# import from rust library
-from .h3ronpy import version, Polygon
+from . import h3ronpy as _native
+from .h3ronpy import version, DEFAULT_CELL_COLUMN_NAME  # noqa: F401
 
 __version__ = version()
 
 H3_CRS = "EPSG:4326"
+
+if not _native.is_release_build():
+    import warnings
+
+    warnings.warn("h3ronpy has not been compiled in release mode. Performance will be degraded.", RuntimeWarning)
diff --git a/h3ronpy/arrow/__init__.py b/h3ronpy/arrow/__init__.py
new file mode 100644
index 0000000..445c879
--- /dev/null
+++ b/h3ronpy/arrow/__init__.py
@@ -0,0 +1,157 @@
+from typing import Union
+
+import pyarrow as pa
+from ..h3ronpy import op
+
+try:
+    import polars as pl
+
+    _HAS_POLARS = True
+except ImportError:
+    _HAS_POLARS = False
+
+
+def _to_arrow_array(arr, dtype) -> pa.Array:
+    converted = None
+    if _HAS_POLARS:
+        if isinstance(arr, pl.Series):
+            converted = arr.to_arrow()
+
+    if converted is None:
+        converted = pa.array(arr, type=dtype)
+
+    if isinstance(arr, pa.ChunkedArray):
+        converted = converted.combine_chunks()
+    return converted
+
+
+def _to_uint64_array(arr) -> pa.Array:
+    return _to_arrow_array(arr, pa.uint64())
+
+
+def change_resolution(arr, resolution: int) -> pa.Array:
+    return op.change_resolution(_to_uint64_array(arr), resolution)
+
+
+def change_resolution_paired(arr, resolution: int) -> pa.Table:
+    """
+    Returns a table/dataframe with two columns: `cell_before` and `cell_after`
+    with the cells h3index before and after the resolution change.
+
+    This can be helpful when joining data in different resolutions via
+    dataframe libraries
+    """
+    return op.change_resolution_paired(_to_uint64_array(arr), resolution)
+
+
+def cells_resolution(arr) -> pa.Array:
+    """
+    Generates a new array containing the resolution of each cell of the
+    input array.
+
+    :param arr:
+    :return:
+    """
+    return op.cells_resolution(_to_uint64_array(arr))
+
+
+def cells_parse(arr, set_failing_to_invalid: bool = False) -> pa.Array:
+    """
+    Parse H3 cells from string arrays
+
+    Setting `set_failing_to_invalid` to true will trigger setting the validity bitmap according
+    the successful parsing of an individual element. Having this set to false will cause the
+    method to fail upon encountering the first unparsable value.
+    """
+    return op.cells_parse(_to_arrow_array(arr, pa.utf8()), set_failing_to_invalid=set_failing_to_invalid)
+
+
+def compact(arr, mixed_resolutions: bool = False) -> pa.Array:
+    """
+    Compact the given cells
+
+    The cells are expected to be of the same resolution, otherwise this operation will fail unless
+    `mixed_resolutions` is set to True. Setting this may lead to slight slow-downs.
+    """
+    return op.compact(_to_uint64_array(arr), mixed_resolutions=mixed_resolutions)
+
+
+def uncompact(arr, target_resolution: int) -> pa.Array:
+    """
+    Uncompact the given cells to the resolution `target_resolution`.
+
+    All higher resolution cells contained in the input array than the given `target_resolution` will
+    be omitted from the output.
+    """
+    return op.uncompact(_to_uint64_array(arr), target_resolution)
+
+
+def _make_h3index_valid_wrapper(fn, h3index_name, wrapper_name):
+    def valid_wrapper(arr, booleanarray: bool = False) -> pa.Array:
+        return fn(_to_uint64_array(arr), booleanarray=booleanarray)
+
+    valid_wrapper.__doc__ = f"""
+    Validate an array of potentially invalid {h3index_name} values by returning a new
+    UInt64 array with the validity mask set accordingly.
+    
+    If `booleanarray` is set to True, a boolean array describing the validity will be
+    returned instead.
+    """
+    valid_wrapper.__name__ = wrapper_name
+    return valid_wrapper
+
+
+cells_valid = _make_h3index_valid_wrapper(op.cells_valid, "cell", "cells_valid")
+vertexes_valid = _make_h3index_valid_wrapper(op.cells_valid, "vertex", "vertexes_valid")
+directededges_valid = _make_h3index_valid_wrapper(op.cells_valid, "directed edge", "directededges_valid")
+
+
+def grid_disk(cellarray, k: int, flatten: bool = False) -> Union[pa.ListArray, pa.Array]:
+    return op.grid_disk(_to_uint64_array(cellarray), k, flatten=flatten)
+
+
+def grid_disk_distances(cellarray, k: int, flatten: bool = False) -> pa.Table:
+    return op.grid_disk_distances(_to_uint64_array(cellarray), k, flatten=flatten)
+
+
+def grid_disk_aggregate_k(cellarray, k: int, aggregation_method: str) -> pa.Table:
+    """
+    Valid values for `aggregation_method` are `"min"` and `"max"`.
+    """
+    return op.grid_disk_aggregate_k(_to_uint64_array(cellarray), k, aggregation_method)
+
+
+def grid_ring_distances(cellarray, k_min: int, k_max: int, flatten: bool = False) -> pa.Table:
+    return op.grid_ring_distances(_to_uint64_array(cellarray), k_min, k_max, flatten=flatten)
+
+
+def cells_area_m2(cellarray) -> pa.Array:
+    return op.cells_area_m2(_to_uint64_array(cellarray))
+
+
+def cells_area_km2(cellarray) -> pa.Array:
+    return op.cells_area_km2(_to_uint64_array(cellarray))
+
+
+def cells_area_rads2(cellarray) -> pa.Array:
+    return op.cells_area_rads2(_to_uint64_array(cellarray))
+
+
+__all__ = [
+    change_resolution.__name__,
+    change_resolution_paired.__name__,
+    cells_resolution.__name__,
+    cells_parse.__name__,
+    compact.__name__,
+    uncompact.__name__,
+    cells_valid.__name__,
+    vertexes_valid.__name__,
+    directededges_valid.__name__,
+    grid_disk.__name__,
+    grid_disk_distances.__name__,
+    grid_ring_distances.__name__,
+    grid_disk_aggregate_k.__name__,
+    cells_area_m2.__name__,
+    cells_area_km2.__name__,
+    cells_area_rads2.__name__,
+]
diff --git a/h3ronpy/raster.py b/h3ronpy/arrow/raster.py
similarity index 55%
rename from h3ronpy/raster.py
rename to h3ronpy/arrow/raster.py
index d93dc4f..a0c5635 100644
--- a/h3ronpy/raster.py
+++ b/h3ronpy/arrow/raster.py
@@ -1,22 +1,18 @@
 """
-Conversion of raster numpy arrays to h3.
+Conversion of raster `numpy` arrays to H3 cells
 
 Resolution search modes
 -----------------------
 
-* "min_diff": chose the h3 resolution where the difference in the area of a pixel and the h3index is as small as possible.
-* "smaller_than_pixel":  chose the h3 resolution where the area of the h3index is smaller than the area of a pixel.
+* "min_diff": chose the H3 resolution where the difference in the area of a pixel and the h3index is as small as possible.
+* "smaller_than_pixel":  chose the H3 resolution where the area of the h3index is smaller than the area of a pixel.
+
 """
 
-import geopandas as gp
-import h3.api.numpy_int as h3
+from ..h3ronpy import raster
+from .. import DEFAULT_CELL_COLUMN_NAME
 import numpy as np
-import pandas as pd
-from shapely.geometry import Polygon
-
-from . import H3_CRS
-from .h3ronpy import raster
-from .util import _array_nditer
+import pyarrow as pa
 
 try:
     # affine library is used by rasterio
@@ -41,9 +37,9 @@ def _get_transform(t):
     raise ValueError("unsupported object for transform")
 
 
-def nearest_h3_resolution(shape, transform, axis_order="yx", search_mode="min_diff"):
+def nearest_h3_resolution(shape, transform, axis_order="yx", search_mode="min_diff") -> int:
     """
-    find the h3 resolution closed to the size of a pixel in an array
+    Find the H3 resolution closest to the size of a pixel in an array
     of the given shape with the given transform
 
     :param shape: dimensions of the 2d array
@@ -55,23 +51,29 @@ def nearest_h3_resolution(shape, transform, axis_order="yx", search_mode="min_di
     return raster.nearest_h3_resolution(shape, _get_transform(transform), axis_order, search_mode)
 
 
-def raster_to_dataframe(in_raster: np.array, transform, h3_resolution: int, nodata_value=None, axis_order: str = "yx",
-                        compacted: bool = True, geo: bool = False):
+def raster_to_dataframe(
+    in_raster: np.array,
+    transform,
+    h3_resolution: int,
+    nodata_value=None,
+    axis_order: str = "yx",
+    compact: bool = True,
+) -> pa.Table:
     """
-    convert a raster/array to a pandas dataframe containing H3 indexes
+    Convert a raster/array to a pandas `DataFrame` containing H3 cell indexes
 
     This function is parallelized and uses the available CPUs by distributing tiles to a thread pool.
 
     The input geometry must be in WGS84.
 
-    :param in_raster: input 2-d array
-    :param transform:  the affine transformation
-    :param nodata_value: the nodata value. For these cells of the array there will be no h3 indexes generated
-    :param axis_order: axis order of the 2d array. Either "xy" or "yx"
-    :param h3_resolution: target h3 resolution
-    :param compacted: return compacted h3 indexes (see H3 docs)
-    :param geo: return a geopandas geodataframe with geometries. increases the memory usage.
-    :return: pandas dataframe or geodataframe
+    :param in_raster: Input 2D array
+    :param transform:  The affine transformation
+    :param nodata_value: The nodata value. For these cells of the array there will be no h3 indexes generated
+    :param axis_order: Axis order of the 2d array. Either "xy" or "yx"
+    :param h3_resolution: Target h3 resolution
+    :param compact: Return compacted h3 indexes (see H3 docs). This results in mixed H3 resolutions, but also can
+            reduce the amount of required memory.
+    :return: Tuple of arrow arrays
     """
 
     dtype = in_raster.dtype
@@ -99,25 +101,7 @@ def raster_to_dataframe(in_raster: np.array, transform, h3_resolution: int, noda
     else:
         raise NotImplementedError(f"no raster_to_h3 implementation for dtype {dtype.name}")
 
-    values, indexes = func(in_raster, _get_transform(transform), h3_resolution, axis_order, compacted, nodata_value)
-    if geo:
-        return gp.GeoDataFrame({
-            "h3index": indexes,
-            "value": values,
-            "geometry": [Polygon(h3.h3_to_geo_boundary(h, geo_json=True)) for h in _array_nditer(indexes)],
-        }, crs=H3_CRS)
-    else:
-        return pd.DataFrame({
-            "h3index": indexes,
-            "value": values
-        })
-
-
-def raster_to_geodataframe(*a, **kw):
-    """
-    convert to a geodataframe
-
-    Uses the same parameters as array_to_dataframe
-    """
-    kw["geo"] = True
-    return raster_to_dataframe(*a, **kw)
+    return pa.Table.from_arrays(
+        arrays=func(in_raster, _get_transform(transform), h3_resolution, axis_order, compact, nodata_value),
+        names=["value", DEFAULT_CELL_COLUMN_NAME],
+    )
diff --git a/h3ronpy/arrow/util.py b/h3ronpy/arrow/util.py
new file mode 100644
index 0000000..a004369
--- /dev/null
+++ b/h3ronpy/arrow/util.py
@@ -0,0 +1,15 @@
+import pyarrow.compute as pc
+import pyarrow as pa
+
+
+# from https://issues.apache.org/jira/browse/ARROW-12099
+def explode_table_include_null(table: pa.Table, column: str) -> pa.Table:
+    other_columns = list(table.schema.names)
+    other_columns.remove(column)
+    indices = pc.list_parent_indices(pc.fill_null(table[column], [None]))
+    result = table.select(other_columns).take(indices)
+    result = result.append_column(
+        pa.field(column, table.schema.field(column).type.value_type),
+        pc.list_flatten(pc.fill_null(table[column], [None])),
+    )
+    return result
diff --git a/h3ronpy/arrow/vector.py b/h3ronpy/arrow/vector.py
new file mode 100644
index 0000000..5101cea
--- /dev/null
+++ b/h3ronpy/arrow/vector.py
@@ -0,0 +1,85 @@
+from ..h3ronpy import vector
+from . import _to_uint64_array, _HAS_POLARS
+from typing import Optional, Tuple, Union
+import pyarrow as pa
+
+
+def cells_to_coordinates(arr, radians: bool = False) -> pa.Table:
+    """
+    convert to point coordinates in degrees
+    """
+    return vector.cells_to_coordinates(_to_uint64_array(arr), radians=radians)
+
+
+def cells_bounds(arr) -> Optional[Tuple]:
+    """
+    Bounds of the complete array as a tuple `(minx, miny, maxx, maxy)`.
+    """
+    return vector.cells_bounds(_to_uint64_array(arr))
+
+
+def cells_bounds_arrays(arr) -> pa.Table:
+    """
+    Build a table/dataframe with the columns `minx`, `miny`, `maxx` and `maxy` containing the bounds of the individual
+    cells from the input array.
+    """
+    return vector.cells_bounds_arrays(_to_uint64_array(arr))
+
+
+def cells_to_wkb_polygons(arr, radians: bool = False, link_cells: bool = False) -> pa.Array:
+    return vector.cells_to_wkb_polygons(_to_uint64_array(arr), radians=radians, link_cells=link_cells)
+
+
+def cells_to_wkb_points(arr, radians: bool = False) -> pa.Array:
+    return vector.cells_to_wkb_points(_to_uint64_array(arr), radians=radians)
+
+
+def vertexes_to_wkb_points(arr, radians: bool = False) -> pa.Array:
+    return vector.vertexes_to_wkb_points(_to_uint64_array(arr), radians=radians)
+
+
+def directededges_to_wkb_lines(arr, radians: bool = False) -> pa.Array:
+    return vector.directededges_to_wkb_lines(_to_uint64_array(arr), radians=radians)
+
+
+def directededges_to_wkb_linestrings(arr, radians: bool = False) -> pa.Array:
+    return vector.directededges_to_wkb_linestrings(_to_uint64_array(arr), radians=radians)
+
+
+def wkb_to_cells(
+    arr, resolution: int, compact: bool = False, all_intersecting: bool = True, flatten: bool = False
+) -> Union[pa.Array, pa.ListArray]:
+    """Convert a Series/Array/List of WKB values to H3 cells"""
+    if _HAS_POLARS:
+        import polars as pl
+
+        if isinstance(arr, pl.Series):
+            arr = arr.to_arrow()
+
+    if not isinstance(arr, pa.LargeBinaryArray):
+        arr = pa.array(arr, type=pa.large_binary())
+
+    if isinstance(arr, pa.ChunkedArray):
+        arr = arr.combine_chunks()
+    return vector.wkb_to_cells(arr, resolution, compact=compact, all_intersecting=all_intersecting, flatten=flatten)
+
+
+def geometry_to_cells(obj, resolution: int, compact: bool = False, all_intersecting: bool = True) -> pa.Array:
+    """
+    Convert a single object which supports the python `__geo_interface__` protocol to H3 cells
+    """
+    return vector.geometry_to_cells(obj, resolution, compact=compact, all_intersecting=all_intersecting)
+
+
+__all__ = [
+    cells_to_coordinates.__name__,
+    cells_bounds.__name__,
+    cells_bounds_arrays.__name__,
+    cells_to_wkb_polygons.__name__,
+    cells_to_wkb_points.__name__,
+    vertexes_to_wkb_points.__name__,
+    directededges_to_wkb_linestrings.__name__,
+    directededges_to_wkb_lines.__name__,
+    wkb_to_cells.__name__,
+    geometry_to_cells.__name__,
+]
diff --git a/h3ronpy/op.py b/h3ronpy/op.py
deleted file mode 100644
index ffa34bb..0000000
--- a/h3ronpy/op.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import typing
-
-import numpy as np
-import pandas as pd
-
-from .h3ronpy import op as native_op
-
-
-def grid_disk_distances(h3indexes: np.ndarray, k_max: int, k_min: int = 0) -> pd.DataFrame:
-    """
-    Vectorized grid-disk building
-
-    Returns a dataframe with the columns `h3index` (the ring center), `ring_h3index` and `ring_k`.
-
-    :param h3indexes:
-    :param k_max:
-    :param k_min:
-    :return:
-    """
-    h3index, ring_h3index, ring_k = native_op.grid_disk_distances(h3indexes, k_min, k_max)
-    return pd.DataFrame({
-        "h3index": h3index,
-        "ring_h3index": ring_h3index,
-        "ring_k": ring_k
-    })
-
-
-def grid_disk_distances_agg_np(h3indexes: np.ndarray, k_max: int, k_min: int = 0, aggregation_method: str = 'min') -> \
-        typing.Tuple[np.ndarray, np.ndarray]:
-    """
-    Vectorized grid-disk building, with the k-values of the rings being aggregated to their `min` or
-    `max` value for each cell.
-
-    :param h3indexes:
-    :param k_max:
-    :param k_min:
-    :param aggregation_method:
-    :return:
-    """
-    return native_op.grid_disk_distances_agg(h3indexes, k_min, k_max, aggregation_method)
-
-
-def grid_disk_distances_agg(h3indexes: np.ndarray, k_max: int, k_min: int = 0,
-                            aggregation_method: str = 'min') -> pd.DataFrame:
-    h3indexes_out, k_out = grid_disk_distances_agg_np(h3indexes, k_max, k_min=k_min, aggregation_method=aggregation_method)
-    return pd.DataFrame({"h3index": h3indexes_out, "k": k_out})
-
-
-def change_resolution(h3indexes: np.ndarray, h3_resolution: int) -> np.ndarray:
-    """Change the resolution of the given `h3indexes` to `h3_resolution`"""
-    return native_op.change_resolution(h3indexes, h3_resolution)
-
-
-def change_resolution_paired(h3indexes: np.ndarray, h3_resolution: int) -> pd.DataFrame:
-    """Change the resolution of the given `h3indexes` to `h3_resolution`
-
-    Returns a dataframe with a column for the `h3index_before` and a `h3index_after` the
-    resolution change.
-    """
-    h3index_before, h3index_after = native_op.change_resolution_paired(h3indexes, h3_resolution)
-    return pd.DataFrame({
-        "h3index_before": h3index_before,
-        "h3index_after": h3index_after,
-    })
diff --git a/h3ronpy/pandas/__init__.py b/h3ronpy/pandas/__init__.py
new file mode 100644
index 0000000..19ccbc8
--- /dev/null
+++ b/h3ronpy/pandas/__init__.py
@@ -0,0 +1,69 @@
+"""
+API to use `h3ronpy` with the `pandas dataframe library <https://pandas.pydata.org>`_ including `geopandas <https://geopandas.org>`_.
+
+.. warning::
+
+    To avoid pulling in unused dependencies, `h3ronpy` does not declare a dependency to `pandas` and `geopandas`. These
+    packages need to be installed separately.
+
+"""
+from .. import arrow as _arrow
+import pyarrow as pa
+from functools import wraps
+import pandas as pd
+
+
+def _wrap(func, ret_type=None):
+    @wraps(func)
+    def wrapper(*args, **kw):
+        result = func(*args, **kw)
+        if isinstance(result, pa.Table):
+            return result.to_pandas(split_blocks=True, self_destruct=True)
+        elif isinstance(result, pa.Array):
+            return result.to_pandas()
+        return result
+
+    if ret_type:
+        # create a copy to avoid modifying the dict of the wrapped function
+        wrapper.__annotations__ = dict(**wrapper.__annotations__)
+        wrapper.__annotations__["return"] = ret_type
+    return wrapper
+
+
+change_resolution = _wrap(_arrow.change_resolution, ret_type=pd.Series)
+change_resolution.__annotations__["return"] = pd.Series
+change_resolution_paired = _wrap(_arrow.change_resolution_paired, ret_type=pd.DataFrame)
+cells_resolution = _wrap(_arrow.cells_resolution, ret_type=pd.Series)
+cells_parse = _wrap(_arrow.cells_parse, ret_type=pd.Series)
+compact = _wrap(_arrow.compact, ret_type=pd.Series)
+uncompact = _wrap(_arrow.uncompact, ret_type=pd.Series)
+cells_valid = _wrap(_arrow.cells_valid, ret_type=pd.Series)
+vertexes_valid = _wrap(_arrow.vertexes_valid, ret_type=pd.Series)
+directededges_valid = _wrap(_arrow.directededges_valid, ret_type=pd.Series)
+grid_disk = _wrap(_arrow.grid_disk, ret_type=pd.Series)
+grid_disk_distances = _wrap(_arrow.grid_disk_distances, ret_type=pd.DataFrame)
+grid_ring_distances = _wrap(_arrow.grid_ring_distances, ret_type=pd.DataFrame)
+grid_disk_aggregate_k = _wrap(_arrow.grid_disk_aggregate_k, ret_type=pd.DataFrame)
+cells_area_m2 = _wrap(_arrow.cells_area_m2, ret_type=pd.Series)
+cells_area_km2 = _wrap(_arrow.cells_area_km2, ret_type=pd.Series)
+cells_area_rads2 = _wrap(_arrow.cells_area_rads2, ret_type=pd.Series)
+
+
+__all__ = [
+    change_resolution.__name__,
+    change_resolution_paired.__name__,
+    cells_resolution.__name__,
+    cells_parse.__name__,
+    compact.__name__,
+    uncompact.__name__,
+    cells_valid.__name__,
+    vertexes_valid.__name__,
+    directededges_valid.__name__,
+    grid_disk.__name__,
+    grid_disk_distances.__name__,
+    grid_ring_distances.__name__,
+    grid_disk_aggregate_k.__name__,
+    cells_area_m2.__name__,
+    cells_area_km2.__name__,
+    cells_area_rads2.__name__,
+]
diff --git a/h3ronpy/pandas/raster.py b/h3ronpy/pandas/raster.py
new file mode 100644
index 0000000..158946a
--- /dev/null
+++ b/h3ronpy/pandas/raster.py
@@ -0,0 +1,58 @@
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import typing
+
+from ..arrow import raster as arrow_raster
+from .vector import cells_dataframe_to_geodataframe
+
+__doc__ = arrow_raster.__doc__
+
+nearest_h3_resolution = arrow_raster.nearest_h3_resolution
+
+
+def raster_to_dataframe(
+    in_raster: np.array,
+    transform,
+    h3_resolution: int,
+    nodata_value=None,
+    axis_order: str = "yx",
+    compact: bool = True,
+    geo: bool = False,
+) -> typing.Union[gpd.GeoDataFrame, pd.DataFrame]:
+    """
+    Convert a raster/array to a pandas `DataFrame` containing H3 indexes
+
+    This function is parallelized and uses the available CPUs by distributing tiles to a thread pool.
+
+    The input geometry must be in WGS84.
+
+    :param in_raster: input 2-d array
+    :param transform:  the affine transformation
+    :param nodata_value: the nodata value. For these cells of the array there will be no h3 indexes generated
+    :param axis_order: axis order of the 2d array. Either "xy" or "yx"
+    :param h3_resolution: target h3 resolution
+    :param compact: Return compacted h3 indexes (see H3 docs). This results in mixed H3 resolutions, but also can
+            reduce the amount of required memory.
+    :param geo: Return a geopandas `GeoDataFrame` with geometries. increases the memory usage.
+    :return: pandas `DataFrame` or `GeoDataFrame`
+    """
+
+    df = arrow_raster.raster_to_dataframe(
+        in_raster, transform, h3_resolution, nodata_value=nodata_value, axis_order=axis_order, compact=compact
+    ).to_pandas()
+
+    if geo:
+        return cells_dataframe_to_geodataframe(df)
+    else:
+        return df
+
+
+def raster_to_geodataframe(*a, **kw) -> gpd.GeoDataFrame:
+    """
+    convert to a geodataframe
+
+    Uses the same parameters as array_to_dataframe
+    """
+    kw["geo"] = True
+    return raster_to_dataframe(*a, **kw)
diff --git a/h3ronpy/pandas/vector.py b/h3ronpy/pandas/vector.py
new file mode 100644
index 0000000..60ffedd
--- /dev/null
+++ b/h3ronpy/pandas/vector.py
@@ -0,0 +1,108 @@
+from . import _wrap
+from ..arrow import vector as _av
+from ..arrow import util as _arrow_util
+import pyarrow as pa
+import pandas as pd
+import geopandas as gpd
+from functools import wraps
+from .. import H3_CRS, DEFAULT_CELL_COLUMN_NAME
+
+
+def _geoseries_from_wkb(func):
+    @wraps(func)
+    def wrapper(*args, **kw):
+        return gpd.GeoSeries.from_wkb(func(*args, **kw), crs=H3_CRS)
+
+    # create a copy to avoid modifying the dict of the wrapped function
+    wrapper.__annotations__ = dict(**wrapper.__annotations__)
+    wrapper.__annotations__["return"] = gpd.GeoSeries
+
+    return wrapper
+
+
+cells_to_coordinates = _wrap(_av.cells_to_coordinates, ret_type=pd.DataFrame)
+cells_bounds = _av.cells_bounds
+cells_bounds_arrays = _wrap(_av.cells_bounds_arrays, ret_type=pd.DataFrame)
+cells_to_wkb_polygons = _wrap(_av.cells_to_wkb_polygons, ret_type=pd.Series)
+cells_to_polygons = _geoseries_from_wkb(cells_to_wkb_polygons)
+cells_to_wkb_points = _wrap(_av.cells_to_wkb_points, ret_type=pd.Series)
+cells_to_points = _geoseries_from_wkb(cells_to_wkb_points)
+vertexes_to_wkb_points = _wrap(_av.vertexes_to_wkb_points, ret_type=pd.Series)
+vertexes_to_points = _geoseries_from_wkb(vertexes_to_wkb_points)
+directededges_to_wkb_lines = _wrap(_av.directededges_to_wkb_lines, ret_type=pd.Series)
+directededges_to_lines = _geoseries_from_wkb(directededges_to_wkb_lines)
+directededges_to_wkb_linestrings = _wrap(_av.directededges_to_wkb_linestrings, ret_type=pd.Series)
+directededges_to_linestrings = _geoseries_from_wkb(directededges_to_wkb_linestrings)
+wkb_to_cells = _wrap(_av.wkb_to_cells, ret_type=pd.Series)
+geometry_to_cells = _wrap(_av.geometry_to_cells, ret_type=pd.Series)
+
+
+@wraps(wkb_to_cells)
+def geoseries_to_cells(geoseries: gpd.GeoSeries, *args, **kw):
+    return _av.wkb_to_cells(geoseries.to_wkb(), *args, **kw).to_pandas()
+
+
+geoseries_to_cells.__name__ = "geoseries_to_cells"
+
+
+def cells_dataframe_to_geodataframe(
+    df: pd.DataFrame, cell_column_name: str = DEFAULT_CELL_COLUMN_NAME
+) -> gpd.GeoDataFrame:
+    """
+    Convert a dataframe with a column containing cells to a geodataframe
+
+    :param df: input dataframe
+    :param cell_column_name: name of the column containing the h3 indexes
+    :return: GeoDataFrame
+    """
+    return gpd.GeoDataFrame(df, geometry=cells_to_polygons(df[cell_column_name]), crs=H3_CRS)
+
+
+def geodataframe_to_cells(
+    gdf: gpd.GeoDataFrame,
+    resolution: int,
+    compact: bool = False,
+    cell_column_name: str = DEFAULT_CELL_COLUMN_NAME,
+    all_intersecting: bool = False,
+) -> pd.DataFrame:
+    """
+    Convert a `GeoDataFrame` to H3 cells while exploding all other columns according to the number of cells derived
+    from the rows geometry.
+
+    The conversion of GeoDataFrames is parallelized using the available CPUs.
+
+    The duplication of all non-cell columns leads to increased memory requirements. Depending on the use-case
+    some of the more low-level conversion functions should be preferred.
+
+    :param gdf:
+    :param resolution:
+    :param compact:
+    :param cell_column_name:
+    :param all_intersecting:
+    :return:
+    """
+    cells = _av.wkb_to_cells(gdf.geometry.to_wkb(), resolution, compact=compact, all_intersecting=all_intersecting)
+    table = pa.Table.from_pandas(pd.DataFrame(gdf.drop(columns="geometry"))).append_column(cell_column_name, cells)
+    return _arrow_util.explode_table_include_null(table, cell_column_name).to_pandas().reset_index(drop=True)
+
+
+__all__ = [
+    cells_to_coordinates.__name__,
+    cells_bounds.__name__,
+    cells_bounds_arrays.__name__,
+    cells_to_wkb_polygons.__name__,
+    cells_to_polygons.__name__,
+    cells_to_wkb_points.__name__,
+    cells_to_points.__name__,
+    vertexes_to_wkb_points.__name__,
+    vertexes_to_points.__name__,
+    directededges_to_wkb_linestrings.__name__,
+    directededges_to_linestrings.__name__,
+    directededges_to_wkb_lines.__name__,
+    directededges_to_lines.__name__,
+    cells_dataframe_to_geodataframe.__name__,
+    wkb_to_cells.__name__,
+    geometry_to_cells.__name__,
+    geoseries_to_cells.__name__,
+    geodataframe_to_cells.__name__,
+]
diff --git a/h3ronpy/polars/__init__.py b/h3ronpy/polars/__init__.py
new file mode 100644
index 0000000..3a50337
--- /dev/null
+++ b/h3ronpy/polars/__init__.py
@@ -0,0 +1,66 @@
+"""
+API to use `h3ronpy` with the `polars dataframe library <https://www.pola.rs/>`_.
+
+.. warning::
+
+    To avoid pulling in unused dependencies, `h3ronpy` does not declare a dependency to `polars`. This
+    package needs to be installed separately.
+
+"""
+
+from functools import wraps
+import polars as pl
+import pyarrow as pa
+from .. import arrow as _arrow
+
+
+def _wrap(func, ret_type=None):
+    @wraps(func, updated=())
+    def wrapper(*args, **kw):
+        result = func(*args, **kw)
+        if isinstance(result, pa.Table) or isinstance(result, pa.Array):
+            return pl.from_arrow(result)
+        return result
+
+    if ret_type:
+        # create a copy to avoid modifying the dict of the wrapped function
+        wrapper.__annotations__ = dict(**wrapper.__annotations__)
+        wrapper.__annotations__["return"] = ret_type
+    return wrapper
+
+
+change_resolution = _wrap(_arrow.change_resolution, ret_type=pl.Series)
+change_resolution_paired = _wrap(_arrow.change_resolution_paired, ret_type=pl.DataFrame)
+cells_resolution = _wrap(_arrow.cells_resolution, ret_type=pl.Series)
+cells_parse = _wrap(_arrow.cells_parse, ret_type=pl.Series)
+compact = _wrap(_arrow.compact, ret_type=pl.Series)
+uncompact = _wrap(_arrow.uncompact, ret_type=pl.Series)
+cells_valid = _wrap(_arrow.cells_valid, ret_type=pl.Series)
+vertexes_valid = _wrap(_arrow.vertexes_valid, ret_type=pl.Series)
+directededges_valid = _wrap(_arrow.directededges_valid, ret_type=pl.Series)
+grid_disk = _wrap(_arrow.grid_disk, ret_type=pl.Series)
+grid_disk_distances = _wrap(_arrow.grid_disk_distances, ret_type=pl.DataFrame)
+grid_ring_distances = _wrap(_arrow.grid_ring_distances, ret_type=pl.DataFrame)
+grid_disk_aggregate_k = _wrap(_arrow.grid_disk_aggregate_k, ret_type=pl.DataFrame)
+cells_area_m2 = _wrap(_arrow.cells_area_m2, ret_type=pl.Series)
+cells_area_km2 = _wrap(_arrow.cells_area_km2, ret_type=pl.Series)
+cells_area_rads2 = _wrap(_arrow.cells_area_rads2, ret_type=pl.Series)
+
+__all__ = [
+    change_resolution.__name__,
+    change_resolution_paired.__name__,
+    cells_resolution.__name__,
+    cells_parse.__name__,
+    compact.__name__,
+    uncompact.__name__,
+    cells_valid.__name__,
+    vertexes_valid.__name__,
+    directededges_valid.__name__,
+    grid_disk.__name__,
+    grid_disk_distances.__name__,
+    grid_ring_distances.__name__,
+    grid_disk_aggregate_k.__name__,
+    cells_area_m2.__name__,
+    cells_area_km2.__name__,
+    cells_area_rads2.__name__,
+]
diff --git a/h3ronpy/polars/raster.py b/h3ronpy/polars/raster.py
new file mode 100644
index 0000000..b362c63
--- /dev/null
+++ b/h3ronpy/polars/raster.py
@@ -0,0 +1,10 @@
+from ..arrow import raster as arrow_raster
+from . import _wrap
+import polars as pl
+
+nearest_h3_resolution = arrow_raster.nearest_h3_resolution
+raster_to_dataframe = _wrap(arrow_raster.raster_to_dataframe, ret_type=pl.DataFrame)
+
+__doc__ = arrow_raster.__doc__
+
+__all__ = [nearest_h3_resolution.__name__, raster_to_dataframe.__name__]
diff --git a/h3ronpy/polars/vector.py b/h3ronpy/polars/vector.py
new file mode 100644
index 0000000..6e3c809
--- /dev/null
+++ b/h3ronpy/polars/vector.py
@@ -0,0 +1,27 @@
+from . import _wrap
+from ..arrow import vector as _av
+import polars as pl
+
+cells_to_coordinates = _wrap(_av.cells_to_coordinates, ret_type=pl.DataFrame)
+cells_bounds = _av.cells_bounds
+cells_bounds_arrays = _wrap(_av.cells_bounds_arrays, ret_type=pl.DataFrame)
+cells_to_wkb_polygons = _wrap(_av.cells_to_wkb_polygons, ret_type=pl.Series)
+cells_to_wkb_points = _wrap(_av.cells_to_wkb_points, ret_type=pl.Series)
+vertexes_to_wkb_points = _wrap(_av.vertexes_to_wkb_points, ret_type=pl.Series)
+directededges_to_wkb_lines = _wrap(_av.directededges_to_wkb_lines, ret_type=pl.Series)
+directededges_to_wkb_linestrings = _wrap(_av.directededges_to_wkb_linestrings, ret_type=pl.Series)
+wkb_to_cells = _wrap(_av.wkb_to_cells, ret_type=pl.Series)
+geometry_to_cells = _wrap(_av.geometry_to_cells, ret_type=pl.Series)
+
+__all__ = [
+    cells_to_coordinates.__name__,
+    cells_bounds.__name__,
+    cells_bounds_arrays.__name__,
+    cells_to_wkb_polygons.__name__,
+    cells_to_wkb_points.__name__,
+    vertexes_to_wkb_points.__name__,
+    directededges_to_wkb_lines.__name__,
+    directededges_to_wkb_linestrings.__name__,
+    wkb_to_cells.__name__,
+    geometry_to_cells.__name__,
+]
diff --git a/h3ronpy/util.py b/h3ronpy/util.py
deleted file mode 100644
index 16ff303..0000000
--- a/h3ronpy/util.py
+++ /dev/null
@@ -1,54 +0,0 @@
-import warnings
-from typing import Union, List, Iterable
-
-import geopandas as gpd
-import h3.api.numpy_int as h3
-import numpy as np
-import pandas as pd
-from shapely.geometry import Polygon
-
-from . import H3_CRS
-
-
-def h3indexes_to_geodataframe(
-        h3indexes: Union[np.array, List[int], Iterable[int]]) -> gpd.GeoDataFrame:
-    """
-    convert a numpy-array with uint64 h3 indexes to a geodataframe
-    """
-    i_arr = np.asarray(h3indexes)
-    return gpd.GeoDataFrame({
-        "geometry": [Polygon(h3.h3_to_geo_boundary(h, geo_json=True)) for h in _array_nditer(i_arr)],
-        "h3resolution": [h3.h3_get_resolution(h) for h in _array_nditer(i_arr)],
-        "h3index": [h3.h3_to_string(h) for h in _array_nditer(i_arr)]
-    }, crs=H3_CRS)
-
-
-def h3index_column_to_geodataframe(df: pd.DataFrame, column_name: str = "h3index") -> gpd.GeoDataFrame:
-    """
-    convert a dataframe with a column containing h3indexes to a geodataframe
-
-    :param df: input dataframe
-    :param column_name: name of the column containing the h3 indexes
-    :return: GeoDataFrame
-    """
-    warnings.warn("h3index_column_to_geodataframe has been deprecated in favor of dataframe_to_geodataframe",
-                  DeprecationWarning)
-    return dataframe_to_geodataframe(df, column_name=column_name)
-
-
-def dataframe_to_geodataframe(df: pd.DataFrame, column_name: str = "h3index") -> gpd.GeoDataFrame:
-    """
-    convert a dataframe with a column containing h3indexes to a geodataframe
-
-    :param df: input dataframe
-    :param column_name: name of the column containing the h3 indexes
-    :return: GeoDataFrame
-    """
-    return gpd.GeoDataFrame(df,
-                            geometry=[Polygon(h3.h3_to_geo_boundary(h, geo_json=True)) for h in
-                                      _array_nditer(df[column_name].to_numpy())],
-                            crs=H3_CRS)
-
-
-def _array_nditer(a: np.array):
-    return np.nditer(a, flags=["zerosize_ok", ])
diff --git a/h3ronpy/vector.py b/h3ronpy/vector.py
deleted file mode 100644
index 16e44b2..0000000
--- a/h3ronpy/vector.py
+++ /dev/null
@@ -1,74 +0,0 @@
-from typing import Tuple, Generator
-
-import geopandas as gpd
-import math
-import numpy as np
-import pandas as pd
-
-from .h3ronpy import vector
-
-
-def geometries_to_h3_generator(geometries: np.array, ids: np.array, h3_resolution: int, do_compact: bool = False,
-                               chunk_size: int = 1000, intersecting: bool = False) -> Generator[
-    Tuple[np.array, np.array], None, None]:
-    """
-    Generator to convert shapely geometries and ids to two numpy arrays with h3indexes and the correlated Ids.
-    Yields (ids, h3indexes)-tuples.
-
-    This function is parallelized and uses the available CPUs by distributing the geometries to a thread pool.
-
-    :param geometries: numpy-array with shapely geometries (in WGS84).
-    :param ids: numpy array with ids. must be unit64
-    :param h3_resolution: H3 resolution to use.
-    :param do_compact: Compacts the h3index when this is set to `True`. Default is `False`
-    :param chunk_size: Number of geometries to include in the yielded chunks
-    :param intersecting: Also include cells which are only intersecting with the geometry, but whose centroid is not contained in the geometry.
-    :return: None
-    """
-    num_chunks = math.ceil(len(geometries) / chunk_size)
-    for (chunk_geometries, chunk_ids) in zip(np.array_split(geometries, num_chunks),
-                                             np.array_split(ids, num_chunks)):
-        (ids, h3indexes) = vector.wkbbytes_with_ids_to_h3(chunk_ids, chunk_geometries, h3_resolution, do_compact,
-                                                          intersecting)
-        yield ids, h3indexes
-
-
-def geodataframe_to_h3(df: gpd.GeoDataFrame, h3_resolution: int, do_compact: bool = False, geometry_column="geometry",
-                       index_column_name="to_h3_idx", chunk_size=1000, intersecting: bool = False):
-    """
-    convert the geometries of a geodataframe to h3 indexes.
-
-    This function is parallelized and uses the available CPUs by distributing the geometries to a thread pool.
-
-    The geometry column will be removed from the resulting dataframe. The input geometry must be in WGS84.
-
-    :param df: The input geodataframe
-    :param h3_resolution: H3 resolution to use.
-    :param do_compact: Compacts the h3index when this is set to `True`. Default is `False`
-    :param geometry_column: The name of the column containing the geometry. Defaults to `geometry`.
-    :param index_column_name: The name for a temporary column used to join the H3 data to the input dataframe
-    :param chunk_size:
-    :param intersecting: Also include cells which are only intersecting with the geometry, but whose centroid is not contained in the geometry.
-    :return:
-    """
-    # add a column with a sequence to merge later with
-    df.insert(loc=0, column=index_column_name, value=np.arange(len(df), dtype=np.uint64))
-    dataframes = []
-    for (ids, h3indexes) in geometries_to_h3_generator(df[geometry_column].to_numpy(), df[index_column_name].to_numpy(),
-                                                       h3_resolution,
-                                                       do_compact=do_compact, chunk_size=chunk_size,
-                                                       intersecting=intersecting):
-        dataframes.append(pd.DataFrame({
-            index_column_name: ids,
-            "h3index": h3indexes
-        }))
-
-    if not dataframes:
-        return pd.DataFrame({})
-    output_df = pd.DataFrame(df.drop(columns=geometry_column)) \
-        .merge(pd.concat(dataframes), on=index_column_name)
-
-    # remove the column used for the merge again
-    del output_df[index_column_name]
-
-    return output_df
diff --git a/install-dev-dependencies.py b/install-dev-dependencies.py
index 967a482..b12e27c 100644
--- a/install-dev-dependencies.py
+++ b/install-dev-dependencies.py
@@ -18,7 +18,7 @@ def install(packages, upgrade=False):
 
 
 if __name__ == "__main__":
-    install(["pip", "toml", "black", "patchelf"], upgrade=True)  # always upgrade pip
+    install(["pip", "toml", "black", "ruff"], upgrade=True)  # always upgrade pip
 
     packages = []
 
diff --git a/justfile b/justfile
index dbc5509..aeb7ed3 100644
--- a/justfile
+++ b/justfile
@@ -1,6 +1,7 @@
 # to be executed via https://github.com/casey/just
 
-readme:
-    jupyter nbconvert --to Markdown --execute README.ipynb
-    git add README_files/*
-    git add README.md
+black:
+    black -l 120 h3ronpy tests *.py docs/source/*.py
+
+ruff:
+    ruff check h3ronpy tests *.py docs/source/*.py
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 7c53363..9e08357 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,17 +1,7 @@
 [build-system]
 requires = [
-    "maturin==0.14.7",
-]
-requires-dist = [
-    "numpy",
-    "geopandas>=0.10",
-    "h3>=3.7",
-    "Shapely>=1.7",
-]
-classifiers = [
-    "Programming Language :: Python :: 3",
-    "Topic :: Scientific/Engineering :: GIS",
-    "License :: OSI Approved :: MIT License",
+    "maturin==1.1",
+    "patchelf",
 ]
 build-backend = "maturin"
 
@@ -22,3 +12,30 @@ testpaths = [
     "h3ronpy",
     "tests"
 ]
+
+[tool.ruff]
+# Never enforce `E501` (line length violations).
+ignore = ["E501"]
+
+[project]
+name = "h3ronpy"
+
+dependencies = [
+    "numpy",
+    "pyarrow>=10.0",
+    "Shapely>=1.7",
+]
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "Topic :: Scientific/Engineering :: GIS",
+    "License :: OSI Approved :: MIT License",
+]
+
+[project.optional-dependencies]
+test = [
+    "rasterio",
+    "polars",
+    "geopandas>=0.10",
+    "pytest>=6",
+    "h3>=3.7",
+]
\ No newline at end of file
diff --git a/requirements.documentation.txt b/requirements.documentation.txt
deleted file mode 100644
index deea74f..0000000
--- a/requirements.documentation.txt
+++ /dev/null
@@ -1,9 +0,0 @@
-notebook 
-rasterio 
-ipympl
-numpy
-geopandas
-h3
-descartes
-scipy
-nbconvert
diff --git a/src/arrow_interop.rs b/src/arrow_interop.rs
new file mode 100644
index 0000000..bd3e0b8
--- /dev/null
+++ b/src/arrow_interop.rs
@@ -0,0 +1,105 @@
+use std::any::{type_name, Any};
+
+use h3arrow::array::{CellIndexArray, DirectedEdgeIndexArray, VertexIndexArray};
+use h3arrow::export::arrow2::array::{Array, UInt64Array};
+use h3arrow::export::arrow2::datatypes::Field;
+use h3arrow::export::arrow2::ffi;
+use pyo3::exceptions::PyValueError;
+use pyo3::ffi::Py_uintptr_t;
+use pyo3::prelude::*;
+use pyo3::Python;
+
+use crate::error::{IntoPyErr, IntoPyResult};
+
+pub(crate) fn with_pyarrow<F, O>(f: F) -> PyResult<O>
+where
+    F: FnOnce(Python, &PyModule) -> PyResult<O>,
+{
+    Python::with_gil(|py| {
+        let pyarrow = py.import("pyarrow")?;
+        f(py, pyarrow)
+    })
+}
+
+/// Arrow array to Python.
+/// from https://github.com/pola-rs/polars/blob/d1e5b1062c6872cd030b04b96505d2fac36b5376/py-polars/src/arrow_interop/to_py.rs
+pub(crate) fn native_to_pyarray(
+    array: Box<dyn Array>,
+    py: Python,
+    pyarrow: &PyModule,
+) -> PyResult<PyObject> {
+    let schema = Box::new(ffi::export_field_to_c(&Field::new(
+        "",
+        array.data_type().clone(),
+        true,
+    )));
+    let array = Box::new(ffi::export_array_to_c(array));
+
+    let schema_ptr: *const ffi::ArrowSchema = &*schema;
+    let array_ptr: *const ffi::ArrowArray = &*array;
+
+    let array = pyarrow.getattr("Array")?.call_method1(
+        "_import_from_c",
+        (array_ptr as Py_uintptr_t, schema_ptr as Py_uintptr_t),
+    )?;
+
+    Ok(array.to_object(py))
+}
+
+pub(crate) fn pyarray_to_native<T: Any + Array + Clone>(obj: &PyAny) -> PyResult<T> {
+    // prepare a pointer to receive the Array struct
+    let array = Box::new(ffi::ArrowArray::empty());
+    let schema = Box::new(ffi::ArrowSchema::empty());
+
+    let array_ptr = &*array as *const ffi::ArrowArray;
+    let schema_ptr = &*schema as *const ffi::ArrowSchema;
+
+    // make the conversion through PyArrow's private API
+    // this changes the pointer's memory and is thus unsafe. In particular, `_export_to_c` can go out of bounds
+    obj.call_method1(
+        "_export_to_c",
+        (array_ptr as Py_uintptr_t, schema_ptr as Py_uintptr_t),
+    )?;
+
+    let array = unsafe {
+        let field = ffi::import_field_from_c(schema.as_ref()).into_pyresult()?;
+        ffi::import_array_from_c(*array, field.data_type).into_pyresult()?
+    };
+
+    // downcast to the concrete type
+    Ok(array
+        .as_any()
+        .downcast_ref::<T>()
+        .ok_or_else(|| {
+            PyValueError::new_err(format!(
+                "Expected {}, found arrow array of type {:?}",
+                type_name::<T>(),
+                array.data_type(),
+            ))
+        })?
+        .clone())
+}
+
+pub(crate) fn pyarray_to_cellindexarray(obj: &PyAny) -> PyResult<CellIndexArray> {
+    pyarray_to_h3array::<CellIndexArray>(obj)
+}
+
+pub(crate) fn pyarray_to_vertexindexarray(obj: &PyAny) -> PyResult<VertexIndexArray> {
+    pyarray_to_h3array::<VertexIndexArray>(obj)
+}
+
+pub(crate) fn pyarray_to_directededgeindexarray(obj: &PyAny) -> PyResult<DirectedEdgeIndexArray> {
+    pyarray_to_h3array::<DirectedEdgeIndexArray>(obj)
+}
+
+pub(crate) fn pyarray_to_uint64array(obj: &PyAny) -> PyResult<UInt64Array> {
+    pyarray_to_native::<UInt64Array>(obj)
+}
+
+fn pyarray_to_h3array<T>(obj: &PyAny) -> PyResult<T>
+where
+    T: TryFrom<UInt64Array>,
+    <T as TryFrom<UInt64Array>>::Error: IntoPyErr,
+{
+    T::try_from(pyarray_to_uint64array(obj)?).into_pyresult()
+}
diff --git a/src/collections.rs b/src/collections.rs
deleted file mode 100644
index 677247d..0000000
--- a/src/collections.rs
+++ /dev/null
@@ -1,64 +0,0 @@
-use pyo3::prelude::*;
-
-use crate::cells_to_h3indexes;
-use crate::error::IntoPyResult;
-use h3ron::error::check_valid_h3_resolution;
-use h3ron::{collections as h3c, H3Cell, Index};
-use numpy::{IntoPyArray, PyArray1};
-
-#[pyclass]
-pub struct H3CompactedVec {
-    pub(crate) inner: h3c::CompactedCellVec,
-}
-
-#[pymethods]
-impl H3CompactedVec {
-    fn len(&self) -> usize {
-        self.inner.len()
-    }
-
-    fn len_resolutions(&self) -> Vec<usize> {
-        self.inner.len_resolutions()
-    }
-
-    #[getter]
-    fn get_is_empty(&self) -> bool {
-        self.inner.is_empty()
-    }
-
-    /// check if the stack contains the index or any of its parents
-    ///
-    /// This function is pretty inefficient.
-    fn contains(&self, h3index: u64) -> bool {
-        self.inner.contains(H3Cell::new(h3index))
-    }
-
-    fn compacted_indexes(&self) -> Py<PyArray1<u64>> {
-        let indexes: Vec<_> = self.inner.iter_compacted_cells().collect();
-        return_h3index_array(indexes)
-    }
-
-    fn compacted_indexes_at_resolution(&self, h3_resolution: u8) -> PyResult<Py<PyArray1<u64>>> {
-        check_valid_h3_resolution(h3_resolution).into_pyresult()?;
-        let cells = self
-            .inner
-            .get_compacted_cells_at_resolution(h3_resolution)
-            .to_vec();
-        Ok(return_h3index_array(cells))
-    }
-
-    fn uncompacted_indexes_at_resolution(&self, h3_resolution: u8) -> PyResult<Py<PyArray1<u64>>> {
-        check_valid_h3_resolution(h3_resolution).into_pyresult()?;
-        let cells = self
-            .inner
-            .iter_uncompacted_cells(h3_resolution)
-            .collect::<Result<Vec<_>, _>>()
-            .into_pyresult()?;
-        Ok(return_h3index_array(cells))
-    }
-}
-
-#[inline]
-fn return_h3index_array(cells: Vec<H3Cell>) -> Py<PyArray1<u64>> {
-    Python::with_gil(|py| cells_to_h3indexes(cells).into_pyarray(py).to_owned())
-}
diff --git a/src/error.rs b/src/error.rs
index 4dc9aaa..0cd7bb9 100644
--- a/src/error.rs
+++ b/src/error.rs
@@ -1,62 +1,85 @@
+use h3arrow::error::Error as A3Error;
 use pyo3::exceptions::{PyRuntimeError, PyValueError};
 use pyo3::{PyErr, PyResult};
+use rasterh3::Error;
 
 pub trait IntoPyResult<T> {
     fn into_pyresult(self) -> PyResult<T>;
 }
 
-trait IntoPyErr {
+pub trait IntoPyErr {
     fn into_pyerr(self) -> PyErr;
 }
 
-impl IntoPyErr for h3ron::Error {
+impl IntoPyErr for h3arrow::export::arrow2::error::Error {
+    fn into_pyerr(self) -> PyErr {
+        PyRuntimeError::new_err(self.to_string())
+    }
+}
+
+impl IntoPyErr for A3Error {
     fn into_pyerr(self) -> PyErr {
         match self {
-            Self::Domain
-            | Self::LatLonDomain
-            | Self::ResDomain
-            | Self::CellInvalid
-            | Self::DirectedEdgeInvalid
-            | Self::UndirectedEdgeInvalid
-            | Self::VertexInvalid
-            | Self::Pentagon
-            | Self::DuplicateInput
-            | Self::NotNeighbors
-            | Self::ResMismatch
-            | Self::MemoryBounds
-            | Self::OptionInvalid
-            | Self::DirectionInvalid(_) => PyValueError::new_err(self.to_string()),
-
-            Self::Failed
-            | Self::MemoryAlloc
-            | Self::UnknownError(_)
-            | Self::DecompressionError(_) => PyRuntimeError::new_err(self.to_string()),
+            A3Error::InvalidCellIndex(e) => e.into_pyerr(),
+            A3Error::InvalidVertexIndex(e) => e.into_pyerr(),
+            A3Error::InvalidDirectedEdgeIndex(e) => e.into_pyerr(),
+            A3Error::InvalidResolution(e) => e.into_pyerr(),
+            A3Error::InvalidLatLng(e) => e.into_pyerr(),
+            A3Error::InvalidGeometry(e) => e.into_pyerr(),
+            A3Error::CompactionError(e) => e.into_pyerr(),
+            A3Error::OutlinerError(e) => e.into_pyerr(),
+            A3Error::Arrow2(e) => e.into_pyerr(),
+            A3Error::NotAPrimitiveArrayU64
+            | A3Error::NonParsableCellIndex
+            | A3Error::InvalidWKB => PyValueError::new_err(self.to_string()),
         }
     }
 }
 
-impl<T> IntoPyResult<T> for Result<T, h3ron::Error> {
-    fn into_pyresult(self) -> PyResult<T> {
+macro_rules! impl_h3o_value_err {
+    ($($err_type:ty,)*) => {
+        $(
+            impl IntoPyErr for $err_type {
+                fn into_pyerr(self) -> PyErr {
+                    PyValueError::new_err(
+                        self.to_string()
+                    )
+                }
+            }
+        )*
+    }
+}
+
+impl_h3o_value_err!(
+    h3arrow::export::h3o::error::CompactionError,
+    h3arrow::export::h3o::error::InvalidCellIndex,
+    h3arrow::export::h3o::error::InvalidDirectedEdgeIndex,
+    h3arrow::export::h3o::error::InvalidGeometry,
+    h3arrow::export::h3o::error::InvalidLatLng,
+    h3arrow::export::h3o::error::InvalidResolution,
+    h3arrow::export::h3o::error::InvalidVertexIndex,
+    h3arrow::export::h3o::error::OutlinerError,
+);
+
+impl IntoPyErr for rasterh3::Error {
+    fn into_pyerr(self) -> PyErr {
         match self {
-            Ok(v) => Ok(v),
-            Err(err) => Err(err.into_pyerr()),
+            Error::TransformNotInvertible | Error::EmptyArray => {
+                PyValueError::new_err(self.to_string())
+            }
+            Error::InvalidLatLng(e) => e.into_pyerr(),
+            Error::InvalidGeometry(e) => e.into_pyerr(),
+            Error::InvalidResolution(e) => e.into_pyerr(),
+            Error::CompactionError(e) => e.into_pyerr(),
         }
     }
 }
 
-impl<T> IntoPyResult<T> for Result<T, h3ron_ndarray::Error> {
+impl<T, E> IntoPyResult<T> for Result<T, E>
+where
+    E: IntoPyErr,
+{
     fn into_pyresult(self) -> PyResult<T> {
-        match self {
-            Ok(v) => Ok(v),
-            Err(err) => match err {
-                h3ron_ndarray::Error::EmptyArray | h3ron_ndarray::Error::UnsupportedArrayShape => {
-                    Err(PyValueError::new_err(err.to_string()))
-                }
-                h3ron_ndarray::Error::TransformNotInvertible => {
-                    Err(PyRuntimeError::new_err(err.to_string()))
-                }
-                h3ron_ndarray::Error::H3ron(h3ron_e) => Err(h3ron_e.into_pyerr()),
-            },
-        }
+        self.map_err(IntoPyErr::into_pyerr)
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index 5e759a3..22224e2 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -8,41 +8,43 @@
     nonstandard_style
 )]
 
-use pyo3::{prelude::*, wrap_pyfunction, PyNativeType, Python};
+use pyo3::{prelude::*, wrap_pyfunction, Python};
 
 use crate::op::init_op_submodule;
 use crate::raster::init_raster_submodule;
 use crate::vector::init_vector_submodule;
-use crate::{collections::H3CompactedVec, polygon::Polygon};
-use h3ron::{H3Cell, Index};
 
-mod collections;
+mod arrow_interop;
 mod error;
 mod op;
-mod polygon;
 mod raster;
 mod transform;
 mod vector;
 
+pub(crate) const DEFAULT_CELL_COLUMN_NAME: &str = "cell";
+
 /// version of the module
 #[pyfunction]
 fn version() -> String {
     env!("CARGO_PKG_VERSION").to_string()
 }
 
-/// h3ron python bindings
+/// indicates if this extension has been compiled in release-mode
+#[pyfunction]
+fn is_release_build() -> bool {
+    #[cfg(debug_assertions)]
+    return false;
+
+    #[cfg(not(debug_assertions))]
+    return true;
+}
+
 #[pymodule]
 fn h3ronpy(py: Python<'_>, m: &PyModule) -> PyResult<()> {
     env_logger::init(); // run with the environment variable RUST_LOG set to "debug" for log output
 
-    m.add("H3CompactedVec", m.py().get_type::<H3CompactedVec>())?;
-    m.add("Polygon", m.py().get_type::<Polygon>())?;
-
     m.add_function(wrap_pyfunction!(version, m)?)?;
-
-    let vector_submod = PyModule::new(py, "vector")?;
-    init_vector_submodule(vector_submod)?;
-    m.add_submodule(vector_submod)?;
+    m.add_function(wrap_pyfunction!(is_release_build, m)?)?;
 
     let raster_submod = PyModule::new(py, "raster")?;
     init_raster_submodule(raster_submod)?;
@@ -52,9 +54,11 @@ fn h3ronpy(py: Python<'_>, m: &PyModule) -> PyResult<()> {
     init_op_submodule(op_submod)?;
     m.add_submodule(op_submod)?;
 
-    Ok(())
-}
+    let vector_submod = PyModule::new(py, "vector")?;
+    init_vector_submodule(vector_submod)?;
+    m.add_submodule(vector_submod)?;
 
-pub fn cells_to_h3indexes(cells: Vec<H3Cell>) -> Vec<u64> {
-    cells.into_iter().map(|cell| cell.h3index()).collect()
+    m.add("DEFAULT_CELL_COLUMN_NAME", DEFAULT_CELL_COLUMN_NAME)?;
+
+    Ok(())
 }
diff --git a/src/op.rs b/src/op.rs
deleted file mode 100644
index 3ed56f2..0000000
--- a/src/op.rs
+++ /dev/null
@@ -1,214 +0,0 @@
-use std::iter::once;
-use std::str::FromStr;
-
-use h3ron::collections::{HashMap, RandomState};
-use h3ron::error::check_valid_h3_resolution;
-use h3ron::iter::GridDiskBuilder;
-use h3ron::{H3Cell, Index};
-use ndarray::ArrayView1;
-use numpy::{PyArray1, PyReadonlyArray1, ToPyArray};
-use pyo3::exceptions::PyValueError;
-use pyo3::prelude::PyModule;
-use pyo3::{pyfunction, wrap_pyfunction, Py, PyErr, PyResult, Python};
-
-use crate::error::IntoPyResult;
-
-fn check_k_min_max(k_min: u32, k_max: u32) -> PyResult<()> {
-    if k_max < k_min {
-        return Err(PyValueError::new_err(format!(
-            "k_min must be smaller or equal to k_max. (k_min={}, k_max={})",
-            k_min, k_max
-        )));
-    }
-    Ok(())
-}
-
-#[inline]
-fn k_min_max_capacity(num_h3indexes: usize, k_min: u32, k_max: u32) -> PyResult<usize> {
-    Ok(num_h3indexes
-        * (h3ron::max_grid_disk_size(k_max).into_pyresult()?
-            - h3ron::max_grid_disk_size(k_min.saturating_sub(1)).into_pyresult()?))
-}
-
-#[allow(clippy::type_complexity)]
-#[pyfunction]
-fn grid_disk_distances(
-    h3index_arr: PyReadonlyArray1<u64>,
-    k_min: u32,
-    k_max: u32,
-) -> PyResult<(Py<PyArray1<u64>>, Py<PyArray1<u64>>, Py<PyArray1<u32>>)> {
-    check_k_min_max(k_min, k_max)?;
-    let h3indexes = h3index_arr.as_array();
-
-    let capacity = k_min_max_capacity(h3indexes.len(), k_min, k_max)?;
-    let mut center_h3indexes = Vec::with_capacity(capacity);
-    let mut ring_h3indexes = Vec::with_capacity(capacity);
-    let mut ks = Vec::with_capacity(capacity);
-
-    let mut builder = GridDiskBuilder::create(k_min, k_max).into_pyresult()?;
-
-    for h3index in h3indexes.iter() {
-        let cell = H3Cell::try_from(*h3index).into_pyresult()?;
-        builder.build_grid_disk(&cell).into_pyresult()?;
-        for (ring_cell, ring_cell_k) in &mut builder {
-            center_h3indexes.push(*h3index);
-            ring_h3indexes.push(ring_cell.h3index());
-            ks.push(ring_cell_k);
-        }
-    }
-
-    Ok(Python::with_gil(|py| {
-        (
-            center_h3indexes.to_pyarray(py).to_owned(),
-            ring_h3indexes.to_pyarray(py).to_owned(),
-            ks.to_pyarray(py).to_owned(),
-        )
-    }))
-}
-
-enum KAggregationMode {
-    Min,
-    Max,
-}
-
-impl FromStr for KAggregationMode {
-    type Err = PyErr;
-
-    fn from_str(s: &str) -> Result<Self, Self::Err> {
-        match s.to_lowercase().as_str() {
-            "min" => Ok(Self::Min),
-            "max" => Ok(Self::Max),
-            _ => Err(PyValueError::new_err("unknown way to aggregate k")),
-        }
-    }
-}
-
-#[allow(clippy::type_complexity)]
-fn kring_distances_agg_internal<A: Fn(&mut u32, u32)>(
-    h3indexes: ArrayView1<u64>,
-    k_min: u32,
-    k_max: u32,
-    agg_closure: A,
-) -> PyResult<(Py<PyArray1<u64>>, Py<PyArray1<u32>>)> {
-    check_k_min_max(k_min, k_max)?;
-
-    let mut cellmap = HashMap::with_capacity_and_hasher(
-        k_min_max_capacity(h3indexes.len(), k_min, k_max)?,
-        RandomState::default(),
-    );
-    let mut builder = GridDiskBuilder::create(k_min, k_max).into_pyresult()?;
-
-    for h3index in h3indexes.iter() {
-        let cell = H3Cell::try_from(*h3index).into_pyresult()?;
-        builder.build_grid_disk(&cell).into_pyresult()?;
-        for (ring_cell, ring_cell_k) in &mut builder {
-            cellmap
-                .entry(ring_cell.h3index())
-                .and_modify(|existing_ring_cell_k| agg_closure(existing_ring_cell_k, ring_cell_k))
-                .or_insert(ring_cell_k);
-        }
-    }
-
-    let capacity = cellmap.len();
-    let (h3indexes_out, k_out) = cellmap.into_iter().fold(
-        (Vec::with_capacity(capacity), Vec::with_capacity(capacity)),
-        |mut vecs, (h3index, k)| {
-            vecs.0.push(h3index);
-            vecs.1.push(k);
-            vecs
-        },
-    );
-    Ok(Python::with_gil(|py| {
-        (
-            h3indexes_out.to_pyarray(py).to_owned(),
-            k_out.to_pyarray(py).to_owned(),
-        )
-    }))
-}
-
-/// Vectorized grid-disk building, with the k-values of the rings being aggregated to their `min` or
-/// `max` value for each cell.
-#[allow(clippy::type_complexity)]
-#[pyfunction]
-fn grid_disk_distances_agg(
-    h3index_arr: PyReadonlyArray1<u64>,
-    k_min: u32,
-    k_max: u32,
-    aggregation_mode_str: &str,
-) -> PyResult<(Py<PyArray1<u64>>, Py<PyArray1<u32>>)> {
-    let h3indexes = h3index_arr.as_array();
-    let agg_closure = match KAggregationMode::from_str(aggregation_mode_str)? {
-        KAggregationMode::Min => |existing_k: &mut u32, new_k: u32| {
-            if *existing_k > new_k {
-                *existing_k = new_k
-            }
-        },
-        KAggregationMode::Max => |existing_k: &mut u32, new_k: u32| {
-            if *existing_k < new_k {
-                *existing_k = new_k
-            }
-        },
-    };
-    kring_distances_agg_internal(h3indexes, k_min, k_max, agg_closure)
-}
-
-#[pyfunction]
-fn change_resolution(
-    h3index_arr: PyReadonlyArray1<u64>,
-    h3_resolution: u8,
-) -> PyResult<Py<PyArray1<u64>>> {
-    check_valid_h3_resolution(h3_resolution).into_pyresult()?;
-
-    let mut out_vec = Vec::with_capacity(h3index_arr.len());
-    for h3index in h3index_arr.as_array().iter() {
-        // TODO: this will fail once the index-mode is checked
-        //       and the function gets applied to edges or other types.
-        let cell = H3Cell::try_from(*h3index).into_pyresult()?;
-        let change_iter = h3ron::iter::change_resolution(once(cell), h3_resolution);
-        out_vec.reserve(change_iter.size_hint().0);
-        for maybe_cell in change_iter {
-            out_vec.push(maybe_cell.into_pyresult()?.h3index())
-        }
-    }
-    Ok(Python::with_gil(|py| out_vec.to_pyarray(py).to_owned()))
-}
-
-#[allow(clippy::type_complexity)]
-#[pyfunction]
-fn change_resolution_paired(
-    h3index_arr: PyReadonlyArray1<u64>,
-    h3_resolution: u8,
-) -> PyResult<(Py<PyArray1<u64>>, Py<PyArray1<u64>>)> {
-    check_valid_h3_resolution(h3_resolution).into_pyresult()?;
-
-    let mut out_vec_after = Vec::with_capacity(h3index_arr.len());
-    let mut out_vec_before = Vec::with_capacity(h3index_arr.len());
-    for h3index in h3index_arr.as_array().iter() {
-        // TODO: this will fail once the index-mode is checked
-        //       and the function gets applied to edges or other types.
-        let cell = H3Cell::try_from(*h3index).into_pyresult()?;
-        let iter = h3ron::iter::change_resolution_tuple(once(cell), h3_resolution);
-        out_vec_after.reserve(iter.size_hint().0);
-        out_vec_before.reserve(iter.size_hint().0);
-        for pair in iter {
-            let (before_cell, after_cell) = pair.into_pyresult()?;
-            out_vec_before.push(before_cell.h3index());
-            out_vec_after.push(after_cell.h3index());
-        }
-    }
-    Ok(Python::with_gil(|py| {
-        (
-            out_vec_before.to_pyarray(py).to_owned(),
-            out_vec_after.to_pyarray(py).to_owned(),
-        )
-    }))
-}
-
-pub fn init_op_submodule(m: &PyModule) -> PyResult<()> {
-    m.add_function(wrap_pyfunction!(grid_disk_distances, m)?)?;
-    m.add_function(wrap_pyfunction!(grid_disk_distances_agg, m)?)?;
-    m.add_function(wrap_pyfunction!(change_resolution, m)?)?;
-    m.add_function(wrap_pyfunction!(change_resolution_paired, m)?)?;
-
-    Ok(())
-}
diff --git a/src/op/compact.rs b/src/op/compact.rs
new file mode 100644
index 0000000..6783d29
--- /dev/null
+++ b/src/op/compact.rs
@@ -0,0 +1,28 @@
+use h3arrow::algorithm::CompactOp;
+use h3arrow::export::h3o::Resolution;
+use pyo3::prelude::*;
+
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, mixed_resolutions = false))]
+pub(crate) fn compact(cellarray: &PyAny, mixed_resolutions: bool) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let compacted = if mixed_resolutions {
+        cellindexarray.compact_mixed_resolutions()
+    } else {
+        cellindexarray.compact()
+    }
+    .into_pyresult()?;
+
+    with_pyarrow(|py, pyarrow| native_to_pyarray(compacted.into_inner().boxed(), py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, target_resolution))]
+pub(crate) fn uncompact(cellarray: &PyAny, target_resolution: u8) -> PyResult<PyObject> {
+    let target_resolution = Resolution::try_from(target_resolution).into_pyresult()?;
+    let out = pyarray_to_cellindexarray(cellarray)?.uncompact(target_resolution);
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out.into_inner().boxed(), py, pyarrow))
+}
diff --git a/src/op/measure.rs b/src/op/measure.rs
new file mode 100644
index 0000000..f6037a7
--- /dev/null
+++ b/src/op/measure.rs
@@ -0,0 +1,23 @@
+use crate::arrow_interop::{native_to_pyarray, pyarray_to_cellindexarray, with_pyarrow};
+use pyo3::prelude::*;
+
+#[pyfunction]
+#[pyo3(signature = (cellarray))]
+pub(crate) fn cells_area_m2(cellarray: &PyAny) -> PyResult<PyObject> {
+    let out = pyarray_to_cellindexarray(cellarray)?.area_m2();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out.boxed(), py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray))]
+pub(crate) fn cells_area_km2(cellarray: &PyAny) -> PyResult<PyObject> {
+    let out = pyarray_to_cellindexarray(cellarray)?.area_km2();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out.boxed(), py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray))]
+pub(crate) fn cells_area_rads2(cellarray: &PyAny) -> PyResult<PyObject> {
+    let out = pyarray_to_cellindexarray(cellarray)?.area_rads2();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out.boxed(), py, pyarrow))
+}
diff --git a/src/op/mod.rs b/src/op/mod.rs
new file mode 100644
index 0000000..bdd2516
--- /dev/null
+++ b/src/op/mod.rs
@@ -0,0 +1,29 @@
+use pyo3::prelude::*;
+
+mod compact;
+mod measure;
+mod neighbor;
+mod parse;
+mod resolution;
+mod valid;
+
+pub fn init_op_submodule(m: &PyModule) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(resolution::change_resolution, m)?)?;
+    m.add_function(wrap_pyfunction!(resolution::change_resolution_paired, m)?)?;
+    m.add_function(wrap_pyfunction!(resolution::cells_resolution, m)?)?;
+    m.add_function(wrap_pyfunction!(neighbor::grid_disk, m)?)?;
+    m.add_function(wrap_pyfunction!(neighbor::grid_disk_distances, m)?)?;
+    m.add_function(wrap_pyfunction!(neighbor::grid_ring_distances, m)?)?;
+    m.add_function(wrap_pyfunction!(neighbor::grid_disk_aggregate_k, m)?)?;
+    m.add_function(wrap_pyfunction!(parse::cells_parse, m)?)?;
+    m.add_function(wrap_pyfunction!(compact::compact, m)?)?;
+    m.add_function(wrap_pyfunction!(compact::uncompact, m)?)?;
+    m.add_function(wrap_pyfunction!(valid::cells_valid, m)?)?;
+    m.add_function(wrap_pyfunction!(valid::vertexes_valid, m)?)?;
+    m.add_function(wrap_pyfunction!(valid::directededges_valid, m)?)?;
+    m.add_function(wrap_pyfunction!(measure::cells_area_m2, m)?)?;
+    m.add_function(wrap_pyfunction!(measure::cells_area_km2, m)?)?;
+    m.add_function(wrap_pyfunction!(measure::cells_area_rads2, m)?)?;
+
+    Ok(())
+}
diff --git a/src/op/neighbor.rs b/src/op/neighbor.rs
new file mode 100644
index 0000000..9870fce
--- /dev/null
+++ b/src/op/neighbor.rs
@@ -0,0 +1,132 @@
+use h3arrow::algorithm::{GridDiskDistances, GridOp, KAggregationMethod};
+use h3arrow::array::CellIndexArray;
+use h3arrow::export::arrow2::array::{Array, PrimitiveArray};
+use pyo3::exceptions::{PyRuntimeError, PyValueError};
+use pyo3::{PyObject, PyResult};
+use std::str::FromStr;
+
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
+use crate::DEFAULT_CELL_COLUMN_NAME;
+use pyo3::prelude::*;
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, k, flatten = false))]
+pub(crate) fn grid_disk(cellarray: &PyAny, k: u32, flatten: bool) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let listarray = cellindexarray.grid_disk(k).into_pyresult()?;
+    if flatten {
+        let cellindexarray = listarray.into_flattened().into_pyresult()?;
+        with_pyarrow(|py, pyarrow| {
+            native_to_pyarray(cellindexarray.into_inner().boxed(), py, pyarrow)
+        })
+    } else {
+        with_pyarrow(|py, pyarrow| native_to_pyarray(listarray.into_inner().boxed(), py, pyarrow))
+    }
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, k, flatten = false))]
+pub(crate) fn grid_disk_distances(cellarray: &PyAny, k: u32, flatten: bool) -> PyResult<PyObject> {
+    let griddiskdistances = pyarray_to_cellindexarray(cellarray)?
+        .grid_disk_distances(k)
+        .into_pyresult()?;
+
+    return_griddiskdistances_table(griddiskdistances, flatten)
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, k_min, k_max, flatten = false))]
+pub(crate) fn grid_ring_distances(
+    cellarray: &PyAny,
+    k_min: u32,
+    k_max: u32,
+    flatten: bool,
+) -> PyResult<PyObject> {
+    if k_min >= k_max {
+        return Err(PyValueError::new_err("k_min must be less than k_max"));
+    }
+    let griddiskdistances = pyarray_to_cellindexarray(cellarray)?
+        .grid_ring_distances(k_min, k_max)
+        .into_pyresult()?;
+
+    return_griddiskdistances_table(griddiskdistances, flatten)
+}
+
+fn return_griddiskdistances_table(
+    griddiskdistances: GridDiskDistances<CellIndexArray>,
+    flatten: bool,
+) -> PyResult<PyObject> {
+    let (cells, distances) = if flatten {
+        (
+            griddiskdistances
+                .cells
+                .into_flattened()
+                .into_pyresult()?
+                .into_inner()
+                .boxed(),
+            griddiskdistances
+                .distances
+                .values()
+                .as_any()
+                .downcast_ref::<PrimitiveArray<u32>>()
+                .ok_or_else(|| PyRuntimeError::new_err("expected primitivearray<u32>"))
+                .map(|pa| pa.clone().to_boxed())?,
+        )
+    } else {
+        (
+            griddiskdistances.cells.into_inner().boxed(),
+            griddiskdistances.distances.boxed(),
+        )
+    };
+
+    with_pyarrow(|py, pyarrow| {
+        let arrays = [
+            native_to_pyarray(cells, py, pyarrow)?,
+            native_to_pyarray(distances, py, pyarrow)?,
+        ];
+        let table = pyarrow
+            .getattr("Table")?
+            .call_method1("from_arrays", (arrays, [DEFAULT_CELL_COLUMN_NAME, "k"]))?;
+        Ok(table.to_object(py))
+    })
+}
+
+struct KAggregationMethodWrapper(KAggregationMethod);
+
+impl FromStr for KAggregationMethodWrapper {
+    type Err = PyErr;
+
+    fn from_str(s: &str) -> Result<Self, Self::Err> {
+        match s.to_lowercase().as_str() {
+            "min" => Ok(Self(KAggregationMethod::Min)),
+            "max" => Ok(Self(KAggregationMethod::Max)),
+            _ => Err(PyValueError::new_err("unknown way to aggregate k")),
+        }
+    }
+}
+
+#[pyfunction]
+#[pyo3(signature = (cellarray, k, aggregation_method))]
+pub(crate) fn grid_disk_aggregate_k(
+    cellarray: &PyAny,
+    k: u32,
+    aggregation_method: &str,
+) -> PyResult<PyObject> {
+    let aggregation_method = KAggregationMethodWrapper::from_str(aggregation_method)?;
+
+    let griddiskaggk = pyarray_to_cellindexarray(cellarray)?
+        .grid_disk_aggregate_k(k, aggregation_method.0)
+        .into_pyresult()?;
+
+    with_pyarrow(|py, pyarrow| {
+        let arrays = [
+            native_to_pyarray(griddiskaggk.cells.into_inner().to_boxed(), py, pyarrow)?,
+            native_to_pyarray(griddiskaggk.distances.to_boxed(), py, pyarrow)?,
+        ];
+        let table = pyarrow
+            .getattr("Table")?
+            .call_method1("from_arrays", (arrays, [DEFAULT_CELL_COLUMN_NAME, "k"]))?;
+        Ok(table.to_object(py))
+    })
+}
diff --git a/src/op/parse.rs b/src/op/parse.rs
new file mode 100644
index 0000000..135e024
--- /dev/null
+++ b/src/op/parse.rs
@@ -0,0 +1,17 @@
+use h3arrow::algorithm::ParseCellsOp;
+use h3arrow::export::arrow2::array::Utf8Array;
+use pyo3::prelude::*;
+
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
+
+#[pyfunction]
+#[pyo3(signature = (utf8array, set_failing_to_invalid = false))]
+pub(crate) fn cells_parse(utf8array: &PyAny, set_failing_to_invalid: bool) -> PyResult<PyObject> {
+    let utf8array: Utf8Array<i32> = pyarray_to_native(utf8array)?;
+    let cells = utf8array
+        .parse_cells(set_failing_to_invalid)
+        .into_pyresult()?;
+
+    with_pyarrow(|py, pyarrow| native_to_pyarray(cells.into_inner().boxed(), py, pyarrow))
+}
diff --git a/src/op/resolution.rs b/src/op/resolution.rs
new file mode 100644
index 0000000..1573799
--- /dev/null
+++ b/src/op/resolution.rs
@@ -0,0 +1,51 @@
+use h3arrow::algorithm::ChangeResolutionOp;
+use h3arrow::export::h3o::Resolution;
+use pyo3::prelude::*;
+
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
+use crate::DEFAULT_CELL_COLUMN_NAME;
+
+#[pyfunction]
+pub(crate) fn change_resolution(cellarray: &PyAny, h3_resolution: u8) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let h3_resolution = Resolution::try_from(h3_resolution).into_pyresult()?;
+    let out = cellindexarray
+        .change_resolution(h3_resolution)
+        .into_pyresult()?;
+
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out.into_inner().boxed(), py, pyarrow))
+}
+
+#[pyfunction]
+pub(crate) fn change_resolution_paired(cellarray: &PyAny, h3_resolution: u8) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let h3_resolution = Resolution::try_from(h3_resolution).into_pyresult()?;
+    let pair = cellindexarray
+        .change_resolution_paired(h3_resolution)
+        .into_pyresult()?;
+
+    with_pyarrow(|py, pyarrow| {
+        let arrays = [
+            native_to_pyarray(pair.before.into_inner().boxed(), py, pyarrow)?,
+            native_to_pyarray(pair.after.into_inner().boxed(), py, pyarrow)?,
+        ];
+        let table = pyarrow.getattr("Table")?.call_method1(
+            "from_arrays",
+            (
+                arrays,
+                [
+                    format!("{}_before", DEFAULT_CELL_COLUMN_NAME),
+                    format!("{}_after", DEFAULT_CELL_COLUMN_NAME),
+                ],
+            ),
+        )?;
+        Ok(table.to_object(py))
+    })
+}
+
+#[pyfunction]
+pub(crate) fn cells_resolution(cellarray: &PyAny) -> PyResult<PyObject> {
+    let resarray = pyarray_to_cellindexarray(cellarray)?.resolution();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(resarray.into_inner().boxed(), py, pyarrow))
+}
diff --git a/src/op/valid.rs b/src/op/valid.rs
new file mode 100644
index 0000000..29b572f
--- /dev/null
+++ b/src/op/valid.rs
@@ -0,0 +1,51 @@
+use h3arrow::array::{
+    CellIndexArray, DirectedEdgeIndexArray, FromIteratorWithValidity, H3Array, VertexIndexArray,
+};
+use h3arrow::export::arrow2::array::{BooleanArray, PrimitiveArray};
+use h3arrow::export::arrow2::bitmap::{Bitmap, MutableBitmap};
+use h3arrow::export::arrow2::datatypes::DataType;
+use pyo3::prelude::*;
+
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
+
+fn h3index_valid<A>(arr: &PyAny, booleanarray: bool) -> PyResult<PyObject>
+where
+    A: FromIteratorWithValidity<Option<u64>> + Into<PrimitiveArray<u64>> + H3Array,
+{
+    let u64array = pyarray_to_uint64array(arr)?;
+    let validated = A::from_iter_with_validity(u64array.iter().map(|v| v.copied()));
+
+    with_pyarrow(|py, pyarrow| {
+        native_to_pyarray(
+            if booleanarray {
+                let pa = validated.into();
+                let bm: Bitmap = pa
+                    .validity()
+                    .cloned()
+                    .unwrap_or_else(|| MutableBitmap::from_len_set(pa.len()).into());
+                BooleanArray::try_new(DataType::Boolean, bm, None)
+                    .into_pyresult()?
+                    .boxed()
+            } else {
+                validated.into().boxed()
+            },
+            py,
+            pyarrow,
+        )
+    })
+}
+
+macro_rules! impl_h3index_valid {
+    ($name:ident, $arr_type:ty) => {
+        #[pyfunction]
+        #[pyo3(signature = (array, booleanarray = false))]
+        pub(crate) fn $name(array: &PyAny, booleanarray: bool) -> PyResult<PyObject> {
+            h3index_valid::<$arr_type>(array, booleanarray)
+        }
+    };
+}
+
+impl_h3index_valid!(cells_valid, CellIndexArray);
+impl_h3index_valid!(vertexes_valid, VertexIndexArray);
+impl_h3index_valid!(directededges_valid, DirectedEdgeIndexArray);
diff --git a/src/polygon.rs b/src/polygon.rs
deleted file mode 100644
index 403bed0..0000000
--- a/src/polygon.rs
+++ /dev/null
@@ -1,61 +0,0 @@
-use crate::error::IntoPyResult;
-use geo_types as gt;
-use h3ron::{H3Cell, Index, ToAlignedLinkedPolygons, ToLinkedPolygons};
-use numpy::PyReadonlyArray1;
-use py_geo_interface::to_py::AsGeoInterface;
-use pyo3::prelude::*;
-
-#[pyclass]
-pub struct Polygon {
-    inner: gt::Polygon<f64>,
-}
-
-#[pymethods]
-impl Polygon {
-    #[staticmethod]
-    #[pyo3(signature = (h3index_arr, smoothen=false))]
-    fn from_h3indexes(
-        h3index_arr: PyReadonlyArray1<u64>,
-        smoothen: bool,
-    ) -> PyResult<Vec<Polygon>> {
-        let h3indexes: Vec<_> = h3index_arr
-            .as_array()
-            .iter()
-            .map(|hi| H3Cell::new(*hi))
-            .collect();
-
-        Ok(h3indexes
-            .to_linked_polygons(smoothen)
-            .into_pyresult()?
-            .into_iter()
-            .map(|poly| Polygon { inner: poly })
-            .collect())
-    }
-
-    #[staticmethod]
-    #[pyo3(signature = (h3index_arr, align_to_h3_resolution, smoothen=false))]
-    fn from_h3indexes_aligned(
-        h3index_arr: PyReadonlyArray1<u64>,
-        align_to_h3_resolution: u8,
-        smoothen: bool,
-    ) -> PyResult<Vec<Polygon>> {
-        let h3indexes: Vec<_> = h3index_arr
-            .as_array()
-            .iter()
-            .map(|hi| H3Cell::new(*hi))
-            .collect();
-
-        Ok(h3indexes
-            .to_aligned_linked_polygons(align_to_h3_resolution, smoothen)
-            .into_pyresult()?
-            .into_iter()
-            .map(|poly| Polygon { inner: poly })
-            .collect())
-    }
-
-    // python __geo_interface__ spec: https://gist.github.com/sgillies/2217756
-    #[getter]
-    fn __geo_interface__(&self, py: Python) -> PyResult<PyObject> {
-        self.inner.as_geointerface_pyobject(py)
-    }
-}
diff --git a/src/raster.rs b/src/raster.rs
index 414d2e7..08fe663 100644
--- a/src/raster.rs
+++ b/src/raster.rs
@@ -2,21 +2,21 @@ use std::hash::Hash;
 use std::iter::repeat;
 use std::str::FromStr;
 
+use h3arrow::array::CellIndexArray;
+use h3arrow::export::arrow2::array::PrimitiveArray;
+use h3arrow::export::h3o::{CellIndex, Resolution};
 use ndarray::ArrayView2;
-use numpy::{IntoPyArray, Ix1, PyArray, PyReadonlyArray2};
+use numpy::PyReadonlyArray2;
 use ordered_float::OrderedFloat;
 use pyo3::exceptions::PyValueError;
 use pyo3::{prelude::*, wrap_pyfunction, PyNativeType};
 
-use h3ron::error::check_valid_h3_resolution;
-use h3ron_ndarray as h3n;
-
-use crate::cells_to_h3indexes;
+use crate::arrow_interop::{native_to_pyarray, with_pyarrow};
 use crate::error::IntoPyResult;
 use crate::transform::Transform;
 
 pub struct AxisOrder {
-    pub inner: h3n::AxisOrder,
+    pub inner: rasterh3::AxisOrder,
 }
 
 impl FromStr for AxisOrder {
@@ -25,10 +25,10 @@ impl FromStr for AxisOrder {
     fn from_str(s: &str) -> Result<Self, Self::Err> {
         match s {
             "yx" | "YX" => Ok(Self {
-                inner: h3n::AxisOrder::YX,
+                inner: rasterh3::AxisOrder::YX,
             }),
             "xy" | "XY" => Ok(Self {
-                inner: h3n::AxisOrder::XY,
+                inner: rasterh3::AxisOrder::XY,
             }),
             _ => Err(PyValueError::new_err("unknown axis order")),
         }
@@ -36,7 +36,7 @@ impl FromStr for AxisOrder {
 }
 
 pub struct ResolutionSearchMode {
-    pub inner: h3n::ResolutionSearchMode,
+    pub inner: rasterh3::ResolutionSearchMode,
 }
 
 impl FromStr for ResolutionSearchMode {
@@ -45,10 +45,10 @@ impl FromStr for ResolutionSearchMode {
     fn from_str(s: &str) -> Result<Self, Self::Err> {
         match s {
             "min_diff" | "min-diff" => Ok(Self {
-                inner: h3n::ResolutionSearchMode::MinDiff,
+                inner: rasterh3::ResolutionSearchMode::MinDiff,
             }),
             "smaller_than_pixel" | "smaller-than-pixel" => Ok(Self {
-                inner: h3n::ResolutionSearchMode::SmallerThanPixel,
+                inner: rasterh3::ResolutionSearchMode::SmallerThanPixel,
             }),
             _ => Err(PyValueError::new_err("unknown resolution search mode")),
         }
@@ -68,13 +68,11 @@ pub fn nearest_h3_resolution(
     let search_mode = ResolutionSearchMode::from_str(search_mode_str)?;
     let shape: [usize; 2] = shape_any.extract()?;
 
-    h3n::resolution::nearest_h3_resolution(
-        &shape,
-        &transform.inner,
-        &axis_order.inner,
-        search_mode.inner,
-    )
-    .into_pyresult()
+    search_mode
+        .inner
+        .nearest_h3_resolution(shape, &transform.inner, &axis_order.inner)
+        .into_pyresult()
+        .map(Into::into)
 }
 
 #[allow(clippy::type_complexity)]
@@ -84,36 +82,28 @@ fn raster_to_h3<'a, T>(
     nodata_value: &'a Option<T>,
     h3_resolution: u8,
     axis_order_str: &str,
-    compacted: bool,
-) -> PyResult<(Vec<T>, Vec<u64>)>
+    compact: bool,
+) -> PyResult<(Vec<T>, Vec<CellIndex>)>
 where
     T: PartialEq + Sized + Sync + Eq + Hash + Copy,
 {
     let axis_order = AxisOrder::from_str(axis_order_str)?;
-    check_valid_h3_resolution(h3_resolution).into_pyresult()?;
+    let h3_resolution = Resolution::try_from(h3_resolution).into_pyresult()?;
 
-    let conv = h3n::H3Converter::new(arr, nodata_value, &transform.inner, axis_order.inner);
+    let conv = rasterh3::H3Converter::new(arr, nodata_value, &transform.inner, axis_order.inner);
 
     let mut values = vec![];
     let mut cells = vec![];
-    for (value, compacted_vec) in conv.to_h3(h3_resolution, compacted).into_pyresult()? {
-        let num_cells = if compacted {
-            let len_before = cells.len();
-            cells.extend(compacted_vec.iter_compacted_cells());
-            cells.len() - len_before
+    for (value, cell_coverage) in conv.to_h3(h3_resolution, compact).into_pyresult()? {
+        let len_before = cells.len();
+        if compact {
+            cells.extend(cell_coverage.into_compacted_iter());
         } else {
-            let mut value_cells = compacted_vec
-                .iter_uncompacted_cells(h3_resolution)
-                .collect::<Result<Vec<_>, _>>()
-                .into_pyresult()?;
-            let num_cells = value_cells.len();
-            cells.append(&mut value_cells);
-            num_cells
+            cells.extend(cell_coverage.into_uncompacted_iter(h3_resolution));
         };
-        values.extend(repeat(*value).take(num_cells));
+        values.extend(repeat(*value).take(cells.len() - len_before));
     }
-
-    Ok((values, cells_to_h3indexes(cells)))
+    Ok((values, cells))
 }
 
 macro_rules! make_raster_to_h3_variant {
@@ -124,25 +114,30 @@ macro_rules! make_raster_to_h3_variant {
             transform: &Transform,
             h3_resolution: u8,
             axis_order_str: &str,
-            compacted: bool,
+            compact: bool,
             nodata_value: Option<$dtype>,
-        ) -> PyResult<(Py<PyArray<$dtype, Ix1>>, Py<PyArray<u64, Ix1>>)> {
+        ) -> PyResult<(PyObject, PyObject)> {
             let arr = np_array.as_array();
-            raster_to_h3(
+            let (values, cells) = raster_to_h3(
                 &arr,
                 transform,
                 &nodata_value,
                 h3_resolution,
                 axis_order_str,
-                compacted,
-            )
-            .map(|(values, h3indexes)| {
-                Python::with_gil(|py| {
-                    (
-                        values.into_pyarray(py).to_owned(),
-                        h3indexes.into_pyarray(py).to_owned(),
-                    )
-                })
+                compact,
+            )?;
+
+            with_pyarrow(|py, pyarrow| {
+                let values = PrimitiveArray::from_vec(values);
+                let cells = native_to_pyarray(
+                    CellIndexArray::from(cells).into_inner().boxed(),
+                    py,
+                    pyarrow,
+                )?;
+
+                let values = native_to_pyarray(values.boxed(), py, pyarrow)?;
+
+                Ok((values, cells))
             })
         }
     };
@@ -156,31 +151,35 @@ macro_rules! make_raster_to_h3_float_variant {
             transform: &Transform,
             h3_resolution: u8,
             axis_order_str: &str,
-            compacted: bool,
+            compact: bool,
             nodata_value: Option<$dtype>,
-        ) -> PyResult<(Py<PyArray<$dtype, Ix1>>, Py<PyArray<u64, Ix1>>)> {
-            {
-                let arr = np_array.as_array();
-                // create a copy with the values wrapped in ordered floats to
-                // support the internal hashing
-                let of_arr = arr.map(|v| OrderedFloat::from(*v));
-                raster_to_h3(
-                    &of_arr.view(),
-                    transform,
-                    &nodata_value.map(OrderedFloat::from),
-                    h3_resolution,
-                    axis_order_str,
-                    compacted,
-                )
-            }
-            .map(|(values, h3indexes)| {
-                let float_vec: Vec<_> = values.into_iter().map(|v| *v).collect();
-                Python::with_gil(|py| {
-                    (
-                        float_vec.into_pyarray(py).to_owned(),
-                        h3indexes.into_pyarray(py).to_owned(),
-                    )
-                })
+        ) -> PyResult<(PyObject, PyObject)> {
+            let arr = np_array.as_array();
+            // create a copy with the values wrapped in ordered floats to
+            // support the internal hashing
+            let of_arr = arr.map(|v| OrderedFloat::from(*v));
+            let (values, cells) = raster_to_h3(
+                &of_arr.view(),
+                transform,
+                &nodata_value.map(OrderedFloat::from),
+                h3_resolution,
+                axis_order_str,
+                compact,
+            )?;
+
+            with_pyarrow(|py, pyarrow| {
+                let values = PrimitiveArray::<$dtype>::from_vec(
+                    values.into_iter().map(|v| v.into_inner()).collect(),
+                );
+                let cells = native_to_pyarray(
+                    CellIndexArray::from(cells).into_inner().boxed(),
+                    py,
+                    pyarrow,
+                )?;
+
+                let values = native_to_pyarray(values.boxed(), py, pyarrow)?;
+
+                Ok((values, cells))
             })
         }
     };
diff --git a/src/transform.rs b/src/transform.rs
index 7dfa9d1..364f8e0 100644
--- a/src/transform.rs
+++ b/src/transform.rs
@@ -1,14 +1,12 @@
-use pyo3::prelude::*;
-
-use h3ron_ndarray as h3n;
 use pyo3::basic::CompareOp;
 use pyo3::exceptions::PyNotImplementedError;
+use pyo3::prelude::*;
 
 /// affine geotransform
 #[pyclass]
 #[derive(Clone)]
 pub struct Transform {
-    pub inner: h3n::Transform,
+    pub(crate) inner: rasterh3::Transform,
 }
 
 #[pymethods]
@@ -17,33 +15,34 @@ impl Transform {
     #[new]
     pub fn new(a: f64, b: f64, c: f64, d: f64, e: f64, f: f64) -> Self {
         Self {
-            inner: h3n::Transform::new(a, b, c, d, e, f),
+            inner: rasterh3::Transform::new(a, b, c, d, e, f),
         }
     }
 
     /// construct a Transform from a six-values array as used by GDAL
     #[staticmethod]
-    pub fn from_gdal(gdal_transform: &PyAny) -> PyResult<Self> {
-        let gt: [f64; 6] = gdal_transform.extract()?;
-        Ok(Transform {
-            inner: h3n::Transform::from_gdal(&gt),
-        })
+    pub fn from_gdal(gdal_transform: [f64; 6]) -> Self {
+        Transform {
+            inner: rasterh3::Transform::from_gdal(&gdal_transform),
+        }
     }
 
     /// construct a Transform from a six-values array as used by rasterio
     #[staticmethod]
-    pub fn from_rasterio(rio_transform: &PyAny) -> PyResult<Self> {
-        let rt: [f64; 6] = rio_transform.extract()?;
-        Ok(Transform {
-            inner: h3n::Transform::from_rasterio(&rt),
-        })
+    pub fn from_rasterio(rio_transform: [f64; 6]) -> Self {
+        Transform {
+            inner: rasterh3::Transform::from_rasterio(&rio_transform),
+        }
     }
 
     fn __richcmp__(&self, other: Transform, op: CompareOp) -> PyResult<bool> {
         match op {
             CompareOp::Eq => Ok(self.inner == other.inner),
             CompareOp::Ne => Ok(self.inner != other.inner),
-            _ => Err(PyNotImplementedError::new_err("not implemented")),
+            _ => Err(PyNotImplementedError::new_err(format!(
+                "{:?} is not implemented",
+                op
+            ))),
         }
     }
 
diff --git a/src/vector.rs b/src/vector.rs
index 2f5a3d6..82a9b4a 100644
--- a/src/vector.rs
+++ b/src/vector.rs
@@ -1,132 +1,268 @@
-use geo_types::Geometry;
-use std::io::Cursor;
+use geo::BoundingRect;
+use h3arrow::algorithm::ToCoordinatesOp;
+use h3arrow::array::from_geo::{ToCellIndexArray, ToCellListArray, ToCellsOptions};
+use h3arrow::array::to_geoarrow::{ToWKBLineStrings, ToWKBLines, ToWKBPoints, ToWKBPolygons};
+use h3arrow::array::CellIndexArray;
+use h3arrow::export::arrow2::array::{BinaryArray, Float64Array, ListArray};
+use h3arrow::export::arrow2::bitmap::Bitmap;
+use h3arrow::export::geoarrow::{GeometryArrayTrait, WKBArray};
+use h3arrow::export::h3o::geom::ToGeo;
+use h3arrow::export::h3o::Resolution;
+use itertools::multizip;
+use pyo3::prelude::*;
+use pyo3::types::PyTuple;
 
-use geozero::wkb::{FromWkb, WkbDialect};
-use numpy::{IntoPyArray, Ix1, PyArray, PyReadonlyArray1};
-use pyo3::exceptions::PyValueError;
-use pyo3::{prelude::*, wrap_pyfunction};
-use rayon::prelude::*;
+use crate::arrow_interop::*;
+use crate::error::IntoPyResult;
 
-use h3ron::{compact_cells, H3Cell, Index, ToH3Cells, ToIntersectingH3Cells};
-use ndarray::{Array1, Zip};
-use pyo3::types::PyBytes;
+#[pyfunction]
+#[pyo3(signature = (cellarray,))]
+pub(crate) fn cells_bounds(cellarray: &PyAny) -> PyResult<Option<PyObject>> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    if let Some(rect) = cellindexarray.bounding_rect() {
+        Python::with_gil(|py| {
+            Ok(Some(
+                PyTuple::new(py, [rect.min().x, rect.min().y, rect.max().x, rect.max().y])
+                    .to_object(py),
+            ))
+        })
+    } else {
+        Ok(None)
+    }
+}
 
-use crate::error::IntoPyResult;
+#[pyfunction]
+#[pyo3(signature = (cellarray,))]
+pub(crate) fn cells_bounds_arrays(cellarray: &PyAny) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let mut minx_vec = vec![0.0f64; cellindexarray.len()];
+    let mut miny_vec = vec![0.0f64; cellindexarray.len()];
+    let mut maxx_vec = vec![0.0f64; cellindexarray.len()];
+    let mut maxy_vec = vec![0.0f64; cellindexarray.len()];
+    let mut validity_vec = vec![false; cellindexarray.len()];
 
-fn geom_to_h3(
-    geom: &Geometry,
-    h3_resolution: u8,
-    do_compact: bool,
-    partial_intersecting: bool,
-) -> PyResult<Vec<u64>> {
-    let mut cells: Vec<H3Cell> = match (geom, partial_intersecting) {
-        (Geometry::Polygon(poly), true) => poly
-            .to_intersecting_h3_cells(h3_resolution)
-            .into_pyresult()?,
-        (Geometry::MultiPolygon(mpoly), true) => mpoly
-            .to_intersecting_h3_cells(h3_resolution)
-            .into_pyresult()?,
-        _ => geom
-            .to_h3_cells(h3_resolution)
-            .into_pyresult()?
-            .iter()
-            .collect(),
-    };
+    for (cell, minx, miny, maxx, maxy, validity) in multizip((
+        cellindexarray.iter(),
+        minx_vec.iter_mut(),
+        miny_vec.iter_mut(),
+        maxx_vec.iter_mut(),
+        maxy_vec.iter_mut(),
+        validity_vec.iter_mut(),
+    )) {
+        if let Some(cell) = cell {
+            if let Some(rect) = cell
+                .to_geom(true)
+                .ok()
+                .and_then(|poly| poly.bounding_rect())
+            {
+                *validity = true;
+                *minx = rect.min().x;
+                *miny = rect.min().y;
+                *maxx = rect.min().x;
+                *maxy = rect.max().y;
+            };
+        }
+    }
 
-    // deduplicate, in the case of overlaps or lines
-    cells.sort_unstable();
-    cells.dedup();
+    let validity_bm = Bitmap::from(validity_vec.as_slice());
 
-    let cells = if do_compact {
-        compact_cells(&cells)
-            .into_pyresult()?
-            .iter()
-            .map(|i| i.h3index())
-            .collect()
-    } else {
-        cells.into_iter().map(|i| i.h3index()).collect()
-    };
-    Ok(cells)
+    with_pyarrow(|py, pyarrow| {
+        let arrays = [
+            native_to_pyarray(
+                Float64Array::from_vec(minx_vec)
+                    .with_validity(Some(validity_bm.clone()))
+                    .boxed(),
+                py,
+                pyarrow,
+            )?,
+            native_to_pyarray(
+                Float64Array::from_vec(miny_vec)
+                    .with_validity(Some(validity_bm.clone()))
+                    .boxed(),
+                py,
+                pyarrow,
+            )?,
+            native_to_pyarray(
+                Float64Array::from_vec(maxx_vec)
+                    .with_validity(Some(validity_bm.clone()))
+                    .boxed(),
+                py,
+                pyarrow,
+            )?,
+            native_to_pyarray(
+                Float64Array::from_vec(maxy_vec)
+                    .with_validity(Some(validity_bm.clone()))
+                    .boxed(),
+                py,
+                pyarrow,
+            )?,
+        ];
+        let table = pyarrow
+            .getattr("Table")?
+            .call_method1("from_arrays", (arrays, ["minx", "miny", "maxx", "maxy"]))?;
+        Ok(table.to_object(py))
+    })
 }
 
-#[allow(clippy::type_complexity)]
 #[pyfunction]
-fn wkbbytes_with_ids_to_h3(
-    id_array: PyReadonlyArray1<u64>,
-    wkb_array: PyReadonlyArray1<PyObject>,
-    h3_resolution: u8,
-    do_compact: bool,
-    partial_intersecting: bool,
-) -> PyResult<(Py<PyArray<u64, Ix1>>, Py<PyArray<u64, Ix1>>)> {
-    // the solution with the argument typed as list of byte-instances is not great. This
-    // maybe can be improved with https://github.com/PyO3/rust-numpy/issues/175
-
-    if id_array.len() != wkb_array.len() {
-        return Err(PyValueError::new_err(
-            "input Ids and WKBs must be of the same length",
-        ));
+#[pyo3(signature = (cellarray, radians = false))]
+pub(crate) fn cells_to_coordinates(cellarray: &PyAny, radians: bool) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+
+    let coordinate_arrays = if radians {
+        cellindexarray.to_coordinates_radians()
+    } else {
+        cellindexarray.to_coordinates()
     }
+    .into_pyresult()?;
 
-    let geom_array: Array1<_> = extract_geometries(wkb_array)?.into();
+    with_pyarrow(|py, pyarrow| {
+        let arrays = [
+            native_to_pyarray(coordinate_arrays.lat.boxed(), py, pyarrow)?,
+            native_to_pyarray(coordinate_arrays.lng.boxed(), py, pyarrow)?,
+        ];
+        let table = pyarrow
+            .getattr("Table")?
+            .call_method1("from_arrays", (arrays, ["lat", "lng"]))?;
+        Ok(table.to_object(py))
+    })
+}
 
-    let out = Zip::from(id_array.as_array())
-        .and(geom_array.view())
-        .into_par_iter()
-        .map(|(id, geom)| {
-            geom_to_h3(geom, h3_resolution, do_compact, partial_intersecting)
-                .map(|h3indexes| (*id, h3indexes))
-        })
-        .try_fold(
-            || (vec![], vec![]),
-            |mut a, b| match b {
-                Ok((id, mut indexes)) => {
-                    for _ in 0..indexes.len() {
-                        a.0.push(id);
-                    }
-                    a.1.append(&mut indexes);
-                    Ok(a)
-                }
-                Err(err) => Err(err),
-            },
-        )
-        .try_reduce(
-            || (vec![], vec![]),
-            |mut a, mut b| {
-                b.0.append(&mut a.0);
-                b.1.append(&mut a.1);
-                Ok(b)
-            },
-        )?;
-
-    Ok(Python::with_gil(|py| {
-        (
-            out.0.into_pyarray(py).to_owned(),
-            out.1.into_pyarray(py).to_owned(),
+#[pyfunction]
+#[pyo3(signature = (cellarray, radians = false, link_cells = false))]
+pub(crate) fn cells_to_wkb_polygons(
+    cellarray: &PyAny,
+    radians: bool,
+    link_cells: bool,
+) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+
+    let cells = if link_cells {
+        WKBArray::from(
+            cellindexarray
+                .iter()
+                .flatten()
+                .to_geom(radians)
+                .into_pyresult()?
+                .0
+                .into_iter()
+                .map(|poly| Some(geo_types::Geometry::from(poly)))
+                .collect::<Vec<_>>(),
         )
-    }))
+    } else {
+        cellindexarray.to_wkb_polygons(!radians).unwrap()
+    }
+    .into_arrow()
+    .boxed();
+
+    with_pyarrow(|py, pyarrow| native_to_pyarray(cells, py, pyarrow))
 }
 
-pub fn init_vector_submodule(m: &PyModule) -> PyResult<()> {
-    m.add_function(wrap_pyfunction!(wkbbytes_with_ids_to_h3, m)?)?;
-    Ok(())
+#[pyfunction]
+#[pyo3(signature = (cellarray, radians = false))]
+pub(crate) fn cells_to_wkb_points(cellarray: &PyAny, radians: bool) -> PyResult<PyObject> {
+    let cellindexarray = pyarray_to_cellindexarray(cellarray)?;
+    let out = cellindexarray
+        .to_wkb_points(!radians)
+        .unwrap()
+        .into_arrow()
+        .boxed();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out, py, pyarrow))
 }
 
-fn extract_geometries(array: PyReadonlyArray1<PyObject>) -> PyResult<Vec<Geometry>> {
-    let array = array.as_array();
-    Python::with_gil(|py| {
-        array
-            .iter()
-            .map(|obj| match obj.extract::<&PyBytes>(py) {
-                // is WKB
-                Ok(pb) => {
-                    let mut cur = Cursor::new(pb.as_bytes());
-                    Geometry::from_wkb(&mut cur, WkbDialect::Wkb)
-                        .map_err(|e| PyValueError::new_err(format!("unable to parse WKB: {:?}", e)))
-                }
-
-                // is some kind ob object, trying to extract a geometry from it
-                Err(_) => obj.extract::<py_geo_interface::Geometry>(py).map(|gi| gi.0),
-            })
-            .collect::<PyResult<Vec<_>>>()
-    })
+#[pyfunction]
+#[pyo3(signature = (vertexarray, radians = false))]
+pub(crate) fn vertexes_to_wkb_points(vertexarray: &PyAny, radians: bool) -> PyResult<PyObject> {
+    let vertexindexarray = pyarray_to_vertexindexarray(vertexarray)?;
+    let out = vertexindexarray
+        .to_wkb_points(!radians)
+        .unwrap()
+        .into_arrow()
+        .boxed();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out, py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (array, radians = false))]
+pub(crate) fn directededges_to_wkb_lines(array: &PyAny, radians: bool) -> PyResult<PyObject> {
+    let array = pyarray_to_directededgeindexarray(array)?;
+    let out = array.to_wkb_lines(!radians).unwrap().into_arrow().boxed();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out, py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (array, radians = false))]
+pub(crate) fn directededges_to_wkb_linestrings(array: &PyAny, radians: bool) -> PyResult<PyObject> {
+    let array = pyarray_to_directededgeindexarray(array)?;
+    let out = array
+        .to_wkb_linestrings(!radians)
+        .unwrap()
+        .into_arrow()
+        .boxed();
+    with_pyarrow(|py, pyarrow| native_to_pyarray(out, py, pyarrow))
+}
+
+#[pyfunction]
+#[pyo3(signature = (array, resolution, compact = false, all_intersecting = true, flatten = false))]
+pub(crate) fn wkb_to_cells(
+    array: &PyAny,
+    resolution: u8,
+    compact: bool,
+    all_intersecting: bool,
+    flatten: bool,
+) -> PyResult<PyObject> {
+    let options = ToCellsOptions {
+        resolution: Resolution::try_from(resolution).into_pyresult()?,
+        compact,
+        all_intersecting,
+    };
+    let wkbarray = WKBArray::new(pyarray_to_native::<BinaryArray<i64>>(array)?);
+
+    if flatten {
+        let array = wkbarray
+            .to_cellindexarray(&options)
+            .into_pyresult()?
+            .into_inner();
+        with_pyarrow(|py, pyarrow| native_to_pyarray(array.boxed(), py, pyarrow))
+    } else {
+        let listarray: ListArray<_> = wkbarray
+            .to_celllistarray(&options)
+            .into_pyresult()?
+            .into_inner();
+        with_pyarrow(|py, pyarrow| native_to_pyarray(listarray.boxed(), py, pyarrow))
+    }
+}
+
+#[pyfunction]
+#[pyo3(signature = (obj, resolution, compact = false, all_intersecting = true))]
+pub(crate) fn geometry_to_cells(
+    obj: py_geo_interface::Geometry,
+    resolution: u8,
+    compact: bool,
+    all_intersecting: bool,
+) -> PyResult<PyObject> {
+    let options = ToCellsOptions {
+        resolution: Resolution::try_from(resolution).into_pyresult()?,
+        compact,
+        all_intersecting,
+    };
+
+    let cellindexarray = CellIndexArray::from(
+        h3arrow::array::from_geo::geometry_to_cells(&obj.0, &options).into_pyresult()?,
+    );
+    with_pyarrow(|py, pyarrow| native_to_pyarray(cellindexarray.into_inner().boxed(), py, pyarrow))
+}
+
+pub fn init_vector_submodule(m: &PyModule) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(cells_to_coordinates, m)?)?;
+    m.add_function(wrap_pyfunction!(cells_bounds, m)?)?;
+    m.add_function(wrap_pyfunction!(cells_bounds_arrays, m)?)?;
+    m.add_function(wrap_pyfunction!(cells_to_wkb_polygons, m)?)?;
+    m.add_function(wrap_pyfunction!(cells_to_wkb_points, m)?)?;
+    m.add_function(wrap_pyfunction!(vertexes_to_wkb_points, m)?)?;
+    m.add_function(wrap_pyfunction!(directededges_to_wkb_linestrings, m)?)?;
+    m.add_function(wrap_pyfunction!(directededges_to_wkb_lines, m)?)?;
+    m.add_function(wrap_pyfunction!(wkb_to_cells, m)?)?;
+    m.add_function(wrap_pyfunction!(geometry_to_cells, m)?)?;
+    Ok(())
 }


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
index 6dba15f..5802dcd 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -1,3 +1,9 @@
 from pathlib import Path
+import geopandas as gpd
 
 TESTDATA_PATH = Path(__file__).parent.parent / "data"
+
+
+def load_africa() -> gpd.GeoDataFrame:
+    world = gpd.read_file(TESTDATA_PATH / "naturalearth_110m_admin_0_countries.fgb")
+    return world[world["CONTINENT"] == "Africa"]
diff --git a/tests/arrow/__init__.py b/tests/arrow/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/arrow/test_vector.py b/tests/arrow/test_vector.py
new file mode 100644
index 0000000..bf825f9
--- /dev/null
+++ b/tests/arrow/test_vector.py
@@ -0,0 +1,11 @@
+from h3ronpy.arrow.vector import geometry_to_cells
+import pyarrow as pa
+import shapely
+
+
+def test_geometry_to_cells():
+    geom = shapely.Polygon(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)))
+    cells = geometry_to_cells(geom, 5, all_intersecting=True)
+    assert isinstance(cells, pa.Array)
+    assert cells.type == pa.uint64()
+    assert len(cells) > 10
diff --git a/tests/pandas/__init__.py b/tests/pandas/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_raster.py b/tests/pandas/test_raster.py
similarity index 77%
rename from tests/test_raster.py
rename to tests/pandas/test_raster.py
index f803527..07006ff 100644
--- a/tests/test_raster.py
+++ b/tests/pandas/test_raster.py
@@ -1,25 +1,27 @@
 try:
     import rasterio
+
     HAS_RASTERIO = True
 except ImportError:
-    # rasterio provides no pre-compiled windows wheels, so it is an
-    # optional dependency to simplify the CI setup of the github actions workflows
+    # rasterio is an optional dependency
     HAS_RASTERIO = False
 
 import numpy as np
 import pytest
 
-from h3ronpy.raster import raster_to_dataframe
+from h3ronpy.pandas.raster import raster_to_dataframe
+from h3ronpy import DEFAULT_CELL_COLUMN_NAME
+
+from tests import TESTDATA_PATH
 
-from . import TESTDATA_PATH
 
 @pytest.mark.skipif(not HAS_RASTERIO, reason="requires rasterio")
 def test_r_tiff():
     dataset = rasterio.open(TESTDATA_PATH / "r.tiff")
     band = dataset.read(1)
-    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0, compacted=True, geo=False)
+    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0, compact=True, geo=False)
     assert len(df) > 100
-    assert df.dtypes["h3index"] == "uint64"
+    assert df.dtypes[DEFAULT_CELL_COLUMN_NAME] == "uint64"
     assert df.dtypes["value"] == "uint8"
 
 
@@ -27,9 +29,9 @@ def test_r_tiff():
 def test_r_tiff_float32():
     dataset = rasterio.open(TESTDATA_PATH / "r.tiff")
     band = dataset.read(1).astype(np.float32)
-    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0.0, compacted=True, geo=False)
+    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0.0, compact=True, geo=False)
     assert len(df) > 100
-    assert df.dtypes["h3index"] == "uint64"
+    assert df.dtypes[DEFAULT_CELL_COLUMN_NAME] == "uint64"
     assert df.dtypes["value"] == "float32"
 
 
diff --git a/tests/pandas/test_resolution.py b/tests/pandas/test_resolution.py
new file mode 100644
index 0000000..bc44705
--- /dev/null
+++ b/tests/pandas/test_resolution.py
@@ -0,0 +1,45 @@
+from h3ronpy.pandas import change_resolution, change_resolution_paired, cells_resolution
+import numpy as np
+import math
+import h3.api.numpy_int as h3
+
+
+def test_change_resolution_up():
+    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 5), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
+    out_res = 9
+    changed = change_resolution(h3indexes, out_res)
+    assert changed.shape[0] == (int(math.pow(7, 4)) + 7)
+    for i in range(len(changed)):
+        assert h3.h3_get_resolution(changed[i]) == out_res
+
+
+def test_change_resolution_paired_up():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    out_res = 9
+    changed_df = change_resolution_paired(h3indexes, out_res)
+    assert len(changed_df) == 7
+    for i in range(len(changed_df)):
+        assert h3.h3_get_resolution(changed_df["cell_before"][i]) == 8
+        assert h3.h3_get_resolution(changed_df["cell_after"][i]) == out_res
+
+
+def test_change_resolution_down():
+    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 5), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
+    out_res = 4
+    changed = change_resolution(h3indexes, out_res)
+    assert changed.shape[0] == 2
+    assert h3.h3_get_resolution(changed[0]) == out_res
+    assert h3.h3_get_resolution(changed[1]) == out_res
+
+
+def test_cells_resolution():
+    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 5), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
+    res = cells_resolution(h3indexes)
+    assert len(res) == 2
+    assert res[0] == 5
+    assert res[1] == 8
diff --git a/tests/pandas/test_vector.py b/tests/pandas/test_vector.py
new file mode 100644
index 0000000..df5c73a
--- /dev/null
+++ b/tests/pandas/test_vector.py
@@ -0,0 +1,111 @@
+import shapely
+
+import pandas as pd
+from shapely.geometry import Point, GeometryCollection
+import pytest
+from h3ronpy.pandas import change_resolution
+from h3ronpy.pandas.vector import (
+    cells_to_points,
+    cells_to_polygons,
+    cells_dataframe_to_geodataframe,
+    geodataframe_to_cells,
+    geoseries_to_cells,
+)
+from h3ronpy import DEFAULT_CELL_COLUMN_NAME
+import geopandas as gpd
+from .. import load_africa
+
+
+def test_cells_to_points():
+    gs = cells_to_points(
+        [
+            0x8009FFFFFFFFFFF,
+        ]
+    )
+    assert isinstance(gs, gpd.GeoSeries)
+    assert gs.geom_type[0] == "Point"
+
+
+def test_cells_to_polygons():
+    cells = change_resolution(
+        [
+            0x8009FFFFFFFFFFF,
+        ],
+        3,
+    )
+    gs = cells_to_polygons(cells)
+    assert isinstance(gs, gpd.GeoSeries)
+    assert gs.geom_type[0] == "Polygon"
+    assert len(gs) == 286
+
+    linked_gs = cells_to_polygons(cells, link_cells=True)
+    assert isinstance(linked_gs, gpd.GeoSeries)
+    assert linked_gs.geom_type[0] == "Polygon"
+    assert len(linked_gs) == 1
+    assert shapely.get_num_coordinates(linked_gs[0]) > 120
+
+
+def test_cells_dataframe_to_geodataframe():
+    df = pd.DataFrame(
+        {
+            DEFAULT_CELL_COLUMN_NAME: [
+                0x8009FFFFFFFFFFF,
+            ],
+            "id": [
+                5,
+            ],
+        }
+    )
+    gdf = cells_dataframe_to_geodataframe(df)
+    assert isinstance(gdf, gpd.GeoDataFrame)
+    assert len(gdf) == len(df)
+    assert (gdf[DEFAULT_CELL_COLUMN_NAME] == df[DEFAULT_CELL_COLUMN_NAME]).all()
+    assert (gdf["id"] == df["id"]).all()
+    assert gdf.geometry.geom_type[0] == "Polygon"
+
+
+def test_cells_dataframe_to_geodataframe_empty():
+    # https://github.com/nmandery/h3ron/issues/17
+    df = pd.DataFrame({DEFAULT_CELL_COLUMN_NAME: []})
+    gdf = cells_dataframe_to_geodataframe(df)  # should not raise an ValueError.
+    assert gdf.empty
+
+
+def test_cells_geodataframe_to_cells():
+    africa = load_africa()
+    df = geodataframe_to_cells(africa, 4)
+    assert len(df) > len(africa)
+    assert df.dtypes[DEFAULT_CELL_COLUMN_NAME] == "uint64"
+
+
+def test_geoseries_to_cells_flatten():
+    africa = load_africa()
+    cells = geoseries_to_cells(africa.geometry, 4, flatten=True)
+    assert len(cells) >= len(africa)
+    assert cells.dtype == "uint64"
+
+
+def test_empty_geometrycollection_omitted():
+    gdf = gpd.GeoDataFrame(
+        {
+            "geometry": [
+                GeometryCollection(),
+            ]
+        },
+        crs="epsg:4326",
+    )
+    df = geodataframe_to_cells(gdf, 4)
+    assert len(df) == 0
+
+
+def test_fail_on_empty_point():
+    gdf = gpd.GeoDataFrame(
+        {
+            "geometry": [
+                Point(),
+            ]
+        },
+        crs="epsg:4326",
+    )
+    with pytest.raises(ValueError):
+        geodataframe_to_cells(gdf, 4)
diff --git a/tests/polars/__init__.py b/tests/polars/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/polars/test_compact.py b/tests/polars/test_compact.py
new file mode 100644
index 0000000..fbd9176
--- /dev/null
+++ b/tests/polars/test_compact.py
@@ -0,0 +1,46 @@
+import pytest
+
+from h3ronpy.polars import compact, change_resolution, uncompact
+import numpy as np
+import h3.api.numpy_int as h3
+
+
+def compact_to_one(expected_cell, input_cells, **kw):
+    compacted = compact(input_cells, **kw)
+    assert len(compacted) == 1
+    assert compacted[0] == expected_cell
+
+
+def test_compact():
+    cell = h3.geo_to_h3(10.3, 45.1, 8)
+    h3indexes = change_resolution(
+        np.array(
+            [
+                cell,
+            ],
+            dtype=np.uint64,
+        ),
+        10,
+    )
+    compact_to_one(cell, h3indexes)
+
+
+def test_compact_mixed_fail():
+    cell = h3.geo_to_h3(10.3, 45.1, 8)
+    with pytest.raises(ValueError, match="heterogen"):
+        compact_to_one(cell, [cell, h3.h3_to_parent(cell, 4)])
+
+
+def test_compact_mixed():
+    cell = h3.geo_to_h3(10.3, 45.1, 8)
+    compact_to_one(cell, [cell, h3.geo_to_h3(10.3, 45.1, 9)], mixed_resolutions=True)
+
+
+def test_uncompact():
+    cells = uncompact(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+        ],
+        9,
+    )
+    assert len(cells) == 7
diff --git a/tests/polars/test_coordinates.py b/tests/polars/test_coordinates.py
new file mode 100644
index 0000000..e390dbc
--- /dev/null
+++ b/tests/polars/test_coordinates.py
@@ -0,0 +1,50 @@
+from h3ronpy.polars.vector import cells_to_coordinates, cells_bounds, cells_bounds_arrays
+import polars as pl
+import numpy as np
+import h3.api.numpy_int as h3
+
+
+def test_cells_to_coordinates():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    coords = cells_to_coordinates(h3indexes)
+    assert len(coords) == 1
+    assert 10.0 < coords["lat"][0] < 11.0
+    assert 45.0 < coords["lng"][0] < 46.0
+
+
+def test_cells_bounds():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    bounds = cells_bounds(h3indexes)
+    assert bounds is not None
+    assert type(bounds) == tuple
+    assert len(bounds) == 4
+    assert bounds[0] < bounds[2]
+    assert bounds[1] < bounds[3]
+
+
+def test_cells_bounds_arrays():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(10.3, 45.1, 5),
+        ],
+        dtype=np.uint64,
+    )
+    bounds_df = cells_bounds_arrays(h3indexes)
+    assert bounds_df is not None
+    assert isinstance(bounds_df, pl.DataFrame)
+    assert len(bounds_df) == 2
+    assert "minx" in bounds_df
+    assert "maxx" in bounds_df
+    assert "miny" in bounds_df
+    assert "maxy" in bounds_df
diff --git a/tests/polars/test_measure.py b/tests/polars/test_measure.py
new file mode 100644
index 0000000..9e1b000
--- /dev/null
+++ b/tests/polars/test_measure.py
@@ -0,0 +1,21 @@
+import numpy as np
+import h3.api.numpy_int as h3
+from h3ronpy.polars import cells_area_km2
+import polars as pl
+
+
+def test_cells_area_km2():
+    cells = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(10.3, 45.1, 5),
+            h3.geo_to_h3(10.3, 45.1, 3),
+        ],
+        dtype=np.uint64,
+    )
+    areas = cells_area_km2(cells)
+    assert isinstance(areas, pl.Series)
+    assert len(areas) == 3
+    assert int(areas[0] * 100) == 62
+    assert int(areas[1]) == 213
+    assert int(areas[2]) == 10456
diff --git a/tests/polars/test_neighbor.py b/tests/polars/test_neighbor.py
new file mode 100644
index 0000000..dca835b
--- /dev/null
+++ b/tests/polars/test_neighbor.py
@@ -0,0 +1,89 @@
+from h3ronpy.polars import grid_disk, grid_disk_distances, grid_ring_distances, grid_disk_aggregate_k
+import numpy as np
+import h3.api.numpy_int as h3
+import polars as pl
+
+
+def test_grid_disk():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(5.3, -5.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    disks = grid_disk(h3indexes, 2)
+    assert len(disks) == 2
+    assert disks.dtype == pl.List(pl.UInt64())
+
+    disks_flat = grid_disk(h3indexes, 2, flatten=True)
+    assert len(disks_flat) > 20
+    assert disks_flat.dtype == pl.UInt64()
+
+
+def test_grid_disk_distances():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(5.3, -5.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    disks = grid_disk_distances(h3indexes, 2)
+    assert type(disks) == pl.DataFrame
+    assert len(disks) == len(h3indexes)
+    assert disks["cell"].dtype == pl.List(pl.UInt64())
+    assert disks["k"].dtype == pl.List(pl.UInt32())
+
+    centers = (
+        grid_disk_distances(h3indexes, 2, flatten=True)
+        .lazy()
+        .filter(pl.col("cell").is_in(pl.Series(h3indexes)))
+        .collect()
+    )
+    assert len(centers) == len(h3indexes)
+    assert len(centers["k"].unique()) == 1
+    assert centers["k"].unique()[0] == 0
+
+    # TODO: check values
+
+
+def test_grid_ring_distances():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(5.3, -5.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    disks = grid_ring_distances(h3indexes, 1, 2)
+    assert type(disks) == pl.DataFrame
+    assert len(disks) == len(h3indexes)
+    assert disks["cell"].dtype == pl.List(pl.UInt64())
+    assert disks["k"].dtype == pl.List(pl.UInt32())
+
+    centers = (
+        grid_ring_distances(h3indexes, 1, 2, flatten=True)
+        .lazy()
+        .filter(pl.col("cell").is_in(pl.Series(h3indexes)))
+        .collect()
+    )
+    assert len(centers) == 0
+
+    # TODO: check values
+
+
+def test_grid_disk_aggregate_k():
+    h3indexes = np.array(
+        [
+            h3.geo_to_h3(10.3, 45.1, 8),
+            h3.geo_to_h3(5.3, -5.1, 8),
+        ],
+        dtype=np.uint64,
+    )
+    disks = grid_disk_aggregate_k(h3indexes, 2, "max")
+    assert type(disks) == pl.DataFrame
+    assert disks["cell"].dtype == pl.UInt64()
+    assert disks["k"].dtype == pl.UInt32()
+
+    # TODO: check values
diff --git a/tests/polars/test_parse.py b/tests/polars/test_parse.py
new file mode 100644
index 0000000..783e6ff
--- /dev/null
+++ b/tests/polars/test_parse.py
@@ -0,0 +1,54 @@
+import pytest
+
+from h3ronpy.polars import cells_parse, cells_valid
+import numpy as np
+import h3.api.numpy_int as h3
+import polars as pl
+
+
+def test_cells_parse():
+    strings = np.array([h3.h3_to_string(h3.geo_to_h3(45.5, 10.2, 5)), "10.2, 45.5, 5"])
+    cells = cells_parse(strings)
+    assert len(cells) == 2
+    assert cells[0] == cells[1]
+
+
+def test_parse_cell_fail():
+    strings = np.array(
+        [
+            "invalid",
+        ]
+    )
+    with pytest.raises(ValueError, match="non-parsable CellIndex"):
+        cells_parse(strings)
+
+
+def test_parse_cell_set_invalid():
+    strings = np.array(
+        [
+            "invalid",
+        ]
+    )
+    cells = cells_parse(strings, set_failing_to_invalid=True)
+    assert len(cells) == 1
+    assert cells[0] is None
+
+
+def test_cells_valid():
+    input = np.array(
+        [45, h3.geo_to_h3(45.5, 10.2, 5)],
+        dtype=np.uint64,
+    )
+    cells = cells_valid(input, booleanarray=False)
+    assert len(cells) == 2
+    assert cells.dtype == pl.datatypes.UInt64()
+    assert cells[0] is None
+    assert cells[1] is not None
+
+    bools = cells_valid(input, booleanarray=True)
+    assert len(bools) == 2
+    assert bools.dtype == pl.datatypes.Boolean()
+    assert bools[0] is False
+    assert bools[1] is True
+
+    assert cells.is_not_null().eq(bools).all()
diff --git a/tests/polars/test_raster.py b/tests/polars/test_raster.py
new file mode 100644
index 0000000..577c4e0
--- /dev/null
+++ b/tests/polars/test_raster.py
@@ -0,0 +1,36 @@
+try:
+    import rasterio
+
+    HAS_RASTERIO = True
+except ImportError:
+    # rasterio is an optional dependency
+    HAS_RASTERIO = False
+
+import numpy as np
+import pytest
+import polars as pl
+
+from h3ronpy.polars.raster import raster_to_dataframe
+from h3ronpy import DEFAULT_CELL_COLUMN_NAME
+
+from tests import TESTDATA_PATH
+
+
+@pytest.mark.skipif(not HAS_RASTERIO, reason="requires rasterio")
+def test_r_tiff():
+    dataset = rasterio.open(TESTDATA_PATH / "r.tiff")
+    band = dataset.read(1)
+    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0, compact=True)
+    assert len(df) > 100
+    assert df[DEFAULT_CELL_COLUMN_NAME].dtype == pl.UInt64
+    assert df["value"].dtype == pl.UInt8
+
+
+@pytest.mark.skipif(not HAS_RASTERIO, reason="requires rasterio")
+def test_r_tiff_float32():
+    dataset = rasterio.open(TESTDATA_PATH / "r.tiff")
+    band = dataset.read(1).astype(np.float32)
+    df = raster_to_dataframe(band, dataset.transform, 8, nodata_value=0.0, compact=True)
+    assert len(df) > 100
+    assert df[DEFAULT_CELL_COLUMN_NAME].dtype == pl.UInt64
+    assert df["value"].dtype == pl.Float32
diff --git a/tests/test_op.py b/tests/test_op.py
deleted file mode 100644
index 4bedd7f..0000000
--- a/tests/test_op.py
+++ /dev/null
@@ -1,55 +0,0 @@
-import h3.api.numpy_int as h3
-import numpy as np
-import math
-
-from h3ronpy.op import grid_disk_distances_agg, grid_disk_distances, change_resolution, change_resolution_paired
-
-
-def test_kring_distances_agg():
-    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 8), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
-    k_max = 4
-    df = grid_disk_distances_agg(h3indexes, k_max, aggregation_method='min')
-    assert len(df) > 100
-    assert df['k'].min() == 0
-    assert df['k'].max() == k_max
-    assert len(np.unique(df["h3index"])) == len(df)
-
-
-def test_kring_distances():
-    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 8), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
-    k_max = 4
-    k_min = 2
-    df = grid_disk_distances(h3indexes, 4, k_min=k_min)
-    assert len(df) > 100
-    assert df['ring_k'].min() == k_min
-    assert df['ring_k'].max() == k_max
-    assert 'h3index' in df
-    assert 'ring_h3index' in df
-
-
-def test_change_resolution_down():
-    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 5), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
-    out_res = 4
-    changed = change_resolution(h3indexes, out_res)
-    assert changed.shape[0] == 2
-    assert h3.h3_get_resolution(changed[0]) == out_res
-    assert h3.h3_get_resolution(changed[1]) == out_res
-
-
-def test_change_resolution_up():
-    h3indexes = np.array([h3.geo_to_h3(10.2, 45.5, 5), h3.geo_to_h3(10.3, 45.1, 8)], dtype=np.uint64)
-    out_res = 9
-    changed = change_resolution(h3indexes, out_res)
-    assert changed.shape[0] == (int(math.pow(7, 4)) + 7)
-    for i in range(len(changed)):
-        assert h3.h3_get_resolution(changed[i]) == out_res
-
-
-def test_change_resolution_paired_up():
-    h3indexes = np.array([h3.geo_to_h3(10.3, 45.1, 8),], dtype=np.uint64)
-    out_res = 9
-    changed_df = change_resolution_paired(h3indexes, out_res)
-    assert len(changed_df) == 7
-    for i in range(len(changed_df)):
-        assert h3.h3_get_resolution(changed_df["h3index_before"][i]) == 8
-        assert h3.h3_get_resolution(changed_df["h3index_after"][i]) == out_res
diff --git a/tests/test_transform.py b/tests/test_transform.py
index 729b4e8..c935c08 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,4 +1,4 @@
-from h3ronpy.raster import Transform
+from h3ronpy.arrow.raster import Transform
 
 
 def test_transform_cmp():
diff --git a/tests/test_util.py b/tests/test_util.py
deleted file mode 100644
index 76f00ec..0000000
--- a/tests/test_util.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import pandas as pd
-
-from h3ronpy.util import dataframe_to_geodataframe
-
-
-def test_empty_to_geo():
-    # https://github.com/nmandery/h3ron/issues/17
-    df = pd.DataFrame({"h3index": []})
-    gdf = dataframe_to_geodataframe(df)  # should not raise an ValueError.
-    assert gdf.empty
diff --git a/tests/test_vector.py b/tests/test_vector.py
deleted file mode 100644
index 177665d..0000000
--- a/tests/test_vector.py
+++ /dev/null
@@ -1,24 +0,0 @@
-import geopandas as gpd
-from shapely.geometry import Point, GeometryCollection
-from h3ronpy.vector import geodataframe_to_h3
-import pytest
-
-
-def test_geodataframe_to_h3():
-    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
-    africa = world[world["continent"] == "Africa"]
-    df = geodataframe_to_h3(africa, 4)
-    assert len(df) > 1000
-    assert df.dtypes["h3index"] == "uint64"
-
-
-def test_empty_geometrycollection_omitted():
-    gdf = gpd.GeoDataFrame({"geometry": [GeometryCollection(),]}, crs="epsg:4326")
-    df = geodataframe_to_h3(gdf, 4)
-    assert len(df) == 0
-
-def test_fail_on_empty_point():
-    gdf = gpd.GeoDataFrame({"geometry": [Point(),]}, crs="epsg:4326")
-    with pytest.raises(Exception):
-        df = geodataframe_to_h3(gdf, 4)
-