Repository: skyportal/healpix-alchemy (language: python)

Issue:
Fast union of regions
We need to be able to rapidly calculate the area or probability contained within the union of a set of ZTF fields. Each ZTF field is represented as a Region, which in turn is represented as rows in the Tiles table.

A Tile is a multi-resolution HEALPix pixel. HEALPix is a tree data structure. At the lowest level (level 0) it has 12 nodes; at all nodes beneath that, it is a quadtree. A multi-resolution HEALPix pixel has several equivalent representations:

1. A level (a nonnegative integer) and a pixel number between 0 and 12 (2^level)^2 - 1 inclusive
2. An interval (start and end) of pixel indices at a common very high level in the tree (level 29 for 64-bit data)

The Tile table currently uses the second representation. If there are N tiles that are each stored as lists sorted by starting pixel index, then finding the union of tiles simply amounts to a merge of N sorted lists followed by a de-duplication step.

Similar algorithms are commonly used for manipulating time interval data (e.g. https://lscsoft.docs.ligo.org/ligo-segments/) and genomics (https://pyranges.readthedocs.io/en/latest/autoapi/pyranges/pyranges/index.html).

This algorithm is not that complicated to implement in Python, but it needs to be done within the database. I have struggled to figure out how to do this in basic SQL. A solution might involve [CTEs](https://www.postgresql.org/docs/current/queries-with.html), [window functions](https://www.postgresql.org/docs/current/functions-window.html), or [procedural SQL](https://www.postgresql.org/docs/current/plpgsql.html).
Factor cone search into separate package
The Cartesian indexing cone search strategy is not going to share any code at all with the HEALPix indexing. Factor the cone search code into a separate package:

https://github.com/skyportal/conesearch-alchemy


Discussion / hints:
PostgreSQL 14 has a new [multiranges](https://www.postgresql.org/docs/14/rangetypes.html#RANGETYPES-BUILTIN) type that may be worth investigating.
See sqlalchemy/sqlalchemy#7156 for multirange support in SQLAlchemy.


Fix patch:
diff --git a/README.md b/README.md
index 0f46a33..6328e79 100644
--- a/README.md
+++ b/README.md
@@ -1,16 +1,532 @@
 # HEALPix Alchemy
 
-The `healpix_alchemy` Python package extends [SQLAlchemy] will provide spatial
-indexing for astronomical sky coordinates, regions, and raster images (e.g.
-LIGO/Virgo and Fermi probability sky maps) in a relational database. It does
-not rely on any database extensions.
+The `healpix_alchemy` Python package is an extension for the [SQLAlchemy]
+object relational mapper. It adds region and image arithmetic to [PostgreSQL]
+(version 14 and newer) databases. It accelerates queries between point clouds,
+regions, and images (sometimes known in the geospatial community as rasters) by
+storing [multi-order] [HEALPix] indices in PostgreSQL's [range types].
 
-This package is a work in progress.
+The `healpix_alchemy` project is designed for astronomy applications,
+particularly for cross-matching galaxy catalogs, observation footprints, and
+all-sky images like [gravitational-wave probability sky maps] or even dust
+maps. However, it could be used in any context in which geometry is embedded on
+a unit sphere.
+
+`healpix_alchemy` is lean and minimalist because it leverages several existing
+projects: it consists of little more than a few lines of glue code to bind
+together [MOCPy], [SQLAlchemy], and PostgreSQL's [range types].
+
+`healpix_alchemy` serves a purpose similar to full-featured astronomy-focused
+spatial extensions like [Q3C], [H3C], and [pg_healpix], and geospatial
+extensions like [PgSphere] and [PostGIS]. What sets `healpix_alchemy` apart
+from these is that it is written in pure Python and requires no server-side
+database extensions. Consequently, `healpix_alchemy` can be used with managed
+PostgreSQL databases in the cloud like [Amazon RDS] and [Google Cloud SQL].
+
+[2MASS Redshift Survey]: https://lweb.cfa.harvard.edu/~dfabricant/huchra/2mass/
+[Aladin]: https://aladin.u-strasbg.fr
+[Amazon RDS]: https://aws.amazon.com/rds/
+[astropy.coordinates.SkyCoord]: https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord
+[astropy.units.Quantity]: https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity
+[custom column types]: https://docs.sqlalchemy.org/en/14/core/custom_types.html
+[Google Cloud SQL]: https://cloud.google.com/sql
+[Górski et al. (2005)]: https://doi.org/10.1086/427976
+[gravitational-wave probability sky maps]: https://emfollow.docs.ligo.org/userguide/tutorial/skymaps.html
+[GW200115_042309]: https://doi.org/10.3847/2041-8213/ac082e
+[H3C]: http://cdsarc.u-strasbg.fr/h3c
+[HEALPix]: https://healpix.sourceforge.io
+[hierarchical progressive surveys (HiPS)]: https://www.ivoa.net/documents/HiPS/
+[ICRS]: https://docs.astropy.org/en/stable/api/astropy.coordinates.builtin_frames.ICRS.html#astropy.coordinates.builtin_frames.ICRS
+[ligo-segments]: https://lscsoft.docs.ligo.org/ligo-segments/
+[MOCPy documentation]: https://cds-astro.github.io/mocpy/
+[MOCPy]: https://github.com/tboch/mocpy
+[multi-order coverage (MOC)]: https://ivoa.net/documents/MOC/
+[multi-order]: https://doi.org/10.1051/0004-6361/201526549
+[pg_healpix]: https://github.com/segasai/pg_healpix
+[PgSphere]: https://pgsphere.github.io
+[PostGIS]: https://postgis.net
+[PostgreSQL]: https://www.postgresql.org
+[pyranges]: https://github.com/biocore-ntnu/pyranges
+[Q3C]: https://github.com/segasai/q3c
+[range types]: https://www.postgresql.org/docs/current/rangetypes.html
+[S200115j]: https://gracedb.ligo.org/superevents/S200115j/view/
+[Singer & Price (2015)]: https://doi.org/10.1103/PhysRevD.93.024013
+[SQLAlchemy bulk insertion]: https://docs.sqlalchemy.org/en/14/orm/persistence_techniques.html#bulk-operations
+[SQLAlchemy]: https://www.sqlalchemy.org
+[UNIQ]: https://healpix.sourceforge.io/doc/html/intro_Geometric_Algebraic_Propert.htm#SECTION420
+[Zwicky Transient Facility]: https://www.ztf.caltech.edu
+
+## Theory of operation
+
+### HEALPix basics
+
+[HEALPix] is a scheme for subdividing and indexing the unit sphere, originally
+described by [Górski et al. (2005)]. Although it was originally designed for
+cosmic microwave background analysis, it has found many uses in astronomy,
+particularly through [multi-order coverage (MOC)] maps and
+[hierarchical progressive surveys (HiPS)] used extensively in the [Aladin]
+astronomical information system. It is also used by LIGO and Virgo to store and
+communicate [gravitational-wave probability sky maps].
+
+HEALPix can be thought of as a tree. At the lowest resolution, level 0, HEALPix
+subdivides the sphere into 12 equal-area base tiles, assigned integer indicies
+0 through 11. At level 1, each of the 12 base tiles is subdivided into 4 tiles.
+Every subsequent level divides each of the preceding level's tiles into 4 new
+tiles. At a given level, each of the base pixels has been divided into
+4<sup>_level_</sup> pixels (_nside_ = 2<sup>_level_</sup> pixels on each side).
+Thus there are _npix_ = 12×4<sup>_level_</sup> pixels at a given resolution,
+assigned integer indices from 0 through (_npix_-1). This is called the NESTED
+indexing scheme. (There is also a RING indexing scheme in which the indicies
+advance from east to west and then from north to south.)
+
+A HEALPix tile, a node in the HEALPix tree, is fully addressed by three pieces
+of information: the indexing scheme (RING or NESTED), the resolution level
+(_level_ or equivalently _nside_), and the pixel index (_ipix_, an integer
+between 0 and _npix_-1).
+
+The image below, reproduced from https://healpix.jpl.nasa.gov, illustrates the
+first 4 levels of refinement of a HEALPix grid.
+
+![The first four levels of HEALPix subdivision of the unit sphere](https://healpix.jpl.nasa.gov/images/healpixGridRefinement.jpg)
+
+### HEALPix interval sets
+
+A region on the sphere can be encoded by a collection of disjoint HEALPix
+tiles, potentially at a mix of different resolution levels. Typically, large
+low-resolution tiles are used on the interior of the region, and small
+high-resolution tiles are used on the boundary. This is called a
+[multi-order coverage (MOC)] map. An example, reproduced from the
+[MOCPy documentation], is shown below.
+
+![An example HEALPix multi-order coverage map](https://cds-astro.github.io/mocpy/_images/plot_SDSS_r.png)
+
+Much like MOCs, LIGO/Virgo/KAGRA gravitational-wave probability sky maps are
+stored as multi-resolution HEALPix data sets, but with a vector of
+floating-point values attached to each tile. A multi-order refinement mesh from
+an example sky map is shown below, reproduced from [Singer & Price (2015)].
+
+![Example multi-order refinement mesh from a gravitational-wave probability sky map](https://emfollow.docs.ligo.org/userguide/_images/healpix-adaptive-mesh.svg)
+
+[Reinecke & Hivon (2015)][multi-order] introduced HEALPix interval sets as an
+alternative encoding of MOCs that enables fast and simple unions,
+intersections, and queries. In an interval set, each HEALPix tile is described
+by the interval of pixel indices at some very high resolution
+_level_<sub>max</sub> that are descendents of that tile. In an interval set, a
+region is encoded as a disjoint collection of such intervals. A tile with a
+NESTED address given by (_level_, _ipix_) may be described as the half-open
+interval
+
+[_ipix_ 4<sup>_level_<sub>max</sub> - _level_</sup>,
+(_ipix_ + 1) 4<sup>_level_<sub>max</sub> - _level_</sup>).
+
+We use _level_<sub>max</sub> = 29 because this is the highest resolution at
+which pixel indices can be stored in a signed 64-bit integer. At this
+resolution, each pixel is scarcely 0.4 milliarcseconds across.
+
+The interval set representation is adventageous because there are simple and
+fast algorithms for interval arithmetic and set operations. Interval analysis
+appears in a suprising variety of scientific contexts from [genomics][pyranges]
+to [gravitational wave data quality][ligo-segments]. Because of the many
+business applications of interval arithmetic, intervals are also supported in
+the [PostgreSQL] database through its [range types].
+
+### Spatial primitives in `healpix_alchemy`
+
+The `healpix_alchemy` package provides two [custom column types] for
+[SQLAlchemy]:
+
+#### `healpix_alchemy.Point`
+
+This class represents a point. A column of this type could store the positions
+of galaxies in a catalog. Under the hood, it is just a `BIGINT`.
+
+Wherever you need to bind a Python value to a `healpix_alchemy.Point`, you may
+provide any one of the following:
+* an instance of [astropy.coordinates.SkyCoord]
+* a sequence of two [astropy.units.Quantity] instances with angle units, which
+  will be interpreted as the right ascension and declination of the point in
+  the [ICRS] frame
+* an integer representing the HEALPix NESTED index of the point at
+  _level_ = _level_<sub>max</sub>
+
+#### `healpix_alchemy.Tile`
+
+This class represents a HEALPix tile. A table containing a column of this type
+and a foreign key could store MOCs or gravitational-wave probability maps.
+Under the hood, it is just an `INT8RANGE`.
+
+Wherever you need to bind a Python value to a `healpix_alchemy.Tile`, you may
+provide any one of the following:
+* A single integer which will be interpreted as a the address of the tile in
+  the [UNIQ] HEALPix indexing scheme
+* A sequence of two integers, which will be interpreted as the lower and upper
+  bounds of the half-open pixel index interval at
+  _level_ = _level_<sub>max</sub>
+* A string like `'[1234,5678)'`
 
 ## Installation
 
-You can install `healpix_alchemy` the Python Package Index:
+You can install `healpix_alchemy` from the Python Package Index using [pip]:
 
-    $ pip install healpix-alchemy
+```console
+$ pip install healpix-alchemy
+```
 
-[SQLAlchemy]: https://www.sqlalchemy.org
+[pip]: https://pip.pypa.io/
+
+## Example
+
+First, some imports:
+
+```pycon
+>>> import sqlalchemy as sa
+>>> from sqlalchemy.ext.declarative import as_declarative, declared_attr
+>>> import healpix_alchemy as ha
+
+```
+
+### Set up tables
+
+This example will use the SQLAlchemy declarative extension for describing table
+schema using Python classes.
+
+SQLAlchemy needs to know the name for each table. You can provide the name by
+setting the `__tablename__` attribute in each model class, or you can
+create a base class that generates the table name automatically from the class
+name.
+
+```pycon
+>>> @as_declarative()
+... class Base:
+...
+...     @declared_attr
+...     def __tablename__(cls):
+...         return cls.__name__.lower()
+
+```
+
+Each row of the `Galaxy` table represents a point in a catalog:
+
+```pycon
+>>> class Galaxy(Base):
+...     id = sa.Column(sa.Text, primary_key=True)
+...     healpix = sa.Column(ha.Point, index=True, nullable=False)
+
+```
+
+Each row of the `Field` table represents a ZTF field:
+
+```pycon
+>>> class Field(Base):
+...     id = sa.Column(sa.Integer, primary_key=True)
+...     tiles = sa.orm.relationship(lambda: FieldTile)
+
+```
+
+Each row of the `FieldTile` represents a multi-resolution HEALPix tile that is
+contained within the corresponding field. There is a one-to-many mapping
+between `Field` and `FieldTile`.
+
+```pycon
+>>> class FieldTile(Base):
+...     id = sa.Column(sa.ForeignKey(Field.id), primary_key=True)
+...     healpix = sa.Column(ha.Tile, primary_key=True, index=True)
+
+```
+
+Each row of the `Localization` table represents a LIGO/Virgo HEALPix
+localization map.
+
+```pycon
+>>> class Localization(Base):
+...     id = sa.Column(sa.Integer, primary_key=True)
+...     tiles = sa.orm.relationship(lambda: LocalizationTile)
+
+```
+
+Each row of the `LocalizationTile` table represents a multi-resolution HEALPix
+tile within a LIGO/Virgo localization map. There is a one-to-many mapping
+between `Localization` and `LocalizationTile`.
+
+```pycon
+>>> class LocalizationTile(Base):
+...     id = sa.Column(sa.ForeignKey(Localization.id), primary_key=True)
+...     healpix = sa.Column(ha.Tile, primary_key=True, index=True)
+...     probdensity = sa.Column(sa.Float, nullable=False)
+
+```
+
+Finally, connect to the database, create all the tables, and start a session.
+
+```pycon
+>>> engine = sa.create_engine('postgresql://user:password@host/database')
+>>> Base.metadata.create_all(engine)
+>>> session = sa.orm.Session(engine)
+
+```
+
+### Populate with sample data
+
+Load the [2MASS Redshift Survey] into the `Galaxy` table. This catalog contains
+44599 galaxies.
+
+It may take up to a minute for this to finish. Advanced users may speed this up
+significantly by vectorizing the conversion from `SkyCoord` to HEALPix indices
+and using [SQLAlchemy bulk insertion].
+
+```pycon
+>>> from astropy.coordinates import SkyCoord
+>>> from astroquery.vizier import Vizier
+>>> vizier = Vizier(columns=['SimbadName', 'RAJ2000', 'DEJ2000'], row_limit=-1)
+>>> data, = vizier.get_catalogs('J/ApJS/199/26/table3')
+>>> data['coord'] = SkyCoord(data['RAJ2000'], data['DEJ2000'])
+>>> for row in data:
+...     session.add(Galaxy(id=row['SimbadName'], healpix=row['coord']))
+>>> session.commit()
+
+```
+
+Load the footprints of the [Zwicky Transient Facility] fields into the `Field`
+and `FieldTile` tables.
+
+It may take up to a minute for this to finish. Advanced users may speed this up
+significantly by using [SQLAlchemy bulk insertion].
+
+```pycon
+>>> from astropy.table import Table
+>>> from astropy.coordinates import SkyCoord
+>>> from astropy import units as u
+>>> url = 'https://raw.githubusercontent.com/ZwickyTransientFacility/ztf_information/master/field_grid/ztf_field_corners.csv'
+>>> for row in Table.read(url):
+...     field_id = int(row['field'])
+...     corners = SkyCoord([row['ra1'], row['ra2'], row['ra3'], row['ra4']],
+...                        [row['dec1'], row['dec2'], row['dec3'], row['dec4']],
+...                        unit=u.deg)
+...     tiles = [FieldTile(healpix=hpx) for hpx in ha.Tile.tiles_from(corners)]
+...     session.add(Field(id=field_id, tiles=tiles))
+>>> session.commit()
+
+```
+
+Load a sky map for LIGO/Virgo event [GW200115_042309] ([S200115j]) into the
+`Localization` and `LocalizationTile` tables.
+
+```pycon
+>>> url = 'https://gracedb.ligo.org/apiweb/superevents/S200115j/files/bayestar.multiorder.fits'
+>>> data = Table.read(url)
+>>> tiles = [LocalizationTile(healpix=row['UNIQ'], probdensity=row['PROBDENSITY']) for row in data]
+>>> session.add(Localization(id=1, tiles=tiles))
+>>> session.commit()
+
+```
+
+### Sample Queries
+
+#### What is the area of each field?
+
+```pycon
+>>> query = sa.select(
+...     FieldTile.id, sa.func.sum(FieldTile.healpix.area)
+... ).group_by(
+...     FieldTile.id
+... ).limit(
+...     5
+... )
+>>> for field_id, area in session.execute(query):
+...     print('Field', field_id, 'has area', area, 'steradians')
+Field 199 has area 0.017380122168993366 steradians
+Field 200 has area 0.017376127427358248 steradians
+Field 201 has area 0.017377126112767028 steradians
+Field 202 has area 0.017376127427358248 steradians
+Field 203 has area 0.017375128741949467 steradians
+
+```
+
+#### How many galaxies are in each field?
+
+```pycon
+>>> count = sa.func.count(Galaxy.id)
+>>> query = sa.select(
+...     FieldTile.id, count
+... ).filter(
+...     FieldTile.healpix.contains(Galaxy.healpix)
+... ).group_by(
+...     FieldTile.id
+... ).order_by(
+...     count.desc()
+... ).limit(
+...     5
+... )
+>>> for field_id, count in session.execute(query):
+...     print('Field', field_id, 'contains', count, 'galaxies')
+Field 1739 contains 343 galaxies
+Field 699 contains 336 galaxies
+Field 700 contains 311 galaxies
+Field 225 contains 303 galaxies
+Field 1740 contains 289 galaxies
+
+```
+
+#### What is the probability density at the position of each galaxy?
+
+```pycon
+>>> query = sa.select(
+...     Galaxy.id, LocalizationTile.probdensity
+... ).filter(
+...     LocalizationTile.id == 1,
+...     LocalizationTile.healpix.contains(Galaxy.healpix)
+... ).order_by(
+...     LocalizationTile.probdensity.desc()
+... ).limit(
+...     5
+... )
+>>> for galaxy_id, probdensity in session.execute(query):
+...     print('Galaxy', galaxy_id, 'has probability density', probdensity, 'per steradian')
+Galaxy 2MASX J02532153+0632222 has probability density 20.70119818753699 per steradian
+Galaxy 2MASX J02530482+0555431 has probability density 20.69509436913713 per steradian
+Galaxy 2MASX J02533119+0628252 has probability density 20.668990360570593 per steradian
+Galaxy 2MASX J02524584+0639206 has probability density 20.656073554123115 per steradian
+Galaxy 2MASX J02534120+0615562 has probability density 20.56675452367264 per steradian
+
+```
+
+#### What is the probability contained within each field?
+
+```pycon
+>>> prob = sa.func.sum(LocalizationTile.probdensity * (FieldTile.healpix * LocalizationTile.healpix).area)
+>>> query = sa.select(
+...     FieldTile.id, prob
+... ).filter(
+...     LocalizationTile.id == 1,
+...     FieldTile.healpix.overlaps(LocalizationTile.healpix)
+... ).group_by(
+...     FieldTile.id
+... ).order_by(
+...     prob.desc()
+... ).limit(
+...     5
+... )
+>>> for field_id, prob in engine.execute(query):
+...     print('Field', field_id, 'containment probability is', prob)
+Field 1499 containment probability is 0.1647961222933018
+Field 1446 containment probability is 0.15593990703399768
+Field 452 containment probability is 0.15420100565087086
+Field 505 containment probability is 0.09911749170919448
+Field 401 containment probability is 0.09621300206793466
+
+```
+
+#### What is the combined area of fields 1000 through 2000?
+
+In the next two examples, we introduce `healpix_alchemy.func.union()` which
+finds the union of a set of tiles. Because it is an aggregate function, it
+should generally be used in a subquery.
+
+```pycon
+>>> union = sa.select(
+...     ha.func.union(FieldTile.healpix).label('healpix')
+... ).filter(
+...     FieldTile.id.between(1000, 2000)
+... ).subquery()
+>>> query = sa.select(
+...     sa.func.sum(union.columns.healpix.area)
+... )
+>>> result = session.execute(query).scalar_one()
+>>> print(result, 'steradians')
+9.332762083260642 steradians
+
+```
+
+#### What is the integrated probability contained within fields 1000 through 2000?
+
+```pycon
+>>> union = sa.select(
+...     ha.func.union(FieldTile.healpix).label('healpix')
+... ).filter(
+...     FieldTile.id.between(1000, 2000)
+... ).subquery()
+>>> prob = sa.func.sum(LocalizationTile.probdensity * (union.columns.healpix * LocalizationTile.healpix).area)
+>>> query = sa.select(
+...     prob
+... ).filter(
+...     LocalizationTile.id == 1,
+...     union.columns.healpix.overlaps(LocalizationTile.healpix)
+... )
+>>> session.execute(query).scalar_one()
+0.8373173527131985
+
+```
+
+#### What is the area of the 90% credible region?
+
+```pycon
+>>> cum_area = sa.func.sum(
+...     LocalizationTile.healpix.area
+... ).over(
+...     order_by=LocalizationTile.probdensity.desc()
+... ).label(
+...     'cum_area'
+... )
+>>> cum_prob = sa.func.sum(
+...     LocalizationTile.probdensity * LocalizationTile.healpix.area
+... ).over(
+...     order_by=LocalizationTile.probdensity.desc()
+... ).label(
+...     'cum_prob'
+... )
+>>> subquery = sa.select(
+...     cum_area,
+...     cum_prob
+... ).filter(
+...     LocalizationTile.id == 1
+... ).subquery()
+>>> query = sa.select(
+...     sa.func.max(subquery.columns.cum_area)
+... ).filter(
+...     subquery.columns.cum_prob <= 0.9
+... )
+>>> result = session.execute(query).scalar_one()
+>>> print(result, 'steradians')
+0.2766278687487106 steradians
+
+```
+
+#### Which galaxies are within the 90% credible region?
+
+```pycon
+>>> cum_prob = sa.func.sum(
+...     LocalizationTile.probdensity * LocalizationTile.healpix.area
+... ).over(
+...     order_by=LocalizationTile.probdensity.desc()
+... ).label(
+...     'cum_prob'
+... )
+>>> subquery1 = sa.select(
+...     LocalizationTile.probdensity,
+...     cum_prob
+... ).filter(
+...     LocalizationTile.id == 1
+... ).subquery()
+>>> min_probdensity = sa.select(
+...     sa.func.min(subquery1.columns.probdensity)
+... ).filter(
+...     subquery1.columns.cum_prob <= 0.9
+... ).scalar_subquery()
+>>> query = sa.select(
+...     Galaxy.id
+... ).filter(
+...     LocalizationTile.id == 1,
+...     LocalizationTile.healpix.contains(Galaxy.healpix),
+...     LocalizationTile.probdensity >= min_probdensity
+... ).limit(
+...     5
+... )
+>>> for galaxy_id, in session.execute(query):
+...     print(galaxy_id)
+2MASX J02424077-0000478
+2MASX J02352772-0921216
+2MASX J02273746-0109226
+2MASX J02414523+0026354
+2MASX J20095408-4822462
+
+```
diff --git a/healpix_alchemy/__init__.py b/healpix_alchemy/__init__.py
index d414eb9..5a627e0 100644
--- a/healpix_alchemy/__init__.py
+++ b/healpix_alchemy/__init__.py
@@ -1,3 +1,4 @@
-from .healpix import Tile, Region
+from .types import Point, Tile
+from . import func  # noqa: F401
 
-__all__ = ('Tile', 'Region')
+__all__ = ('Point', 'Tile')
diff --git a/healpix_alchemy/constants.py b/healpix_alchemy/constants.py
new file mode 100644
index 0000000..1e54eec
--- /dev/null
+++ b/healpix_alchemy/constants.py
@@ -0,0 +1,18 @@
+from astropy.coordinates import ICRS
+from astropy import units as u
+from astropy_healpix import level_to_nside, HEALPix
+from mocpy import IntervalSet
+import sqlalchemy as sa
+
+LEVEL = IntervalSet.HPX_MAX_ORDER
+"""Base HEALPix resolution. This is the maximum HEALPix level that can be
+stored in a signed 8-byte integer data type."""
+
+HPX = HEALPix(nside=level_to_nside(LEVEL), order='nested', frame=ICRS())
+"""HEALPix projection object."""
+
+PIXEL_AREA = HPX.pixel_area.to_value(u.sr)
+"""Native pixel area in steradians."""
+
+PIXEL_AREA_LITERAL = sa.literal(PIXEL_AREA, sa.Float)
+"""Pixel area as an SQL literal."""
diff --git a/healpix_alchemy/func.py b/healpix_alchemy/func.py
new file mode 100644
index 0000000..35b2e9f
--- /dev/null
+++ b/healpix_alchemy/func.py
@@ -0,0 +1,8 @@
+"""SQLAlchemy functions."""
+from sqlalchemy import func as _func
+
+from .types import Tile as _Tile
+
+
+def union(tiles):
+    return _func.unnest(_func.range_agg(tiles), type_=_Tile)
diff --git a/healpix_alchemy/healpix.py b/healpix_alchemy/healpix.py
deleted file mode 100644
index d01c43c..0000000
--- a/healpix_alchemy/healpix.py
+++ /dev/null
@@ -1,117 +0,0 @@
-"""Model base classes for multiresolution HEALPix data."""
-from astropy.coordinates import ICRS
-from astropy import units as u
-from astropy_healpix import level_to_nside, uniq_to_level_ipix, HEALPix
-from mocpy import MOC, IntervalSet
-import numpy as np
-from sqlalchemy import BigInteger, Column
-from sqlalchemy.ext.declarative import declared_attr
-from sqlalchemy.orm import relationship
-from sqlalchemy.orm.mapper import Mapper
-
-LEVEL = IntervalSet.HPX_MAX_ORDER
-"""Base HEALPix resolution. This is the maximum HEALPix level that can be
-stored in a signed 8-byte integer data type."""
-
-HPX = HEALPix(nside=level_to_nside(LEVEL), order='nested', frame=ICRS())
-"""HEALPix projection object."""
-
-PIXEL_AREA = HPX.pixel_area.to_value(u.sr)
-"""Native pixel area in steradians."""
-
-
-class Tile:
-    """Mixin class for a table that stores a HEALPix multiresolution tile."""
-
-    nested_lo = Column(
-        BigInteger, index=True, nullable=False, primary_key=True,
-        doc=('Lower end of range of HEALPix nested indices at',
-             f'nside=2**{LEVEL}'))
-
-    nested_hi = Column(
-        BigInteger, nullable=False,
-        doc=('Upper end of range of HEALPix nested indices at',
-             f'nside=2**{LEVEL}'))
-
-    def __init__(self, *args, uniq=None, **kwargs):
-        super().__init__(*args, **kwargs)
-        if uniq is not None:
-            self.uniq = uniq
-
-    @property
-    def uniq(self):
-        """HEALPix UNIQ pixel index."""
-        # This is the same expression as in astropy_healpix.level_ipix_to_uniq,
-        # but reproduced here so that SQLAlchemy can map it to SQL.
-        level = int(LEVEL - (np.ceil((np.log2(self.nested_hi -
-                                              self.nested_lo)/2))))
-        shift = 2 * (LEVEL - level)
-        ipix = int(self.nested_lo >> shift)
-        return int(ipix + (1 << 2 * (level + 1)))
-
-    @uniq.setter
-    def uniq(self, value):
-        """HEALPix UNIQ pixel index."""
-        level, ipix = uniq_to_level_ipix(value)
-        shift = 2 * (LEVEL - level)
-        self.nested_lo = int(ipix << shift)
-        self.nested_hi = int(((ipix + 1) << shift) - 1)
-
-
-def _class_or_lambda(cls):
-    """Shim to allow subclasses to set tile_class to a lambda function."""
-    # Adopted from sqlalchemy.orm.relationships.RelationshipProperty.entity
-    if callable(cls) and not isinstance(cls, (type, Mapper)):
-        cls = cls()
-    return cls
-
-
-class Region:
-    """Mixin class for a HEALPix multiresolution region (e.g., multi-order
-    coverage map)."""
-
-    @declared_attr
-    def tiles(cls):
-        return relationship(cls.tile_class, backref=cls.__tablename__)
-
-    @classmethod
-    def from_moc(cls, moc, tile_args={}, *args, **kwargs):
-        """Create a new instance from a multi-order coverage map.
-
-        Parameters
-        ----------
-        moc : mocpy.MOC
-            The multi-order coverage map.
-
-        Returns
-        -------
-        list
-            A list of Tile instances.
-        """
-        tile_class = _class_or_lambda(cls.tile_class)
-        nested_ranges = moc._interval_set.nested
-        # FIXME: MOCpy should return an array of size(0, 2) 0 for an empty MOC,
-        # but it actually returns an array of size(1, 0).
-        if nested_ranges.size == 0:
-            nested_ranges = nested_ranges.reshape(-1, 2)
-        tiles = [tile_class(nested_lo=int(lo), nested_hi=int(hi),
-                            **tile_args)
-                 for lo, hi in nested_ranges]
-        return cls(*args, tiles=tiles, **kwargs)
-
-    @classmethod
-    def from_polygon(cls, polygon, *args, **kwargs):
-        """Create a new instance from a polygon.
-
-        Parameters
-        ----------
-        polygon : astropy.coordinates.SkyCoord
-            The vertices of the polygon.
-
-        Returns
-        -------
-        list
-            A list of Tile instances.
-        """
-        moc = MOC.from_polygon_skycoord(polygon)
-        return cls.from_moc(moc, *args, **kwargs)
diff --git a/healpix_alchemy/types.py b/healpix_alchemy/types.py
new file mode 100644
index 0000000..4f6455c
--- /dev/null
+++ b/healpix_alchemy/types.py
@@ -0,0 +1,96 @@
+"""SQLAlchemy types for multiresolution HEALPix data."""
+from collections.abc import Sequence
+from numbers import Integral
+
+from astropy.coordinates import SkyCoord
+from astropy_healpix import uniq_to_level_ipix
+from mocpy import MOC
+import numpy as np
+import sqlalchemy as sa
+from sqlalchemy.dialects.postgresql import INT8RANGE
+
+from .constants import HPX, LEVEL, PIXEL_AREA_LITERAL
+
+__all__ = ('Point', 'Tile')
+
+
+class Point(sa.TypeDecorator):
+
+    cache_ok = True
+    impl = sa.BigInteger
+
+    def process_bind_param(self, value, dialect):
+        if isinstance(value, SkyCoord):
+            value = HPX.skycoord_to_healpix(value)
+        elif isinstance(value, Sequence) and len(value) == 2:
+            value = HPX.lonlat_to_healpix(*value)
+        if isinstance(value, np.int64):
+            value = int(value)
+        return value
+
+
+class Tile(sa.TypeDecorator):
+
+    cache_ok = True
+    impl = INT8RANGE
+
+    def process_bind_param(self, value, dialect):
+        if isinstance(value, Integral):
+            level, ipix = uniq_to_level_ipix(value)
+            shift = 2 * (LEVEL - level)
+            value = (ipix << shift, (ipix + 1) << shift)
+        if isinstance(value, Sequence) and len(value) == 2:
+            value = f'[{value[0]},{value[1]})'
+        return value
+
+    class comparator_factory(INT8RANGE.comparator_factory):
+
+        @property
+        def lower(self):
+            return sa.func.lower(self, type_=Point)
+
+        @property
+        def upper(self):
+            return sa.func.upper(self, type_=Point)
+
+        @property
+        def length(self):
+            return self.upper - self.lower
+
+        @property
+        def area(self):
+            return self.length * PIXEL_AREA_LITERAL
+
+    @classmethod
+    def tiles_from(cls, obj):
+        if isinstance(obj, MOC):
+            return cls.tiles_from_moc(obj)
+        elif isinstance(obj, SkyCoord):
+            return cls.tiles_from_polygon_skycoord(obj)
+        else:
+            raise TypeError('Unknown type')
+
+    @classmethod
+    def tiles_from_polygon_skycoord(cls, polygon):
+        return cls.tiles_from_moc(
+            MOC.from_polygon_skycoord(
+                polygon.transform_to(HPX.frame)))
+
+    @classmethod
+    def tiles_from_moc(cls, moc):
+        return [tuple(interval) for interval in moc._interval_set.nested]
+
+
+@sa.event.listens_for(sa.Index, 'after_parent_attach')
+def _create_indices(index, parent):
+    """Set index method to SP-GiST_ for any indexed Tile or Region columns.
+
+    .. _SP-GiST: https://www.postgresql.org/docs/current/spgist.html
+    """
+    if (
+        index._column_flag and
+        len(index.expressions) == 1 and
+        isinstance(index.expressions[0], sa.Column) and
+        isinstance(index.expressions[0].type, Tile)
+    ):
+        index.dialect_options['postgresql']['using'] = 'spgist'
diff --git a/poetry.lock b/poetry.lock
index cd9dc2a..9af6ca2 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -32,6 +32,30 @@ numpy = "*"
 docs = ["sphinx-astropy", "matplotlib"]
 test = ["pytest-astropy", "hypothesis (<4.42)"]
 
+[[package]]
+name = "astroquery"
+version = "0.4.3"
+description = "Functions and classes to access online astronomical data resources"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[package.dependencies]
+astropy = ">=3.1.2"
+beautifulsoup4 = ">=4.3.2"
+html5lib = ">=0.999"
+keyring = ">=4.0"
+numpy = ">=1.15.0"
+pyvo = ">=1.1"
+requests = ">=2.4.3"
+six = "*"
+
+[package.extras]
+all = ["mocpy (>=0.5.2)", "astropy-healpix", "boto3", "regions"]
+all_lt_39 = ["mocpy (>=0.5.2)", "regions", "pyregion", "astropy-healpix", "aplpy", "boto3"]
+docs = ["sphinx-astropy"]
+test = ["pytest-astropy", "matplotlib", "jinja2", "flask", "pytest-dependency"]
+
 [[package]]
 name = "atomicwrites"
 version = "1.4.0"
@@ -54,6 +78,21 @@ docs = ["furo", "sphinx", "zope.interface", "sphinx-notfound-page"]
 tests = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface"]
 tests_no_zope = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins"]
 
+[[package]]
+name = "beautifulsoup4"
+version = "4.10.0"
+description = "Screen-scraping library"
+category = "dev"
+optional = false
+python-versions = ">3.0.0"
+
+[package.dependencies]
+soupsieve = ">1.2"
+
+[package.extras]
+html5lib = ["html5lib"]
+lxml = ["lxml"]
+
 [[package]]
 name = "cdshealpix"
 version = "0.6.1"
@@ -65,6 +104,36 @@ python-versions = "*"
 [package.dependencies]
 astropy = "*"
 
+[[package]]
+name = "certifi"
+version = "2021.10.8"
+description = "Python package for providing Mozilla's CA Bundle."
+category = "dev"
+optional = false
+python-versions = "*"
+
+[[package]]
+name = "cffi"
+version = "1.15.0"
+description = "Foreign Function Interface for Python calling C code."
+category = "dev"
+optional = false
+python-versions = "*"
+
+[package.dependencies]
+pycparser = "*"
+
+[[package]]
+name = "charset-normalizer"
+version = "2.0.7"
+description = "The Real First Universal Charset Detector. Open, modern and actively maintained alternative to Chardet."
+category = "dev"
+optional = false
+python-versions = ">=3.5.0"
+
+[package.extras]
+unicode_backport = ["unicodedata2"]
+
 [[package]]
 name = "colorama"
 version = "0.4.4"
@@ -73,6 +142,25 @@ category = "dev"
 optional = false
 python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
 
+[[package]]
+name = "cryptography"
+version = "35.0.0"
+description = "cryptography is a package which provides cryptographic recipes and primitives to Python developers."
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.dependencies]
+cffi = ">=1.12"
+
+[package.extras]
+docs = ["sphinx (>=1.6.5,!=1.8.0,!=3.1.0,!=3.1.1)", "sphinx-rtd-theme"]
+docstest = ["doc8", "pyenchant (>=1.6.11)", "twine (>=1.12.0)", "sphinxcontrib-spelling (>=4.0.1)"]
+pep8test = ["black", "flake8", "flake8-import-order", "pep8-naming"]
+sdist = ["setuptools_rust (>=0.11.4)"]
+ssh = ["bcrypt (>=3.1.5)"]
+test = ["pytest (>=6.2.0)", "pytest-cov", "pytest-subtests", "pytest-xdist", "pretend", "iso8601", "pytz", "hypothesis (>=1.11.4,!=3.79.2)"]
+
 [[package]]
 name = "cycler"
 version = "0.11.0"
@@ -113,6 +201,32 @@ python-versions = ">=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*"
 [package.extras]
 docs = ["sphinx"]
 
+[[package]]
+name = "html5lib"
+version = "1.1"
+description = "HTML parser based on the WHATWG HTML specification"
+category = "dev"
+optional = false
+python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
+
+[package.dependencies]
+six = ">=1.9"
+webencodings = "*"
+
+[package.extras]
+all = ["genshi", "chardet (>=2.2)", "lxml"]
+chardet = ["chardet (>=2.2)"]
+genshi = ["genshi"]
+lxml = ["lxml"]
+
+[[package]]
+name = "idna"
+version = "3.3"
+description = "Internationalized Domain Names in Applications (IDNA)"
+category = "dev"
+optional = false
+python-versions = ">=3.5"
+
 [[package]]
 name = "importlib-metadata"
 version = "4.8.2"
@@ -138,6 +252,36 @@ category = "dev"
 optional = false
 python-versions = "*"
 
+[[package]]
+name = "jeepney"
+version = "0.7.1"
+description = "Low-level, pure Python DBus protocol wrapper."
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.extras]
+test = ["pytest", "pytest-trio", "pytest-asyncio", "testpath", "trio", "async-timeout"]
+trio = ["trio", "async-generator"]
+
+[[package]]
+name = "keyring"
+version = "23.2.1"
+description = "Store and access your passwords safely."
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.dependencies]
+importlib-metadata = ">=3.6"
+jeepney = {version = ">=0.4.2", markers = "sys_platform == \"linux\""}
+pywin32-ctypes = {version = "<0.1.0 || >0.1.0,<0.1.1 || >0.1.1", markers = "sys_platform == \"win32\""}
+SecretStorage = {version = ">=3.2", markers = "sys_platform == \"linux\""}
+
+[package.extras]
+docs = ["sphinx", "jaraco.packaging (>=8.2)", "rst.linker (>=1.9)"]
+testing = ["pytest (>=4.6)", "pytest-checkdocs (>=2.4)", "pytest-flake8", "pytest-cov", "pytest-enabler (>=1.0.1)", "pytest-black (>=0.3.7)", "pytest-mypy"]
+
 [[package]]
 name = "kiwisolver"
 version = "1.3.2"
@@ -178,6 +322,14 @@ pyparsing = ">=2.2.1"
 python-dateutil = ">=2.7"
 setuptools_scm = ">=4"
 
+[[package]]
+name = "mimeparse"
+version = "0.1.3"
+description = "A module provides basic functions for parsing mime-type names and matching them against a list of media-ranges."
+category = "dev"
+optional = false
+python-versions = "*"
+
 [[package]]
 name = "mirakuru"
 version = "2.4.1"
@@ -326,6 +478,14 @@ category = "dev"
 optional = false
 python-versions = "*"
 
+[[package]]
+name = "pycparser"
+version = "2.21"
+description = "C parser in Python"
+category = "dev"
+optional = false
+python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
+
 [[package]]
 name = "pyerfa"
 version = "2.0.0.1"
@@ -388,6 +548,21 @@ aspect = ["aspectlib"]
 elasticsearch = ["elasticsearch"]
 histogram = ["pygal", "pygaljs"]
 
+[[package]]
+name = "pytest-doctestplus"
+version = "0.11.0"
+description = "Pytest plugin with advanced doctest features."
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[package.dependencies]
+packaging = ">=17.0"
+pytest = ">=4.6"
+
+[package.extras]
+test = ["pytest-remotedata (>=0.3.2)", "sphinx"]
+
 [[package]]
 name = "pytest-postgresql"
 version = "3.1.2"
@@ -423,6 +598,57 @@ category = "dev"
 optional = false
 python-versions = "*"
 
+[[package]]
+name = "pyvo"
+version = "1.1"
+description = "Astropy affiliated package for accessing Virtual Observatory data and services"
+category = "dev"
+optional = false
+python-versions = ">=3.5"
+
+[package.dependencies]
+astropy = "*"
+mimeparse = "*"
+requests = "*"
+
+[[package]]
+name = "pywin32-ctypes"
+version = "0.2.0"
+description = ""
+category = "dev"
+optional = false
+python-versions = "*"
+
+[[package]]
+name = "requests"
+version = "2.26.0"
+description = "Python HTTP for Humans."
+category = "dev"
+optional = false
+python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*, !=3.5.*"
+
+[package.dependencies]
+certifi = ">=2017.4.17"
+charset-normalizer = {version = ">=2.0.0,<2.1.0", markers = "python_version >= \"3\""}
+idna = {version = ">=2.5,<4", markers = "python_version >= \"3\""}
+urllib3 = ">=1.21.1,<1.27"
+
+[package.extras]
+socks = ["PySocks (>=1.5.6,!=1.5.7)", "win-inet-pton"]
+use_chardet_on_py3 = ["chardet (>=3.0.2,<5)"]
+
+[[package]]
+name = "secretstorage"
+version = "3.3.1"
+description = "Python bindings to FreeDesktop.org Secret Service API"
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.dependencies]
+cryptography = ">=2.0"
+jeepney = ">=0.6"
+
 [[package]]
 name = "setuptools-scm"
 version = "6.3.2"
@@ -446,6 +672,14 @@ category = "main"
 optional = false
 python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*"
 
+[[package]]
+name = "soupsieve"
+version = "2.3.1"
+description = "A modern CSS selector implementation for Beautiful Soup."
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
 [[package]]
 name = "sqlalchemy"
 version = "1.4.27"
@@ -503,6 +737,27 @@ category = "main"
 optional = false
 python-versions = ">=3.6"
 
+[[package]]
+name = "urllib3"
+version = "1.26.7"
+description = "HTTP library with thread-safe connection pooling, file post, and more."
+category = "dev"
+optional = false
+python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*, <4"
+
+[package.extras]
+brotli = ["brotlipy (>=0.6.0)"]
+secure = ["pyOpenSSL (>=0.14)", "cryptography (>=1.3.4)", "idna (>=2.0.0)", "certifi", "ipaddress"]
+socks = ["PySocks (>=1.5.6,!=1.5.7,<2.0)"]
+
+[[package]]
+name = "webencodings"
+version = "0.5.1"
+description = "Character encoding aliases for legacy web content"
+category = "dev"
+optional = false
+python-versions = "*"
+
 [[package]]
 name = "zipp"
 version = "3.6.0"
@@ -518,7 +773,7 @@ testing = ["pytest (>=4.6)", "pytest-checkdocs (>=2.4)", "pytest-flake8", "pytes
 [metadata]
 lock-version = "1.1"
 python-versions = "^3.7"
-content-hash = "65d4ec8a708f7cd25469de3f3c79c84682201dee44501a9e18942bd7dd52ad5e"
+content-hash = "3fc5d006ccfbda4867afdbd8a09965be581f91a7a3b2ed4373be6e83d3b70237"
 
 [metadata.files]
 astropy = [
@@ -570,6 +825,10 @@ astropy-healpix = [
     {file = "astropy_healpix-0.6-cp39-cp39-win_amd64.whl", hash = "sha256:46f7d2918079281d5be5899a286828904a01b65f088ab564f7c1e8636914cea1"},
     {file = "astropy_healpix-0.6.tar.gz", hash = "sha256:409a6621c383641456c074f0f0350a24a4a58e910eaeef14e9bbce3e00ad6690"},
 ]
+astroquery = [
+    {file = "astroquery-0.4.3-py3-none-any.whl", hash = "sha256:31fb96f6fa122384f12c3088cb9b62e2e2489d062f96de3cf3cf5939ffd9b99a"},
+    {file = "astroquery-0.4.3.tar.gz", hash = "sha256:5ee4b47089a5ba106963be22ac3a7f7387f30d932091e340243841fe77463f90"},
+]
 atomicwrites = [
     {file = "atomicwrites-1.4.0-py2.py3-none-any.whl", hash = "sha256:6d1784dea7c0c8d4a5172b6c620f40b6e4cbfdf96d783691f2e1302a7b88e197"},
     {file = "atomicwrites-1.4.0.tar.gz", hash = "sha256:ae70396ad1a434f9c7046fd2dd196fc04b12f9e91ffb859164193be8b6168a7a"},
@@ -578,6 +837,10 @@ attrs = [
     {file = "attrs-21.2.0-py2.py3-none-any.whl", hash = "sha256:149e90d6d8ac20db7a955ad60cf0e6881a3f20d37096140088356da6c716b0b1"},
     {file = "attrs-21.2.0.tar.gz", hash = "sha256:ef6aaac3ca6cd92904cdd0d83f629a15f18053ec84e6432106f7a4d04ae4f5fb"},
 ]
+beautifulsoup4 = [
+    {file = "beautifulsoup4-4.10.0-py3-none-any.whl", hash = "sha256:9a315ce70049920ea4572a4055bc4bd700c940521d36fc858205ad4fcde149bf"},
+    {file = "beautifulsoup4-4.10.0.tar.gz", hash = "sha256:c23ad23c521d818955a4151a67d81580319d4bf548d3d49f4223ae041ff98891"},
+]
 cdshealpix = [
     {file = "cdshealpix-0.6.1-cp310-cp310-macosx_10_7_x86_64.whl", hash = "sha256:547875af55feb6e9e60b2a553eaee8ecd9be652264310b8b36c27ffb86b21fbf"},
     {file = "cdshealpix-0.6.1-cp310-cp310-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:523c607bdcb8274bf2923f9730eaeb78cbba4955fa4d15c37abce5ab43ce876c"},
@@ -601,10 +864,92 @@ cdshealpix = [
     {file = "cdshealpix-0.6.1-cp39-none-win_amd64.whl", hash = "sha256:bdf2a2074dd9f9351eb5cca72f8196c4d234c77c2d4355f65c3950a39a441caa"},
     {file = "cdshealpix-0.6.1.tar.gz", hash = "sha256:72b33dea542ceedb3fb8d2e13921a1b8b5249439758edad15dab7803f38abf67"},
 ]
+certifi = [
+    {file = "certifi-2021.10.8-py2.py3-none-any.whl", hash = "sha256:d62a0163eb4c2344ac042ab2bdf75399a71a2d8c7d47eac2e2ee91b9d6339569"},
+    {file = "certifi-2021.10.8.tar.gz", hash = "sha256:78884e7c1d4b00ce3cea67b44566851c4343c120abd683433ce934a68ea58872"},
+]
+cffi = [
+    {file = "cffi-1.15.0-cp27-cp27m-macosx_10_9_x86_64.whl", hash = "sha256:c2502a1a03b6312837279c8c1bd3ebedf6c12c4228ddbad40912d671ccc8a962"},
+    {file = "cffi-1.15.0-cp27-cp27m-manylinux1_i686.whl", hash = "sha256:23cfe892bd5dd8941608f93348c0737e369e51c100d03718f108bf1add7bd6d0"},
+    {file = "cffi-1.15.0-cp27-cp27m-manylinux1_x86_64.whl", hash = "sha256:41d45de54cd277a7878919867c0f08b0cf817605e4eb94093e7516505d3c8d14"},
+    {file = "cffi-1.15.0-cp27-cp27m-win32.whl", hash = "sha256:4a306fa632e8f0928956a41fa8e1d6243c71e7eb59ffbd165fc0b41e316b2474"},
+    {file = "cffi-1.15.0-cp27-cp27m-win_amd64.whl", hash = "sha256:e7022a66d9b55e93e1a845d8c9eba2a1bebd4966cd8bfc25d9cd07d515b33fa6"},
+    {file = "cffi-1.15.0-cp27-cp27mu-manylinux1_i686.whl", hash = "sha256:14cd121ea63ecdae71efa69c15c5543a4b5fbcd0bbe2aad864baca0063cecf27"},
+    {file = "cffi-1.15.0-cp27-cp27mu-manylinux1_x86_64.whl", hash = "sha256:d4d692a89c5cf08a8557fdeb329b82e7bf609aadfaed6c0d79f5a449a3c7c023"},
+    {file = "cffi-1.15.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:0104fb5ae2391d46a4cb082abdd5c69ea4eab79d8d44eaaf79f1b1fd806ee4c2"},
+    {file = "cffi-1.15.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:91ec59c33514b7c7559a6acda53bbfe1b283949c34fe7440bcf917f96ac0723e"},
+    {file = "cffi-1.15.0-cp310-cp310-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:f5c7150ad32ba43a07c4479f40241756145a1f03b43480e058cfd862bf5041c7"},
+    {file = "cffi-1.15.0-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:00c878c90cb53ccfaae6b8bc18ad05d2036553e6d9d1d9dbcf323bbe83854ca3"},
+    {file = "cffi-1.15.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:abb9a20a72ac4e0fdb50dae135ba5e77880518e742077ced47eb1499e29a443c"},
+    {file = "cffi-1.15.0-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:a5263e363c27b653a90078143adb3d076c1a748ec9ecc78ea2fb916f9b861962"},
+    {file = "cffi-1.15.0-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:f54a64f8b0c8ff0b64d18aa76675262e1700f3995182267998c31ae974fbc382"},
+    {file = "cffi-1.15.0-cp310-cp310-win32.whl", hash = "sha256:c21c9e3896c23007803a875460fb786118f0cdd4434359577ea25eb556e34c55"},
+    {file = "cffi-1.15.0-cp310-cp310-win_amd64.whl", hash = "sha256:5e069f72d497312b24fcc02073d70cb989045d1c91cbd53979366077959933e0"},
+    {file = "cffi-1.15.0-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:64d4ec9f448dfe041705426000cc13e34e6e5bb13736e9fd62e34a0b0c41566e"},
+    {file = "cffi-1.15.0-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:2756c88cbb94231c7a147402476be2c4df2f6078099a6f4a480d239a8817ae39"},
+    {file = "cffi-1.15.0-cp36-cp36m-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:3b96a311ac60a3f6be21d2572e46ce67f09abcf4d09344c49274eb9e0bf345fc"},
+    {file = "cffi-1.15.0-cp36-cp36m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:75e4024375654472cc27e91cbe9eaa08567f7fbdf822638be2814ce059f58032"},
+    {file = "cffi-1.15.0-cp36-cp36m-manylinux_2_5_i686.manylinux1_i686.whl", hash = "sha256:59888172256cac5629e60e72e86598027aca6bf01fa2465bdb676d37636573e8"},
+    {file = "cffi-1.15.0-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:27c219baf94952ae9d50ec19651a687b826792055353d07648a5695413e0c605"},
+    {file = "cffi-1.15.0-cp36-cp36m-win32.whl", hash = "sha256:4958391dbd6249d7ad855b9ca88fae690783a6be9e86df65865058ed81fc860e"},
+    {file = "cffi-1.15.0-cp36-cp36m-win_amd64.whl", hash = "sha256:f6f824dc3bce0edab5f427efcfb1d63ee75b6fcb7282900ccaf925be84efb0fc"},
+    {file = "cffi-1.15.0-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:06c48159c1abed75c2e721b1715c379fa3200c7784271b3c46df01383b593636"},
+    {file = "cffi-1.15.0-cp37-cp37m-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:c2051981a968d7de9dd2d7b87bcb9c939c74a34626a6e2f8181455dd49ed69e4"},
+    {file = "cffi-1.15.0-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:fd8a250edc26254fe5b33be00402e6d287f562b6a5b2152dec302fa15bb3e997"},
+    {file = "cffi-1.15.0-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:91d77d2a782be4274da750752bb1650a97bfd8f291022b379bb8e01c66b4e96b"},
+    {file = "cffi-1.15.0-cp37-cp37m-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:45db3a33139e9c8f7c09234b5784a5e33d31fd6907800b316decad50af323ff2"},
+    {file = "cffi-1.15.0-cp37-cp37m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:263cc3d821c4ab2213cbe8cd8b355a7f72a8324577dc865ef98487c1aeee2bc7"},
+    {file = "cffi-1.15.0-cp37-cp37m-win32.whl", hash = "sha256:17771976e82e9f94976180f76468546834d22a7cc404b17c22df2a2c81db0c66"},
+    {file = "cffi-1.15.0-cp37-cp37m-win_amd64.whl", hash = "sha256:3415c89f9204ee60cd09b235810be700e993e343a408693e80ce7f6a40108029"},
+    {file = "cffi-1.15.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:4238e6dab5d6a8ba812de994bbb0a79bddbdf80994e4ce802b6f6f3142fcc880"},
+    {file = "cffi-1.15.0-cp38-cp38-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:0808014eb713677ec1292301ea4c81ad277b6cdf2fdd90fd540af98c0b101d20"},
+    {file = "cffi-1.15.0-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:57e9ac9ccc3101fac9d6014fba037473e4358ef4e89f8e181f8951a2c0162024"},
+    {file = "cffi-1.15.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:8b6c2ea03845c9f501ed1313e78de148cd3f6cad741a75d43a29b43da27f2e1e"},
+    {file = "cffi-1.15.0-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:10dffb601ccfb65262a27233ac273d552ddc4d8ae1bf93b21c94b8511bffe728"},
+    {file = "cffi-1.15.0-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:786902fb9ba7433aae840e0ed609f45c7bcd4e225ebb9c753aa39725bb3e6ad6"},
+    {file = "cffi-1.15.0-cp38-cp38-win32.whl", hash = "sha256:da5db4e883f1ce37f55c667e5c0de439df76ac4cb55964655906306918e7363c"},
+    {file = "cffi-1.15.0-cp38-cp38-win_amd64.whl", hash = "sha256:181dee03b1170ff1969489acf1c26533710231c58f95534e3edac87fff06c443"},
+    {file = "cffi-1.15.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:45e8636704eacc432a206ac7345a5d3d2c62d95a507ec70d62f23cd91770482a"},
+    {file = "cffi-1.15.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:31fb708d9d7c3f49a60f04cf5b119aeefe5644daba1cd2a0fe389b674fd1de37"},
+    {file = "cffi-1.15.0-cp39-cp39-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:6dc2737a3674b3e344847c8686cf29e500584ccad76204efea14f451d4cc669a"},
+    {file = "cffi-1.15.0-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:74fdfdbfdc48d3f47148976f49fab3251e550a8720bebc99bf1483f5bfb5db3e"},
+    {file = "cffi-1.15.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ffaa5c925128e29efbde7301d8ecaf35c8c60ffbcd6a1ffd3a552177c8e5e796"},
+    {file = "cffi-1.15.0-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:3f7d084648d77af029acb79a0ff49a0ad7e9d09057a9bf46596dac9514dc07df"},
+    {file = "cffi-1.15.0-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:ef1f279350da2c586a69d32fc8733092fd32cc8ac95139a00377841f59a3f8d8"},
+    {file = "cffi-1.15.0-cp39-cp39-win32.whl", hash = "sha256:2a23af14f408d53d5e6cd4e3d9a24ff9e05906ad574822a10563efcef137979a"},
+    {file = "cffi-1.15.0-cp39-cp39-win_amd64.whl", hash = "sha256:3773c4d81e6e818df2efbc7dd77325ca0dcb688116050fb2b3011218eda36139"},
+    {file = "cffi-1.15.0.tar.gz", hash = "sha256:920f0d66a896c2d99f0adbb391f990a84091179542c205fa53ce5787aff87954"},
+]
+charset-normalizer = [
+    {file = "charset-normalizer-2.0.7.tar.gz", hash = "sha256:e019de665e2bcf9c2b64e2e5aa025fa991da8720daa3c1138cadd2fd1856aed0"},
+    {file = "charset_normalizer-2.0.7-py3-none-any.whl", hash = "sha256:f7af805c321bfa1ce6714c51f254e0d5bb5e5834039bc17db7ebe3a4cec9492b"},
+]
 colorama = [
     {file = "colorama-0.4.4-py2.py3-none-any.whl", hash = "sha256:9f47eda37229f68eee03b24b9748937c7dc3868f906e8ba69fbcbdd3bc5dc3e2"},
     {file = "colorama-0.4.4.tar.gz", hash = "sha256:5941b2b48a20143d2267e95b1c2a7603ce057ee39fd88e7329b0c292aa16869b"},
 ]
+cryptography = [
+    {file = "cryptography-35.0.0-cp36-abi3-macosx_10_10_x86_64.whl", hash = "sha256:d57e0cdc1b44b6cdf8af1d01807db06886f10177469312fbde8f44ccbb284bc9"},
+    {file = "cryptography-35.0.0-cp36-abi3-macosx_11_0_arm64.whl", hash = "sha256:ced40344e811d6abba00295ced98c01aecf0c2de39481792d87af4fa58b7b4d6"},
+    {file = "cryptography-35.0.0-cp36-abi3-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:54b2605e5475944e2213258e0ab8696f4f357a31371e538ef21e8d61c843c28d"},
+    {file = "cryptography-35.0.0-cp36-abi3-manylinux_2_17_aarch64.manylinux2014_aarch64.manylinux_2_24_aarch64.whl", hash = "sha256:7b7ceeff114c31f285528ba8b390d3e9cfa2da17b56f11d366769a807f17cbaa"},
+    {file = "cryptography-35.0.0-cp36-abi3-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:2d69645f535f4b2c722cfb07a8eab916265545b3475fdb34e0be2f4ee8b0b15e"},
+    {file = "cryptography-35.0.0-cp36-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:4a2d0e0acc20ede0f06ef7aa58546eee96d2592c00f450c9acb89c5879b61992"},
+    {file = "cryptography-35.0.0-cp36-abi3-manylinux_2_24_x86_64.whl", hash = "sha256:07bb7fbfb5de0980590ddfc7f13081520def06dc9ed214000ad4372fb4e3c7f6"},
+    {file = "cryptography-35.0.0-cp36-abi3-musllinux_1_1_aarch64.whl", hash = "sha256:7eba2cebca600a7806b893cb1d541a6e910afa87e97acf2021a22b32da1df52d"},
+    {file = "cryptography-35.0.0-cp36-abi3-musllinux_1_1_x86_64.whl", hash = "sha256:18d90f4711bf63e2fb21e8c8e51ed8189438e6b35a6d996201ebd98a26abbbe6"},
+    {file = "cryptography-35.0.0-cp36-abi3-win32.whl", hash = "sha256:c10c797ac89c746e488d2ee92bd4abd593615694ee17b2500578b63cad6b93a8"},
+    {file = "cryptography-35.0.0-cp36-abi3-win_amd64.whl", hash = "sha256:7075b304cd567694dc692ffc9747f3e9cb393cc4aa4fb7b9f3abd6f5c4e43588"},
+    {file = "cryptography-35.0.0-pp36-pypy36_pp73-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:a688ebcd08250eab5bb5bca318cc05a8c66de5e4171a65ca51db6bd753ff8953"},
+    {file = "cryptography-35.0.0-pp36-pypy36_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d99915d6ab265c22873f1b4d6ea5ef462ef797b4140be4c9d8b179915e0985c6"},
+    {file = "cryptography-35.0.0-pp36-pypy36_pp73-manylinux_2_24_x86_64.whl", hash = "sha256:928185a6d1ccdb816e883f56ebe92e975a262d31cc536429041921f8cb5a62fd"},
+    {file = "cryptography-35.0.0-pp37-pypy37_pp73-macosx_10_10_x86_64.whl", hash = "sha256:ebeddd119f526bcf323a89f853afb12e225902a24d29b55fe18dd6fcb2838a76"},
+    {file = "cryptography-35.0.0-pp37-pypy37_pp73-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:22a38e96118a4ce3b97509443feace1d1011d0571fae81fc3ad35f25ba3ea999"},
+    {file = "cryptography-35.0.0-pp37-pypy37_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:eb80e8a1f91e4b7ef8b33041591e6d89b2b8e122d787e87eeb2b08da71bb16ad"},
+    {file = "cryptography-35.0.0-pp37-pypy37_pp73-manylinux_2_24_x86_64.whl", hash = "sha256:abb5a361d2585bb95012a19ed9b2c8f412c5d723a9836418fab7aaa0243e67d2"},
+    {file = "cryptography-35.0.0-pp37-pypy37_pp73-win_amd64.whl", hash = "sha256:1ed82abf16df40a60942a8c211251ae72858b25b7421ce2497c2eb7a1cee817c"},
+    {file = "cryptography-35.0.0.tar.gz", hash = "sha256:9933f28f70d0517686bd7de36166dda42094eac49415459d9bdf5e7df3e0086d"},
+]
 cycler = [
     {file = "cycler-0.11.0-py3-none-any.whl", hash = "sha256:3a27e95f763a428a739d2add979fa7494c912a32c17c4c38c4d5f082cad165a3"},
     {file = "cycler-0.11.0.tar.gz", hash = "sha256:9c87405839a19696e837b3b818fed3f5f69f16f1eec1a1ad77e043dcea9c772f"},
@@ -665,6 +1010,14 @@ greenlet = [
     {file = "greenlet-1.1.2-cp39-cp39-win_amd64.whl", hash = "sha256:013d61294b6cd8fe3242932c1c5e36e5d1db2c8afb58606c5a67efce62c1f5fd"},
     {file = "greenlet-1.1.2.tar.gz", hash = "sha256:e30f5ea4ae2346e62cedde8794a56858a67b878dd79f7df76a0767e356b1744a"},
 ]
+html5lib = [
+    {file = "html5lib-1.1-py2.py3-none-any.whl", hash = "sha256:0d78f8fde1c230e99fe37986a60526d7049ed4bf8a9fadbad5f00e22e58e041d"},
+    {file = "html5lib-1.1.tar.gz", hash = "sha256:b2e5b40261e20f354d198eae92afc10d750afb487ed5e50f9c4eaf07c184146f"},
+]
+idna = [
+    {file = "idna-3.3-py3-none-any.whl", hash = "sha256:84d9dd047ffa80596e0f246e2eab0b391788b0503584e8945f2368256d2735ff"},
+    {file = "idna-3.3.tar.gz", hash = "sha256:9d643ff0a55b762d5cdb124b8eaa99c66322e2157b69160bc32796e824360e6d"},
+]
 importlib-metadata = [
     {file = "importlib_metadata-4.8.2-py3-none-any.whl", hash = "sha256:53ccfd5c134223e497627b9815d5030edf77d2ed573922f7a0b8f8bb81a1c100"},
     {file = "importlib_metadata-4.8.2.tar.gz", hash = "sha256:75bdec14c397f528724c1bfd9709d660b33a4d2e77387a3358f20b848bb5e5fb"},
@@ -673,6 +1026,14 @@ iniconfig = [
     {file = "iniconfig-1.1.1-py2.py3-none-any.whl", hash = "sha256:011e24c64b7f47f6ebd835bb12a743f2fbe9a26d4cecaa7f53bc4f35ee9da8b3"},
     {file = "iniconfig-1.1.1.tar.gz", hash = "sha256:bc3af051d7d14b2ee5ef9969666def0cd1a000e121eaea580d4a313df4b37f32"},
 ]
+jeepney = [
+    {file = "jeepney-0.7.1-py3-none-any.whl", hash = "sha256:1b5a0ea5c0e7b166b2f5895b91a08c14de8915afda4407fb5022a195224958ac"},
+    {file = "jeepney-0.7.1.tar.gz", hash = "sha256:fa9e232dfa0c498bd0b8a3a73b8d8a31978304dcef0515adc859d4e096f96f4f"},
+]
+keyring = [
+    {file = "keyring-23.2.1-py3-none-any.whl", hash = "sha256:bd2145a237ed70c8ce72978b497619ddfcae640b6dcf494402d5143e37755c6e"},
+    {file = "keyring-23.2.1.tar.gz", hash = "sha256:6334aee6073db2fb1f30892697b1730105b5e9a77ce7e61fca6b435225493efe"},
+]
 kiwisolver = [
     {file = "kiwisolver-1.3.2-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:1d819553730d3c2724582124aee8a03c846ec4362ded1034c16fb3ef309264e6"},
     {file = "kiwisolver-1.3.2-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:8d93a1095f83e908fc253f2fb569c2711414c0bfd451cab580466465b235b470"},
@@ -760,6 +1121,9 @@ matplotlib = [
     {file = "matplotlib-3.5.0-pp37-pypy37_pp73-win_amd64.whl", hash = "sha256:a0ea10faa3bab0714d3a19c7e0921279a68d57552414d6eceaea99f97d7735db"},
     {file = "matplotlib-3.5.0.tar.gz", hash = "sha256:38892a254420d95594285077276162a5e9e9c30b6da08bdc2a4d53331ad9a6fa"},
 ]
+mimeparse = [
+    {file = "mimeparse-0.1.3.tar.gz", hash = "sha256:534ff6feefe1cd03984f444e6415aacc79c0a85f3b868ec41a2fd5003004c31e"},
+]
 mirakuru = [
     {file = "mirakuru-2.4.1-py3-none-any.whl", hash = "sha256:097324abe7479b3e6a8b745d0e3980664c8f6d3aec06cdeff8fc38cb2c9abc93"},
     {file = "mirakuru-2.4.1.tar.gz", hash = "sha256:7025d121d2f04e957bd6ae3239531d60ebba787839c89d6052479beb58b0cd0b"},
@@ -982,6 +1346,10 @@ py = [
 py-cpuinfo = [
     {file = "py-cpuinfo-8.0.0.tar.gz", hash = "sha256:5f269be0e08e33fd959de96b34cd4aeeeacac014dd8305f70eb28d06de2345c5"},
 ]
+pycparser = [
+    {file = "pycparser-2.21-py2.py3-none-any.whl", hash = "sha256:8ee45429555515e1f6b185e78100aea234072576aa43ab53aefcae078162fca9"},
+    {file = "pycparser-2.21.tar.gz", hash = "sha256:e644fdec12f7872f86c58ff790da456218b10f863970249516d60a5eaca77206"},
+]
 pyerfa = [
     {file = "pyerfa-2.0.0.1-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:278832de7803f2fb0ef4b14263200f98dfdb3eaa78dc63835d93796fd8fc42c6"},
     {file = "pyerfa-2.0.0.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:629248cebc8626a52e80f69d4e2f30cc6e751f57803f5ba7ec99edd09785d181"},
@@ -1025,6 +1393,10 @@ pytest-benchmark = [
     {file = "pytest-benchmark-3.4.1.tar.gz", hash = "sha256:40e263f912de5a81d891619032983557d62a3d85843f9a9f30b98baea0cd7b47"},
     {file = "pytest_benchmark-3.4.1-py2.py3-none-any.whl", hash = "sha256:36d2b08c4882f6f997fd3126a3d6dfd70f3249cde178ed8bbc0b73db7c20f809"},
 ]
+pytest-doctestplus = [
+    {file = "pytest-doctestplus-0.11.0.tar.gz", hash = "sha256:efcb24a366d3b033d343298c2a60eae418b06803c85373450f10306815c7c632"},
+    {file = "pytest_doctestplus-0.11.0-py3-none-any.whl", hash = "sha256:90e5aff08072aaea1221745482087b02d4851221a7342e539faa2370ed02b1bd"},
+]
 pytest-postgresql = [
     {file = "pytest-postgresql-3.1.2.tar.gz", hash = "sha256:f3582a51506b0aa2dd292bfd12bac5f27982ef9235b516140437d7af8e4b7287"},
     {file = "pytest_postgresql-3.1.2-py3-none-any.whl", hash = "sha256:32a15caa8cbebb7e8cbccb1c03bda55a00e7de9a2409add8b04dfbfa60f6432b"},
@@ -1037,6 +1409,22 @@ pytz = [
     {file = "pytz-2021.3-py2.py3-none-any.whl", hash = "sha256:3672058bc3453457b622aab7a1c3bfd5ab0bdae451512f6cf25f64ed37f5b87c"},
     {file = "pytz-2021.3.tar.gz", hash = "sha256:acad2d8b20a1af07d4e4c9d2e9285c5ed9104354062f275f3fcd88dcef4f1326"},
 ]
+pyvo = [
+    {file = "pyvo-1.1-py3-none-any.whl", hash = "sha256:a7a9b4f93b62d23e3d82837aa95877d8dd05e58b31078e08650fd4fa38972491"},
+    {file = "pyvo-1.1.tar.gz", hash = "sha256:e14f527ecbec48b71d57f48f301ec6a809f886096114529cf2d1fdcbaf6deb2b"},
+]
+pywin32-ctypes = [
+    {file = "pywin32-ctypes-0.2.0.tar.gz", hash = "sha256:24ffc3b341d457d48e8922352130cf2644024a4ff09762a2261fd34c36ee5942"},
+    {file = "pywin32_ctypes-0.2.0-py2.py3-none-any.whl", hash = "sha256:9dc2d991b3479cc2df15930958b674a48a227d5361d413827a4cfd0b5876fc98"},
+]
+requests = [
+    {file = "requests-2.26.0-py2.py3-none-any.whl", hash = "sha256:6c1246513ecd5ecd4528a0906f910e8f0f9c6b8ec72030dc9fd154dc1a6efd24"},
+    {file = "requests-2.26.0.tar.gz", hash = "sha256:b8aa58f8cf793ffd8782d3d8cb19e66ef36f7aba4353eec859e74678b01b07a7"},
+]
+secretstorage = [
+    {file = "SecretStorage-3.3.1-py3-none-any.whl", hash = "sha256:422d82c36172d88d6a0ed5afdec956514b189ddbfb72fefab0c8a1cee4eaf71f"},
+    {file = "SecretStorage-3.3.1.tar.gz", hash = "sha256:fd666c51a6bf200643495a04abb261f83229dcb6fd8472ec393df7ffc8b6f195"},
+]
 setuptools-scm = [
     {file = "setuptools_scm-6.3.2-py3-none-any.whl", hash = "sha256:4c64444b1d49c4063ae60bfe1680f611c8b13833d556fd1d6050c0023162a119"},
     {file = "setuptools_scm-6.3.2.tar.gz", hash = "sha256:a49aa8081eeb3514eb9728fa5040f2eaa962d6c6f4ec9c32f6c1fba88f88a0f2"},
@@ -1045,6 +1433,10 @@ six = [
     {file = "six-1.16.0-py2.py3-none-any.whl", hash = "sha256:8abb2f1d86890a2dfb989f9a77cfcfd3e47c2a354b01111771326f8aa26e0254"},
     {file = "six-1.16.0.tar.gz", hash = "sha256:1e61c37477a1626458e36f7b1d82aa5c9b094fa4802892072e49de9c60c4c926"},
 ]
+soupsieve = [
+    {file = "soupsieve-2.3.1-py3-none-any.whl", hash = "sha256:1a3cca2617c6b38c0343ed661b1fa5de5637f257d4fe22bd9f1338010a1efefb"},
+    {file = "soupsieve-2.3.1.tar.gz", hash = "sha256:b8d49b1cd4f037c7082a9683dfa1801aa2597fb11c3a1155b7a5b94829b4f1f9"},
+]
 sqlalchemy = [
     {file = "SQLAlchemy-1.4.27-cp27-cp27m-macosx_10_14_x86_64.whl", hash = "sha256:6afa9e4e63f066e0fd90a21db7e95e988d96127f52bfb298a0e9bec6999357a9"},
     {file = "SQLAlchemy-1.4.27-cp27-cp27m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:ec1c908fa721f2c5684900cc8ff75555b1a5a2ae4f5a5694eb0e37a5263cea44"},
@@ -1095,6 +1487,14 @@ typing-extensions = [
     {file = "typing_extensions-4.0.0-py3-none-any.whl", hash = "sha256:829704698b22e13ec9eaf959122315eabb370b0884400e9818334d8b677023d9"},
     {file = "typing_extensions-4.0.0.tar.gz", hash = "sha256:2cdf80e4e04866a9b3689a51869016d36db0814d84b8d8a568d22781d45d27ed"},
 ]
+urllib3 = [
+    {file = "urllib3-1.26.7-py2.py3-none-any.whl", hash = "sha256:c4fdf4019605b6e5423637e01bc9fe4daef873709a7973e195ceba0a62bbc844"},
+    {file = "urllib3-1.26.7.tar.gz", hash = "sha256:4987c65554f7a2dbf30c18fd48778ef124af6fab771a377103da0585e2336ece"},
+]
+webencodings = [
+    {file = "webencodings-0.5.1-py2.py3-none-any.whl", hash = "sha256:a0af1213f3c2226497a97e2b3aa01a7e4bee4f403f95be16fc9acd2947514a78"},
+    {file = "webencodings-0.5.1.tar.gz", hash = "sha256:b36a1c245f2d304965eb4e0a82848379241dc04b865afcc4aab16748587e1923"},
+]
 zipp = [
     {file = "zipp-3.6.0-py3-none-any.whl", hash = "sha256:9fe5ea21568a0a70e50f273397638d39b03353731e6cbbb3fd8502a33fec40bc"},
     {file = "zipp-3.6.0.tar.gz", hash = "sha256:71c644c5369f4a6e07636f0aa966270449561fcea2e3d6747b8d23efaa9d7832"},
diff --git a/pyproject.toml b/pyproject.toml
index 6565a79..d072f97 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -16,17 +16,31 @@ classifiers = [
 python = "^3.7"
 astropy = "*"
 astropy_healpix = "*"
-mocpy = ">=0.9"
+mocpy = ">=0.10.0"
 sqlalchemy = "*"
 
 [tool.poetry.dev-dependencies]
+astroquery = "*"
 numpy = "*"
 pandas = "*"
 psycopg2-binary = "*"
 pytest = "*"
 pytest-benchmark = "*"
+pytest-doctestplus = "*"
 pytest-postgresql = "<4.0.0"  # pytest-postgresql 4.0.0 uses psycopg3, but sqlalchemy does not yet support it
 
+[tool.pytest.ini_options]
+addopts = "--doctest-glob *.md"
+doctest_plus = true
+doctest_optionflags = [
+    "FLOAT_CMP"
+]
+filterwarnings = [
+    "error::sqlalchemy.exc.SAWarning",
+    "error::DeprecationWarning",
+    "error::PendingDeprecationWarning"
+]
+
 [build-system]
 requires = ["poetry>=0.12"]
 build-backend = "poetry.masonry.api"


Test patch:
diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml
index 0b02ea6..b0a230c 100644
--- a/.github/workflows/build-and-test.yml
+++ b/.github/workflows/build-and-test.yml
@@ -50,17 +50,17 @@ jobs:
     - name: Test with pytest
       run: |
         poetry run pytest --benchmark-group-by=fullfunc --benchmark-sort=fullname --benchmark-json benchmark_results.json
-    - name: Create runtime plots
-      if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'}}
-      run: |
-        mkdir plots
-        poetry run python .github/workflows/runtime_plots.py
-    - name: Upload plots to github artifacts
-      if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'}}
-      uses: actions/upload-artifact@v2
-      with:
-        name: runtime-plots
-        path: plots
+    # - name: Create runtime plots
+    #   if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'}}
+    #   run: |
+    #     mkdir plots
+    #     poetry run python .github/workflows/runtime_plots.py
+    # - name: Upload plots to github artifacts
+    #   if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'}}
+    #   uses: actions/upload-artifact@v2
+    #   with:
+    #     name: runtime-plots
+    #     path: plots
     - name: Upload benchmark results to github artifacts
       if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.9'}}
       uses: actions/upload-artifact@v2
diff --git a/conftest.py b/conftest.py
new file mode 100644
index 0000000..edaf1f3
--- /dev/null
+++ b/conftest.py
@@ -0,0 +1,19 @@
+"""Pytest configuration for running the doctests in README.md."""
+from unittest.mock import Mock
+import sqlalchemy as sa
+import pytest
+
+
+@pytest.fixture
+def engine(postgresql):
+    """Create an SQLAlchemy engine with a disposable PostgreSQL database."""
+    return sa.create_engine('postgresql://', poolclass=sa.pool.StaticPool,
+                            creator=lambda: postgresql)
+
+
+@pytest.fixture(autouse=True)
+def add_mock_create_engine(monkeypatch, request):
+    """Monkey patch sqlalchemy.create_engine for doctests in README.md."""
+    if request.node.name == 'README.md':
+        engine = request.getfixturevalue('engine')
+        monkeypatch.setattr(sa, 'create_engine', Mock(return_value=engine))
diff --git a/healpix_alchemy/tests/__init__.py b/healpix_alchemy/tests/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/healpix_alchemy/tests/conftest.py b/healpix_alchemy/tests/conftest.py
deleted file mode 100644
index f817afb..0000000
--- a/healpix_alchemy/tests/conftest.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import pytest
-from sqlalchemy import create_engine
-from sqlalchemy.orm import sessionmaker
-from sqlalchemy.pool import StaticPool
-
-
-@pytest.fixture
-def engine(postgresql):
-    """Create an SQLAlchemy engine with a disposable PostgreSQL database."""
-    return create_engine('postgresql://',
-                         poolclass=StaticPool,
-                         creator=lambda: postgresql)
-
-
-@pytest.fixture
-def session(engine, record_database_size):
-    """Create an SQLAlchemy session with a disposable PostgreSQL database."""
-    Session = sessionmaker(bind=engine)
-    session = Session()
-    yield session
-    record_database_size()
-    session.close()
-
-
-@pytest.fixture
-def record_database_size(record_property, engine):
-
-    def func():
-        (database_size,), = engine.execute(
-            '''SELECT pg_size_pretty(sum(pg_relation_size(C.oid)))
-               FROM pg_class C
-               LEFT JOIN pg_namespace N ON (N.oid = C.relnamespace)
-               WHERE nspname NOT IN ('pg_catalog', 'information_schema')''')
-        record_property('database_size', database_size)
-
-    return func
-
-
-def pytest_terminal_summary(terminalreporter):
-    terminalreporter.section('database size')
-    for report in terminalreporter.getreports(''):
-        try:
-            database_size = dict(report.user_properties)['database_size']
-        except KeyError:
-            pass
-        else:
-            terminalreporter.line(f'{report.nodeid}: {database_size}')
diff --git a/healpix_alchemy/tests/test_regions.py b/healpix_alchemy/tests/test_regions.py
deleted file mode 100644
index b3cf1e8..0000000
--- a/healpix_alchemy/tests/test_regions.py
+++ /dev/null
@@ -1,122 +0,0 @@
-import re
-
-from astropy.coordinates import SkyCoord
-from astropy.table import Table
-from astropy import units as u
-from astropy.utils.data import get_readable_fileobj
-from astropy_healpix import level_ipix_to_uniq
-from mocpy import MOC, IntervalSet
-import numpy as np
-from sqlalchemy.ext.declarative import declarative_base
-from sqlalchemy import Column, Integer, ForeignKey
-from sqlalchemy.orm import relationship
-
-from .. import Region, Tile
-
-
-Base = declarative_base()
-
-
-class Field(Region, Base):
-    __tablename__ = 'fields'
-
-    field_id = Column(
-        Integer,
-        primary_key=True)
-
-    tile_class = lambda: FieldTile  # noqa: E731
-
-    tiles = relationship("FieldTile")
-
-
-class FieldTile(Tile, Base):
-    __tablename__ = 'fieldtiles'
-
-    instrument_field_id = Column(
-        ForeignKey('fields.field_id', ondelete="CASCADE"),
-        nullable=False,
-        doc='Instrument Field ID',
-    )
-
-    field = relationship(
-        "Field",
-        foreign_keys=instrument_field_id
-    )
-
-
-def get_ztf_footprint_corners():
-    """Return the corner offsets of the ZTF footprint."""
-    x = 6.86 / 2
-    return [-x, +x, +x, -x] * u.deg, [-x, -x, +x, +x] * u.deg
-
-
-def get_footprints_grid(lon, lat, offsets):
-    """Get a grid of footprints for an equatorial-mount telescope.
-
-    Parameters
-    ----------
-    lon : astropy.units.Quantity
-        Longitudes of footprint vertices at the standard pointing.
-        Should be an array of length N.
-    lat : astropy.units.Quantity
-        Latitudes of footprint vertices at the standard pointing.
-        Should be an array of length N.
-    offsets : astropy.coordinates.SkyCoord
-        Pointings for the field grid.
-        Should have length M.
-
-    Returns
-    -------
-    astropy.coordinates.SkyCoord
-        Footprints with dimensions (M, N).
-    """
-    lon = np.repeat(lon[np.newaxis, :], len(offsets), axis=0)
-    lat = np.repeat(lat[np.newaxis, :], len(offsets), axis=0)
-    return SkyCoord(lon, lat, frame=offsets[:, np.newaxis].skyoffset_frame())
-
-
-def test_fields(request, session, engine):
-    """Generate ZTF fields and ingest them."""
-    Base.metadata.create_all(engine)
-
-    url = ('https://github.com/ZwickyTransientFacility/ztf_information/raw/'
-           'master/field_grid/ZTF_Fields.txt')
-    with get_readable_fileobj(url, show_progress=False) as f:
-        first_line, *lines = f
-        names = re.split(r'\s\s+', first_line.lstrip('%').strip())
-        table = Table.read(lines, format='ascii', names=names)
-    # keep subset of fields for speed
-    table = table[::100]
-
-    lon, lat = get_ztf_footprint_corners()
-    centers = SkyCoord(table['RA'] * u.deg, table['Dec'] * u.deg)
-    vertices = get_footprints_grid(lon, lat, centers)
-
-    mocs = [MOC.from_polygon_skycoord(verts) for verts in vertices]
-
-    fields = [Field.from_moc(moc, field_id=int(field_id))
-              for field_id, moc in zip(table['ID'], mocs)]
-    for field in fields:
-        session.add(field)
-    session.commit()
-
-    fields = session.query(Field).all()
-    assert len(fields) == 18
-
-
-def test_tile(engine):
-    """Generate an example Tile and confirm index calculation."""
-    Base.metadata.create_all(engine)
-
-    LEVEL = IntervalSet.HPX_MAX_ORDER
-
-    level, ipix = 3, 12
-    value = level_ipix_to_uniq(level, ipix)
-    shift = 2 * (LEVEL - level)
-    nested_lo = int(ipix << shift)
-    nested_hi = int(((ipix + 1) << shift) - 1)
-
-    tile = Tile(uniq=value)
-    assert tile.nested_lo == nested_lo
-    assert tile.nested_hi == nested_hi
-    assert tile.uniq == value