ベクトルタイルを使って星の海を探検する



ガイア計画というものを知りました。

https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A4%E3%82%A2%E8%A8%88%E7%94%BB

https://www.cosmos.esa.int/web/gaia/home

欧州宇宙機関　(ESA)　が打ち上げた専用の探査機で全天球の天体を観測してデータを蓄積するプロジェクトで、現在のデータ量は1.8TB、観測された天体の数は約1,900億個だそうです。このデータを利用して全天球の星の地図を作ってみました。　

### データのダウンロード

ガイア計画のデータは以下でアーカイブされていて、自由に利用できます。日本にもミラーサイトがあり国立天文台のウェブサイト jvo.nao.ac.jp からデータをダウンロードできます。今回は作業時点で最も新しい [EDR3](https://www.cosmos.esa.int/web/gaia/earlydr3) (Early Data Research 3) を使いました。

https://gea.esac.esa.int/archive/

http://jvo.nao.ac.jp/portal/gaia/edr3.do?action=gotoDownloadPage

なお、データの利用にあたっては所定のクレジットを併記する必要があり、[専用ページ](https://gea.esac.esa.int/archive/documentation/GEDR3/Miscellaneous/sec_credit_and_citation_instructions/)で確認できます。

### データの加工

GAIA データは gzip 圧縮された単純な CSV としてダウンロードできます。データには天球上での星の位置が含まれていて、それは黄道面を基準に黄経(`ecl_lon`)、黄緯(`ecl_lat`)として与えられているようです。今回は web 地図の上に天球を投影しようと思いましたので、赤道面を基準にした位置情報（赤経・赤緯）に変換しました。

<img width=600 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/f911e89d-5d05-110b-f1bb-a3534e2dcc73.png">

この変換は単純な球座標上での回転移動とみなして以下の1次変換でOKということにしました。赤経をα、赤緯をβ、黄経δ、黄緯をγ、黄道傾斜角（地球の自転軸の傾き）をεとすると、以下の連立方程式で (δ, γ) から (α, β) が計算できるはずです。左辺は赤道面を基準にしたxyz座標、右辺は回転行列と黄道面を基準にしたxyz座標の積です。

```math
\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*}
```

これを解くと、

```math
\begin{equation}
\alpha=\mathrm{tan}^{-1}\frac{\cos\epsilon\cos\delta\sin\gamma-\sin\epsilon\sin\delta}{\cos\delta\cos\gamma}
\end{equation}
```

```math
\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` プロパティとして保存しています。この値は結局使いませんでした。

```javascript
// 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万個ほどの星のデータが入っているようです。

```shell
$ 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 にダウンロードしたデータを用いて以下のようなコマンドで複数に分けて処理しました。

```shell
# 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 コマンドで結合します。

```shell
$ 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](https://github.com/maptiler/tileserver-gl) でベクトルタイルを配信します。以下のコマンドで　`http://localhost:8080/data/gaia.json` からベクトルタイルが利用できるようになります。

```shell
$ nvm use v10
$ npx tileserver-gl ./gaia.mbtiles
```

### ベクトルタイルをレンダリングする

タイルを読み込んでレンダリングするための html とスタイルを作成します。
私が属している [Geolonia](https://geolonia.com) の [Embed API](https://docs.geolonia.com/embed-api/) を使うと mapboxgl/maplibregl と互換の Map クラスを利用できます。URL に緯度経度をシリアライズする `data-hash="on"` などの属性が使えるので便利です。


```html
<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°の範囲にしか星が分布しませんでした。ただ全単射ではありますので、このまま星の分布を眺めたいと思います。

![スクリーンショット 2021-12-02 13.41.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/68a79abc-6136-2425-c618-301cedb58d0c.png)

北緯66°と南緯66°付近に集まるラグビーボールの縫い目のような線が出現しました。66°というのはおそらく 90° から黄道傾斜角を引いた値からくるもので、線は8つに分割したデータを tile-join することで出現しているもののように思えます。

![スクリーンショット 2021-12-02 13.42.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/28b01a39-04d2-f4b8-f495-025b53f988c7.png)

ところどころ不自然な角張を持つ疎な領域や密な領域があります。これは元データの欠損や重複に由来しているもののように見えます。

![スクリーンショット 2021-12-02 13.42.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/5c26ce64-9444-e5bb-028f-0f8e1945da48.png)

星の集まった銀河のような構造も見えます。

![スクリーンショット 2021-12-02 13.43.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/6987413c-56a6-3612-2e1f-9577a3ec864e.png)

星の密度が薄い場所があります。これは暗黒物質でしょうか。

![スクリーンショット 2021-12-02 13.43.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/125062/4a324185-f1d6-8315-5e3c-250ded364cbb.png)

#### まとめ

今回使った Gaia EDR3 のようなテラバイト級のデータを手軽にブラウジングするのは大変です。 解決方法として、ベクトルタイルで逐次データをロードできるようにするのは有効なのではないでしょうか。

この記事に関連するソースは以下にあります。

https://github.com/kamataryo/gaia-map

##### 謝辞

This work has made use of data from the European Space Agency (ESA) mission [Gaia](https://www.cosmos.esa.int/gaia), processed by the [Gaia Data Processing and Analysis Consortium](https://www.cosmos.esa.int/web/gaia/dpac/consortium) (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. 

このエントリでは、欧州宇宙機関（ESA）のミッション [Gaia](https://www.cosmos.esa.int/gaia) のデータを利用しており、これは、 [Gaia Data Processing and Analysis Consortium](https://www.cosmos.esa.int/web/gaia/dpac/consortium) (DPAC) が処理を行なったものです。 DPAC への資金は、各国の機関、特に Gaia Multilateral Agreement に参加している機関から提供されています。