在学术论文中提及HEALPix Alchemy，请引用以下两个参考文献。

HEALPix Alchemy论文

Singer, L. P., Parazin, B., Coughlin, M. W., et al. (2022). "HEALPix Alchemy: Fast All-Sky Geometry and Image Arithmetic in a Relational Database for Multimessenger Astronomy Brokers." Astronomical Journal 163 209. https://doi.org/10.3847/1538-3881/ac5ab8

代码存档在Zenodo上

Singer, L. P., Parazin, B., Coughlin, M. W., et al. (2022). "skyportal/healpix-alchemy: Version 1.0.1." https://doi.org/10.5281/zenodo.5768564

HEALPix Alchemy

healpix_alchemy 是一个针对 Python 的扩展，用于扩展 SQLAlchemy 对象关系映射。它增加了区域和图像算术功能到 PostgreSQL（版本 14 及更高版本）数据库中。通过在 PostgreSQL 的 范围类型中存储 多级 HEALPix 索引，它加速了点云、区域和图像（有时在地理空间社区中称为栅格）之间的查询。

healpix_alchemy 项目是为天文学应用设计的，特别是用于交叉匹配星系目录、观测足迹以及全天空图像，如 引力波概率天空图 或甚至尘埃图。然而，它也可以在任何嵌入在单位球体上的几何形状的上下文中使用。

healpix_alchemy 精简且简约，因为它利用了几个现有的项目：它主要由几行粘合代码组成，用于绑定 MOCPy、SQLAlchemy 和 PostgreSQL 的 范围类型。

healpix_alchemy 的用途类似于具有完整功能的专注于天文学的空间扩展，如 Q3C、H3C 和 pg_healpix，以及地理空间扩展，如 PgSphere 和 PostGIS。与这些扩展的区别在于，healpix_alchemy 是用纯 Python 编写的，无需服务器端数据库扩展。因此，healpix_alchemy 可以与云中的托管 PostgreSQL 数据库（如 Amazon RDS 和 Google Cloud SQL）一起使用。

操作理论

HEALPix 基础知识

HEALPix 是一种将单位球体细分和索引的方案，最初由 Górski 等人 (2005) 描述。虽然它最初是为宇宙微波背景分析而设计的，但它已在天文学中找到了许多用途，尤其是在 多级覆盖 (MOC) 图和 分层渐进调查 (HiPS) 中，这些在 Aladin 天文信息系统中被广泛使用。它还被 LIGO 和 Virgo 用于存储和通信 引力波概率天空图。

HEALPix 可以被视为一棵树。在最低分辨率级别 0，HEALPix 将球体细分为 12 个相等的基板瓦片，分配整数索引 0 到 11。在级别 1，12 个基板瓦片中的每一个都被细分为 4 个瓦片。每个后续级别都将前一个级别的每个瓦片细分为 4 个新的瓦片。在给定的级别，每个基板像素都被细分为 4^级别 像素（每边的 nside = 2^级别 像素）。因此，在给定的分辨率下，有 npix = 12×4^级别 像素，分配整数索引从 0 到 (npix-1)。这被称为 NESTED 索引方案。（还有一种 RING 索引方案，其中索引从东到西然后从北到南增加。）

HEALPix 瓦片，HEALPix 树的节点，可以通过三部分信息完全识别：索引方案（RING 或 NESTED）、分辨率级别（级别 或等价地 nside）和像素索引（ipix，一个介于 0 和 npix-1 之间的整数）。

下图中，从 https://healpix.jpl.nasa.gov 复制了 HEALPix 网格的第一级细化。

HEALPix 区间集

球面上的一个区域可以通过一组不相交的HEALPix瓦片进行编码，这些瓦片可能处于不同的分辨率级别。通常，在大区域的内部使用大低分辨率瓦片，而在边界处使用小高分辨率瓦片。这被称为多阶覆盖（MOC）图。以下是一个示例，它来自MOCPy文档。

与MOC类似，LIGO/Virgo/KAGRA引力波概率天空地图以多分辨率HEALPix数据集的形式存储，但每个瓦片都附加了一个浮点数值向量。以下是一个天空地图的多阶细化网格示例，它来自Singer & Price (2015)。

Reinecke & Hivon (2015)介绍了HEALPix区间集，作为MOC的替代编码方式，它能够实现快速简单的并集、交集和查询。在区间集中，每个HEALPix瓦片由一些非常高分辨率level_max的像素索引区间描述，这些像素索引是该瓦片的子节点。在区间集中，一个区域被编码为这种区间的不相交集合。具有NESTED地址（level，ipix）的瓦片可以描述为以下半开区间

[ipix 4^{level_max - level}，(ipix + 1) 4^{level_max - level}）。

我们使用level_max = 29，因为这是可以用有符号64位整数存储像素索引的最高分辨率。在这个分辨率下，每个像素的宽度仅为0.4毫角秒。

区间集表示的优点是，有简单的快速算法用于区间算术和集合运算。区间分析出现在从基因组学到引力波数据质量的多种科学背景下。由于区间算术的许多商业应用，区间也支持在PostgreSQL数据库中，通过它的范围类型。

healpix_alchemy中的空间原语

healpix_alchemy包为SQLAlchemy提供了两种自定义列类型

healpix_alchemy.Point

此类表示一个点。此类类型的列可以存储目录中星系的位置。底层上，它只是一个BIGINT。

在任何需要将Python值绑定到healpix_alchemy.Point的地方，您可以提供以下任何一个

• astropy.coordinates.SkyCoord的实例
• 两个具有角度单位的astropy.units.Quantity实例的序列，将解释为点在ICRS框架中的赤经和赤纬
• 代表在level = level_max处的点的HEALPix NESTED索引的整数

healpix_alchemy.Tile

此类表示一个HEALPix瓦片。包含此类类型列和外键的表可以存储MOC或引力波概率图。底层上，它只是一个INT8RANGE。

在任何需要将Python值绑定到healpix_alchemy.Tile的地方，您可以提供以下任何一个

• 一个整数，它将被解释为在UNIQ HEALPix索引方案中的瓦片地址
• 两个整数的序列，它将被解释为在level = level_max处的右半开像素索引区间的上下界
• 类似于'[1234,5678)'的字符串

安装

您可以使用pip从Python包索引安装healpix_alchemy

$ pip install healpix-alchemy

开发安装

欢迎贡献力量！本包使用Poetry打包和依赖工具以及pytest进行单元测试。要在开发环境中安装healpix_alchemy，请按照以下说明操作。

1. 按照官方Poetry安装说明安装Poetry。

2. 克隆此仓库

   $ git clone https://github.com/skyportal/healpix-alchemy.git
   $ cd healpix-alchemy
   

3. 通过运行以下命令初始化由Poetry管理的虚拟环境，其中已安装healpix_alchemy及其所有依赖项：

   $ poetry install
   

   现在，您可以通过运行以下命令进入虚拟环境内的shell：

   $ poetry shell
   

4. 要运行测试套件，包括此README文件中的示例，请在Poetry shell中运行以下命令：

   $ pytest
   

示例

首先，进行一些导入

>>> from sqlalchemy import orm
>>> import sqlalchemy as sa
>>> import healpix_alchemy as ha

设置表格

此示例将使用SQLAlchemy声明式扩展，通过Python类描述表模式。

SQLAlchemy需要知道每个表的名字。您可以在每个模型类中设置__tablename__属性来提供名字，或者创建一个基类，它会自动从类名生成表名。

>>> @orm.as_declarative()
... class Base:
...
...     @orm.declared_attr
...     def __tablename__(cls):
...         return cls.__name__.lower()

Galaxy表中的每一行代表目录中的一个点

>>> class Galaxy(Base):
...     id = sa.Column(sa.Text, primary_key=True)
...     hpx = sa.Column(ha.Point, index=True, nullable=False)

Field表中的每一行代表ZTF字段

>>> class Field(Base):
...     id = sa.Column(sa.Integer, primary_key=True)
...     tiles = orm.relationship(lambda: FieldTile)

FieldTile表中的每一行代表包含在相应字段内的多分辨率HEALPix瓦片。存在一个从Field到FieldTile的一对多映射。

>>> class FieldTile(Base):
...     id = sa.Column(sa.ForeignKey(Field.id), primary_key=True)
...     hpx = sa.Column(ha.Tile, primary_key=True, index=True)

Skymap表中的每一行代表LIGO/Virgo HEALPix定位图。

>>> class Skymap(Base):
...     id = sa.Column(sa.Integer, primary_key=True)
...     tiles = orm.relationship(lambda: SkymapTile)

SkymapTile表中的每一行代表LIGO/Virgo定位图内的多分辨率HEALPix瓦片。存在一个从Skymap到SkymapTile的一对多映射。

>>> class SkymapTile(Base):
...     id = sa.Column(sa.ForeignKey(Skymap.id), primary_key=True)
...     hpx = sa.Column(ha.Tile, primary_key=True, index=True)
...     probdensity = sa.Column(sa.Float, nullable=False)

最后，连接到数据库，创建所有表格，并启动会话。

>>> engine = sa.create_engine('postgresql://user:password@host/database')
>>> Base.metadata.create_all(engine)
>>> session = orm.Session(engine)

填充示例数据

将2MASS红移巡天加载到Galaxy表中。此目录包含44599个星系。

这可能需要长达一分钟的时间才能完成。高级用户可以通过将SkyCoord到HEALPix索引的转换矢量化并使用SQLAlchemy批量插入来显著加快此过程。

>>> 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'], hpx=row['coord']))
>>> session.commit()

将Zwicky瞬变设施字段的足迹加载到Field和FieldTile表中。

这可能需要长达一分钟的时间才能完成。高级用户可以通过使用SQLAlchemy批量插入来显著加快此过程。

>>> from astropy.table import Table
>>> from astropy.coordinates import SkyCoord
>>> from astropy import units as u
>>> url = 'https://raw.githubusercontent.com/ZwickyTransientFacility/ztf_information/9fd0ba8842709f42a134c88827309ccab728fcb7/field_grid/ztf_field_corners.csv'
>>> for row in Table.read(url):
...     field_id = int(row['field'])
...     corners = SkyCoord(row['ra1', 'ra2', 'ra3', 'ra4'],
...                        row['dec1', 'dec2', 'dec3', 'dec4'],
...                        unit=u.deg)
...     tiles = [FieldTile(hpx=hpx) for hpx in ha.Tile.tiles_from(corners)]
...     session.add(Field(id=field_id, tiles=tiles))
>>> session.commit()

将LIGO/Virgo事件GW200115_042309 (S200115j)的星空图加载到Skymap和SkymapTile表中。

>>> url = 'https://gracedb.ligo.org/apiweb/superevents/S200115j/files/bayestar.multiorder.fits'
>>> data = Table.read(url)
>>> tiles = [SkymapTile(hpx=row['UNIQ'], probdensity=row['PROBDENSITY']) for row in data]
>>> session.add(Skymap(id=1, tiles=tiles))
>>> session.commit()

最后，运行ANALYZE以准备数据供使用

>>> session.execute(sa.text('ANALYZE'))
<sqlalchemy.engine.cursor.CursorResult object at 0x...>

示例查询

每个字段的面积是多少？

>>> query = sa.select(
...     FieldTile.id, sa.func.sum(FieldTile.hpx.area)
... ).group_by(
...     FieldTile.id
... ).limit(
...     5
... )
>>> for id, area in session.execute(query):
...     print(f'Field {id} has area {area:.3g} sr')
Field 199 has area 0.0174 sr
Field 200 has area 0.0174 sr
Field 201 has area 0.0174 sr
Field 202 has area 0.0174 sr
Field 203 has area 0.0174 sr

每个字段中有多少星系？

>>> count = sa.func.count(Galaxy.id)
>>> query = sa.select(
...     FieldTile.id, count
... ).filter(
...     FieldTile.hpx.contains(Galaxy.hpx)
... ).group_by(
...     FieldTile.id
... ).order_by(
...     count.desc()
... ).limit(
...     5
... )
>>> for id, n in session.execute(query):
...     print(f'Field {id} contains {n} 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

每个星系位置的概率密度是多少？

>>> query = sa.select(
...     Galaxy.id, SkymapTile.probdensity
... ).filter(
...     SkymapTile.id == 1,
...     SkymapTile.hpx.contains(Galaxy.hpx)
... ).order_by(
...     SkymapTile.probdensity.desc()
... ).limit(
...     5
... )
>>> for id, p in session.execute(query):
...     print(f'{id} has prob. density {p:.5g}/sr')
2MASX J02532153+0632222 has prob. density 20.701/sr
2MASX J02530482+0555431 has prob. density 20.695/sr
2MASX J02533119+0628252 has prob. density 20.669/sr
2MASX J02524584+0639206 has prob. density 20.656/sr
2MASX J02534120+0615562 has prob. density 20.567/sr

每个字段包含的概率是多少？

>>> area = (FieldTile.hpx * SkymapTile.hpx).area
>>> prob = sa.func.sum(SkymapTile.probdensity * area)
>>> query = sa.select(
...     FieldTile.id, prob
... ).filter(
...     SkymapTile.id == 1,
...     FieldTile.hpx.overlaps(SkymapTile.hpx)
... ).group_by(
...     FieldTile.id
... ).order_by(
...     prob.desc()
... ).limit(
...     5
... )
>>> for id, prob in session.execute(query):
...     print(f'Field {id} probability is {prob:.3g}')
Field 1499 probability is 0.165
Field 1446 probability is 0.156
Field 452 probability is 0.154
Field 505 probability is 0.0991
Field 401 probability is 0.0962

字段1000到2000的总面积是多少？

在接下来的两个示例中，我们介绍了healpix_alchemy.func.union()，该函数用于找到一组瓦片的并集。因为它是一个聚合函数，通常应将其用于子查询中。

>>> union = sa.select(
...     ha.func.union(FieldTile.hpx).label('hpx')
... ).filter(
...     FieldTile.id.between(1000, 2000)
... ).subquery()
>>> query = sa.select(
...     sa.func.sum(union.columns.hpx.area)
... )
>>> result = session.execute(query).scalar_one()
>>> print(f'{result:.3g} sr')
9.33 sr

字段1000到2000中包含的积分概率是多少？

>>> union = sa.select(
...     ha.func.union(FieldTile.hpx).label('hpx')
... ).filter(
...     FieldTile.id.between(1000, 2000)
... ).subquery()
>>> prob = sa.func.sum(SkymapTile.probdensity * (union.columns.hpx * SkymapTile.hpx).area)
>>> query = sa.select(
...     prob
... ).filter(
...     SkymapTile.id == 1,
...     union.columns.hpx.overlaps(SkymapTile.hpx)
... )
>>> result = session.execute(query).scalar_one()
>>> print(f'{result:.3g}')
0.837

90%置信区域的面积是多少？

>>> cum_area = sa.func.sum(
...     SkymapTile.hpx.area
... ).over(
...     order_by=SkymapTile.probdensity.desc()
... ).label(
...     'cum_area'
... )
>>> cum_prob = sa.func.sum(
...     SkymapTile.probdensity * SkymapTile.hpx.area
... ).over(
...     order_by=SkymapTile.probdensity.desc()
... ).label(
...     'cum_prob'
... )
>>> subquery = sa.select(
...     cum_area,
...     cum_prob
... ).filter(
...     SkymapTile.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(f'{result:.3g} sr')
0.277 sr

哪些星系位于90%置信区域内？

>>> cum_prob = sa.func.sum(
...     SkymapTile.probdensity * SkymapTile.hpx.area
... ).over(
...     order_by=SkymapTile.probdensity.desc()
... ).label(
...     'cum_prob'
... )
>>> subquery = sa.select(
...     SkymapTile.probdensity,
...     cum_prob
... ).filter(
...     SkymapTile.id == 1
... ).subquery()
>>> min_probdensity = sa.select(
...     sa.func.min(subquery.columns.probdensity)
... ).filter(
...     subquery.columns.cum_prob <= 0.9
... ).scalar_subquery()
>>> query = sa.select(
...     Galaxy.id
... ).filter(
...     SkymapTile.id == 1,
...     SkymapTile.hpx.contains(Galaxy.hpx),
...     SkymapTile.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