ガイア計画というものを知りました。
欧州宇宙機関 (ESA) が打ち上げた専用の探査機で全天球の天体を観測してデータを蓄積するプロジェクトで、現在のデータ量は1.8TB、観測された天体の数は約1,900億個だそうです。このデータを利用して全天球の星の地図を作ってみました。
データのダウンロード
ガイア計画のデータは以下でアーカイブされていて、自由に利用できます。日本にもミラーサイトがあり国立天文台のウェブサイト jvo.nao.ac.jp からデータをダウンロードできます。今回は作業時点で最も新しい EDR3 (Early Data Research 3) を使いました。
なお、データの利用にあたっては所定のクレジットを併記する必要があり、専用ページで確認できます。
データの加工
GAIA データは gzip 圧縮された単純な CSV としてダウンロードできます。データには天球上での星の位置が含まれていて、それは黄道面を基準に黄経(ecl_lon)、黄緯(ecl_lat)として与えられているようです。今回は web 地図の上に天球を投影しようと思いましたので、赤道面を基準にした位置情報（赤経・赤緯）に変換しました。
この変換は単純な球座標上での回転移動とみなして以下の1次変換でOKということにしました。赤経をα、赤緯をβ、黄経δ、黄緯をγ、黄道傾斜角（地球の自転軸の傾き）をεとすると、以下の連立方程式で (δ, γ) から (α, β) が計算できるはずです。左辺は赤道面を基準にしたxyz座標、右辺は回転行列と黄道面を基準にしたxyz座標の積です。
\begin{eqnarray*}
\left(
\begin{array}{c}
\cos\beta\cos\alpha \\
\cos\beta\sin\alpha \\
\sin\beta \\
\end{array}
\right)
=\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos\epsilon & -\sin\epsilon \\
0 & \sin\epsilon & \cos\epsilon \\
\end{array}
\right)
\left(
\begin{array}{c}
\cos\delta\cos\gamma \\
\cos\delta\sin\gamma \\
\sin\delta \\
\end{array}
\right)
\end{eqnarray*}
これを解くと、
\begin{equation}
\alpha=\mathrm{tan}^{-1}\frac{\cos\epsilon\cos\delta\sin\gamma-\sin\epsilon\sin\delta}{\cos\delta\cos\gamma}
\end{equation}
\begin{equation}
\beta=\mathrm{sin}^{-1}(\sin\epsilon\cos\delta\sin\gamma+\cos\epsilon\sin\delta)
\end{equation}
なお、今回は簡便に黄道傾斜角（地球の自転軸の傾き）として ε = 23.4° を使いました。
注意 結果として後述しますが、この射影変換が間違っているようで、意図した通り天球が平面に展開されませんでした。この件は別途検証しようと思います。
ベクトルタイルの生成
データが次のフォーマットの順に加工されるようにパイプして mbtiles を生成します。
1. csv.gz // 元データ。gzip 圧縮された csv
2. csv // csv
3. GeoJSON // 位置情報フォーマット
4. NDGeoJSONs // 3 を 改行区切りで結合したもの
5. mbtiles // ベクトルタイル
今回は 2 の csv から 3 のGeoJSON を生成するする単純なストリームの処理のみを自作しました。csv のカラムには例えば phot_g_mean_mag のような等級（星の明るさ）などの膨大な観測データがあるのですが、地学の素養がなくデータの意味を十分に理解できませんでした..。このストリームでは、 phot_g_mean_mag は m プロパティとして保存しています。この値は結局使いませんでした。
// bin/build.js
const csv = require('csv-parser')
// 定数と計算用ユーティリティ
const ECL_INC = 23.4
const ang2rad = (ang) => ang * Math.PI / 180
const rad2ang = (rad) => 180 * rad / Math.PI
const sin = (ang) => Math.sin(ang2rad(ang))
const cos = (ang) => Math.cos(ang2rad(ang))
const atan = (rad) => rad2ang(Math.atan(rad))
const asin = (rad) => rad2ang(Math.asin(rad))
// 黄道 -> 赤道の座標変換
const ecl2eq = ([lng, lat]) => {
const eq_lng = atan(cos(lng) * cos(lat) / (cos(ECL_INC) * cos(lng) * sin(lat) - sin(ECL_INC) * sin(lng)))
const eq_lat = asin(sin(ECL_INC) * cos(lng) * sin(lat) + cos(ECL_INC) * sin(lng))
return [eq_lng, eq_lat]
}
process.stdin
.pipe(csv())
.on('data', (chunk) => {
const { ecl_lon, ecl_lat, phot_g_mean_mag } = chunk
const m = parseFloat(phot_g_mean_mag)
const ecl_coord = [parseFloat(ecl_lon), parseFloat(ecl_lat)]
if(!ecl_coord[0] || !ecl_coord[1] !m) {
return
}
try {
const eq_coord = ecl2eq(ecl_coord);
const feature = {
type: 'Feature',
properties: { m },
geometry: { type: 'Point', coordinates: eq_coord },
}
process.stdout.write(JSON.stringify(feature) + '\n')
} catch (error) {
process.stderr.write(`[failure] ${JSON.stringify(error)}\n`)
}
})
tippeacanoe は標準入力をサポートしているので、例えば以下のようなコマンドで簡単に mbtiles を生成できます。1つのファイルあたりに約50万個〜90万個ほどの星のデータが入っているようです。
$ find ~/Downloads/GaiaSource_*.csv.gz | \
xargs -I {} sh -c "cat {} | gunzip | node bin/build.js" | \
tippecanoe -zg -o gaia.mbtiles --drop-densest-as-needed
なお、実際には S3 にダウンロードしたデータを用いて以下のようなコマンドで複数に分けて処理しました。
# NUMBER_PREFIX は 0 から 7
$ aws s3api list-objects --bucket $GAIA_DOWNLOAD_BUCKET_NAME | \
jq .Contents | jq -r 'map(.Key)[]' | \
grep "GaiaSource_$NUMBER_PREFIX" | \
xargs -t -I {} sh -c "aws s3 cp s3://$GAIA_DOWNLOAD_BUCKET_NAME/{} - | gunzip | node bin/build.js" | \
tippecanoe -zg -o gaia$NUMBER_PREFIX.mbtiles -l gaia --drop-densest-as-needed
生成した 8 つの mbtiles は tippeacanoe の tile-join コマンドで結合します。
$ tile-join -o gaia.mbtiles \
gaia0.mbtiles \
gaia1.mbtiles \
gaia2.mbtiles \
gaia3.mbtiles \
gaia4.mbtiles \
gaia5.mbtiles \
gaia6.mbtiles \
gaia7.mbtiles --no-tile-size-limit
ベクトルタイルをホストする
tileserver-gl でベクトルタイルを配信します。以下のコマンドで http://localhost:8080/data/gaia.json からベクトルタイルが利用できるようになります。
$ nvm use v10
$ npx tileserver-gl ./gaia.mbtiles
ベクトルタイルをレンダリングする
タイルを読み込んでレンダリングするための html とスタイルを作成します。
私が属している Geolonia の Embed API を使うと mapboxgl/maplibregl と互換の Map クラスを利用できます。URL に緯度経度をシリアライズする data-hash="on" などの属性が使えるので便利です。
<html>
<head>
<style>
html, body, #map {
width: 100%;
height: 100%;
margin: 0;
padding: 0;
}
</style>
</head>
<body>
<div id="map" data-hash="on"></div>
<script src="https://cdn.geolonia.com/v1/embed?geolonia-api-key=YOUR-API-KEY"></script>
<script>
new window.geolonia.Map({
container: document.getElementById('map'),
style: {
version: 8,
sources: {
gaia: {
type: 'vector',
url: 'http://localhost:8080/data/gaia.json',
},
},
layers: [
{
id: 'background',
type: 'background',
paint: {
'background-color': "black"
}
},
{
id: 'stars',
type: 'circle',
source: 'gaia',
'source-layer': 'gaia',
paint: {
'circle-color': '#ffffff',
'circle-radius': 1
}
}
]
}
})
</script>
</body>
</html>
結果と考察
3,384個、613GB の GZIP されたCSV から8つの mbtiles を生成するのに、 AWS EC2 の t3.2xlarge インスタンスを使った処理で概ね24時間かかりました。全ての mbtiles を結合したデータのサイズは 13.5 GB でした。
星の海を探検してみる
どうも射影変換が間違っているようで、 -45° から + 45°の範囲にしか星が分布しませんでした。ただ全単射ではありますので、このまま星の分布を眺めたいと思います。
北緯66°と南緯66°付近に集まるラグビーボールの縫い目のような線が出現しました。66°というのはおそらく 90° から黄道傾斜角を引いた値からくるもので、線は8つに分割したデータを tile-join することで出現しているもののように思えます。
ところどころ不自然な角張を持つ疎な領域や密な領域があります。これは元データの欠損や重複に由来しているもののように見えます。
星の集まった銀河のような構造も見えます。
星の密度が薄い場所があります。これは暗黒物質でしょうか。
まとめ
今回使った Gaia EDR3 のようなテラバイト級のデータを手軽にブラウジングするのは大変です。 解決方法として、ベクトルタイルで逐次データをロードできるようにするのは有効なのではないでしょうか。
この記事に関連するソースは以下にあります。
謝辞
This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
このエントリでは、欧州宇宙機関（ESA）のミッション Gaia のデータを利用しており、これは、 Gaia Data Processing and Analysis Consortium (DPAC) が処理を行なったものです。 DPAC への資金は、各国の機関、特に Gaia Multilateral Agreement に参加している機関から提供されています。