全球データ可視化手法HiPSの地理空間への応用
はじめに
日常的に使用されている球面座標系は北極・南極域に画素が集中するため、一般的なWeb地図で用いられているWebメルカトル（Pseudo-Mercator）法（EPSG 3857）では北南緯85.06度以上を切り捨てています。
天文学の全球データを適切に扱うHEALPix/HiPSを地理空間データに応用することで、地理空間情報を全球で適切に扱うことができるようになります。
HiPSとは
HiPS (Hierarchical Progressive Survey) は、天文学データをWeb地図のように自由に拡大・縮小して可視化するためのデータ配信規格で、HEALPix投影法を使用することで北極・南極域を含む全球の領域を破綻なく投影します。
HiPS実例
Aladin LiteというWeb星図でHiPS星図/地図を見ることができます。
- 月全球地図: "Moon LROC-WAC-100m" progressive survey
- 天の川全天星図: "Gaia DR3 color flux map" progressive survey
- その他HiPSタイル一覧
HEALPixとは
Figure 2: Orthographic view of HEALPix partition of the sphere. - "The HEALPix Primer"
HEALPixは1997年に考案された地図投影法で、全球を H×K（一般的にはH=4, K=3の 12 面）の領域に分割し、その領域内を一辺
この投影法は分割された各画素で面積が等しいという特徴があり、全球のデータを等しく扱うことに適しています。
HEALPix分割
HEALPixは「Nside」と呼ばれる分割サイズと、「画素（Pixel）」と呼ばれる領域番号の二組の数字(Nside, Pix)で領域を特定します。
まず基本分割として、球面を経度方向に4分割、緯度方向に3分割の12面に分割します（H=4,K=3が一般的）。
それぞれの「面（Face）」には[0, 11] の番号が割り振られます（この時、4の左半分は0度線をまたいでいることに注意）。
0 1 2 3
4 5 6 7
8 9 10 11
続いて各面を一辺
この
ある「Nside」での全球の総画素数（総分割数）は 12 * Nside * Nsideで求められ、以下のように変化します。
- Nside=1: 12
- Nside=2: 48
- Nside=4: 192
- Nside=8: 3072
- Nside=16: 12288
- ...
Nside は1 << zoom（zoom=[0, ∞)）の式で求められ、HiPS規格では「Norder」と呼ばれています。これは地理空間のXYZタイルにおける「Z」に相当します。
分割された画素にはZ階数曲線に従って領域番号が振られます。
Figure 3: Cylindrical projection of the HEALPix division - "The HEALPix Primer"
Nside=2 & Nside=4
これにより、ひとつ上位・下位の解像度の領域番号は次のように関連します。
- 下位の領域番号を 2bit 右シフトすると、ひとつ上位の解像度の領域番号が得られる
- 上位の領域番号を 2bit 左シフトし
[0, 3]を加えると、ひとつ下位の解像度の領域番号が得られる
Figure 1: Quadrilateral tree pixel numbering. - "The HEALPix Primer"
HEALPixの領域番号はビットレベルで以下のような構成になっており、下位の解像度の番号からビット演算のみで、任意の上位解像度の領域番号を求めることができます。
Norder=4 (Nside=16)
faceyxyxyxyx
↓ (pix >> 4)
Norder=2 (Nside=4)
faceyxyx
↓（pix >> 2）
Norder=1 (Nside=2)
faceyx
↓（pix >> 2）
Norder=0 (Nside=1)
face
HEALPixの領域番号のビット列（faceを除く）は次のような構成になります。
(Nside=4)
Y
^
+----+----+----+----+
3 |1010|1011|1110|1111|
+----+----+----+----+
2 |1000|1001|1100|1101|
+----+----+----+----+
1 |0010|0011|0110|0111|
+----+----+----+----+
0 |0000|0001|0100|0101|
+----+----+----+----+ > X
O 0 1 2 3
この時、原点Oは各Faceの菱形の左の頂点に対応し、X軸は右下方向の辺、Y軸は右上方向の辺に対応します（上図は左に45度倒した状態に等しい）。
領域番号がXとYの情報が交互に織り交ぜられていることを思い出すと、ビット列を一つ飛ばしに読み出すと、原点（菱形の左の頂点）からのベクトルの長さが得られます。
YXYX X X
0b1101 => 1 1 => 0b11 => 3
Y Y
=> 1 0 => 0b10 => 2
実際に0b1101のタイルを確認すると、(3, 2)の画素であることが確認できます。
HEALPix実装
経緯度からHEALPixの領域番号を求めるオープンソース実装が複数公開されています。
healpix
healpix (GPL): https://healpix.sourceforge.io/
NASA JPLの実装。最適化が行われているらしいですが、ライセンスがGPLであることに注意してください。
const auto zoom = 0;
const auto nside = 1 << zoom;
const auto theta = std::numbers::pi / 2 - radians(latitude);
const auto phi = radians(longitude);
const auto pix = ang2pix_nest(nside, theta, phi);
healpix_bare
healpix_bare (BSD): https://sourceforge.net/projects/healpix/files/healpix_bare_1.0/
オリジナル実装と思われます。BSDライセンスでの提供のため残されているとのこと。
ひどく遅いというわけではないので本実装でも不都合はありません。
const auto zoom = 0;
const auto nside = 1 << zoom;
const auto theta = std::numbers::pi / 2 - radians(latitude);
const auto phi = radians(longitude);
const auto pix = ang2nest(nside, t_ang { theta, phi });
HiPS画像
HEALPixの画素番号から「面」「X」「Y」の情報を抜き出すと「X」「Y」を画像の座標として各面での画像が得られます。
HiPS画像は同様の方式で画像を作成しますが、座標のとり方がHEALPixの座標系とは異なっていることに注意が必要です。
具体的にはHEALPixの「X」をYに、「Y」をXでプロットした画像（つまりHEALPix座標系を転置したもの）がHiPS画像になります。
from "HiPS – Hierarchical Progressive Survey Version 1.0"
この転置の背景はビット列をfaceyxyxyx...と見るか、facexyxyxy...と見るかの違いで、前者はレンダリングした画像が倒立像になりますが、後者は正立像になるようです。
どうやら、生成された画像が違和感のないものにするための転置のように思えます。
本記事でHEALPixの座標系を前者faceyxyxyx...としたのは参照実装healpixのint nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num)や、healpix_bareのstruct t_hpd { int64_t x, y; int32_t f; };をもとにしています。
HiPS 階層構造
HiPS標準では以下のようにファイルを配置します。
Norder{K}/Dir{D}/Npix{N}.{ext}
-
K: Norder値（XYZタイルのZに相当） -
D:int(N / 10000) * 10000した値（1フォルダ1万ファイルに収める意図がある様子） -
N: HEALPix 領域番号 -
ext: 拡張子
地理空間情報への応用の提案
HiPSの考え方は地理空間系で従来から用いられていたXYZタイルの延長線上と見ることができ、現行の地理空間情報技術と非常に親和性が高いです。
また、HEALPixの(face, x, y)はタイル内座標と見ることができ、Mapbox Vector Tileのようなベクトルタイル（HiPS Vector Tile？）の作成に応用することが考えられます。
メリット
地理空間タイルをHiPS、あるいはHiPSを応用したデータ配信手法にすることで次のようなメリットがあります。
- 北極・南極域の破綻がなくなる
- 衛星コンステレーションによる全球地理空間データを極域も含めて適切に処理できる
- 画素間の等面積性を利用したデータ解析が可能になる
- 緯度差によって面積が異ならないので、記録されたデータを素直に解析できる
デメリット
一方デメリットとしては、ズームレベル0（Norder=0, Nside=1）の時点でタイル数が12枚なので、従来のXYZタイルに対して同一ズームレベルで12倍のファイルが生成されてしまいます。
（そのためかHiPSではDir{D}で1フォルダ10000ファイルに収める工夫を凝らしています。）
ラスタータイルであれば、Cloud Optimized GeoTIFFを応用して、一つのTIFFにHiPSタイルを収めることができるでしょう。
Cloud Optimized HiPS
以下がサンプル実装です。
TIFF画像を以下のようにタイリングします。
- ImageWidth: 4×Nside×タイルサイズ
- ImageLength: 3×Nside×タイルサイズ
- TileWidth: タイルサイズ
- TileLength: タイルサイズ
ここでHiPSは画像が転置されていることに注意してNorder=0の表示と同じになるようにタイルを並べると以下のような順に並びます。
z=1 (nside=2)
+----+----+----+----+----+----+----+----+
| 0 | 2 | 4 | 6 | 8 | 10 | 12 | 14 |
+----+----+----+----+----+----+----+----+
| 1 | 3 | 5 | 7 | 9 | 11 | 13 | 15 |
+----+----+----+----+----+----+----+----+
| 16 | 18 | 20 | 22 | 24 | 26 | 28 | 30 |
+----+----+----+----+----+----+----+----+
| 17 | 19 | 21 | 23 | 25 | 27 | 29 | 31 |
+----+----+----+----+----+----+----+----+
| 32 | 34 | 36 | 38 | 40 | 42 | 44 | 46 |
+----+----+----+----+----+----+----+----+
| 33 | 35 | 37 | 39 | 41 | 43 | 45 | 47 |
+----+----+----+----+----+----+----+----+
images: "Gaia DR3 color flux map" progressive survey
※512x512のタイル画像を使用した場合、Norder=5以上でuint16_tの最大値を超えるので、上記サンプル実装はBigTIFFを用いています。
Discussion