MapLibre GL JS で「夜」を描いた話

## はじめに

<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/816106c8-9494-4297-c98a-76aa04ae6f84.png">
<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/03854d08-e3ce-4707-e51f-5ed72f5b90c6.png">

https://github.com/kikuchan/maplibre-gl-nightlayer

[MapLibre GL JS](https://github.com/maplibre/maplibre-gl-js) 用の、[夜を描画するレイヤー](https://github.com/kikuchan/maplibre-gl-nightlayer)を実装したので、その紹介と解説記事です。

## カスタムレイヤー機能

まずは MapLibre GL JS で使えるカスタムレイヤー機能について解説します。
この機能を使うと、地図の上に独自の WebGL レイヤーを重ねることができます。

### 使い方

`CustomLayerInterface` を実装したオブジェクトを `Map` の `addLayer` に渡すことで、カスタムレイヤーとして機能させることができます。

具体的には `type` フィールドが `"custom"` な、(`id` フィールドがある、) `render` 関数を持つオブジェクトです。
```js
{ id: 'my-custom-layer', type: 'custom', render: function ... }
```
レイヤーの描画が必要になるたびに `render` 関数が呼び出されます。
また `onAdd`, `onRemove`, `prerender` という名前の関数があれば、適時呼び出されます。

ところで `render` 関数と`prerender` 関数に渡される引数は下記のような定義になっているのですが、MapLibre GL JS の v4 と v5 で少し異なっています。 詳しくは後で見ます。

:::note info
本記事の執筆時点ではまだ v5 の正式版がリリースされていないため、`5.0.0-pre.10` 時点での定義を使用しています。 また、見やすさのため型を若干省略しています。（以後同様）
:::

```ts
// v4 interface
render(gl: WebGLRenderingContext, matrix: mat4, options: CustomRenderMethodInput): void;
// v5 interface
render(gl: WebGLRenderingContext, options: CustomRenderMethodInput): void;
```


さて、そのようなオブジェクトを素直にクラスとして TypeScript で実装すると以下のような見た目になります。

```ts
import type { CustomLayerInterface, CustomRenderMethodInput, Map as MaplibreMap } from 'maplibre-gl';

class MyCustomLayer implements CustomLayerInterface {
  id = 'my-custom-layer';
  type = 'custom' as const;

      :

  onAdd(map: MaplibreMap, gl: WebGLRenderingContext) {
    // マップに追加されたときの処理
  }

  onRemove(map: MaplibreMap, gl: WebGLRenderingContext) {
    // マップから削除されたときの処理
  }

  // v5 用
  render(gl: WebGLRenderingContext, options: CustomRenderMethodInput) {
    // gl を使って絵を描く
  }
```

最後に `Map` の `addLayer` APIを使ってこれを追加します。

```ts
const map = new maplibregl.Map({
    :
});

map.on('load', () => {
  map.addLayer(new MyCustomLayer())
});
```

これで、作成したクラス[^object]がカスタムレイヤーとして機能します。
[^object]: もちろん素のオブジェクトのままでも良いです。お好みでどうぞ！

### 絵の描き方

先ほどの `render` 関数に WebGL のコンテキストが渡されてくるので、それを使って自由に、伸びやかに、**気合で**絵を描きましょう！

ただし地図として描く際の座標変換が少し自明ではないので、そこを詳しく見ていきます。

#### v4 の場合
```ts
render(gl: WebGLRenderingContext, matrix: mat4, options: CustomRenderMethodInput): void
```

2番目の引数 `matrix` がカメラ変換行列なので、これを用いて座標を変換します。

`MercatorCoordinate.fromLngLat` などを使って緯度経度を内部表現（0～1のメルカトル座標）に変換した後、`matrix` を適用することで画面表示用の座標に変換できます。

```glsl
gl_Position = u_matrix * vec4(a_position, 0., 1.);
```

v4 の場合はこれだけなので、大変シンプルです。


#### v5 の場合
```ts
render(gl: WebGLRenderingContext, options: CustomRenderMethodInput): void
```

地球儀モードなどの場合に単純な変換だけでは表現できなくなるためか、(現時点での最新版では) `matrix` 引数が削除されています。

代わりに、`options` (の型である `CustomRenderMethodInput`) に `shaderData` と `defaultProjectionData` が追加されています。

```ts
      :
    shaderData: {
        variantName: string;
        vertexShaderPrelude: string;
        define: string;
    };
    defaultProjectionData: ProjectionData;
};
```

まず手順としては `shaderData` の内容が変わっていないかを確認します。 確認のためのキーとして `variantName` を使うことが意図されています。

```ts
if (!this.#programCache.has(shaderData.variantName)) {
    // シェーダを作る
}
```
内容が変わっていれば、その内容に応じたシェーダを作成しなければなりません。

具体的には、頂点シェーダに `vertexShaderPrelude` と `define` を直接埋め込みます。 そうすることで、それぞれの**投影方法に固有な** `projectTile` などのヘルパー関数群が自動的に定義されます。
```ts
const vertexSource = `#version 300 es
    precision highp float;
    
    ${shaderData.vertexShaderPrelude}
    ${shaderData.define}
    
      :
`;
```

次に進みます。

通常の手順に加え、ヘルパー関数が内部で使う uniform 変数を、`defaultProjectionData` の値に従って適切に設定する必要があります。

下記は `test/examples/globe-custom-simple.html` から拝借したものです。
```ts
gl.uniformMatrix4fv(
    gl.getUniformLocation(program, 'u_projection_fallback_matrix'),
    false,
    defaultProjectionData.fallbackMatrix // convert mat4 from gl-matrix to a plain array
);
gl.uniformMatrix4fv(
    gl.getUniformLocation(program, 'u_projection_matrix'),
    false,
    defaultProjectionData.mainMatrix // convert mat4 from gl-matrix to a plain array
);
gl.uniform4f(
    gl.getUniformLocation(program, 'u_projection_tile_mercator_coords'),
    ...defaultProjectionData.tileMercatorCoords
);
gl.uniform4f(
    gl.getUniformLocation(program, 'u_projection_clipping_plane'),
    ...defaultProjectionData.clippingPlane
);
gl.uniform1f(
    gl.getUniformLocation(program, 'u_projection_transition'),
    defaultProjectionData.projectionTransition
);
```

これでやっと座標変換のための準備が整いました[^program]。

[^program]: ライブラリに用意されている `Program` クラスなどを使って組み立てると適切にやってくれます。 しかし使い方が少し複雑だったので今回は直接行いました。

（v4 の時と同様）`MercatorCoordinate.fromLngLat` などを使って緯度経度を内部表現（0～1のメルカトル座標）に変換した後、ヘルパー関数 `projectTile` に渡すことで画面表示用の座標に変換できます。

```glsl
gl_Position = projectTile(a_position);
```


:::note info

地球儀モードの場合には `projectTile` が非線形な変換になるため、引き渡す頂点データはあらかじめ細かくメッシュ化[^mesh]しておく必要があります。
:::
[^mesh]: メッシュ化されたタイル座標頂点データを生成する `createTileMesh` 関数が新たに追加されていますので、これを使うと少しだけ楽ができます。

## 「夜」の実装

次に、「**夜**」の実装について解説します。

夜というのは太陽が作る**地球の影**ですから、v5 で導入予定の地球儀モードにおいては、ライティングの設定で比較的簡単に実現できる[^lighting]表現です。

[^lighting]: 現状のベータ版では `sky` を指定するとデフォルトでは何故かほぼ夜の面しか見えませんが、例えば後述の太陽直下点（とこの [`lightPosition`](https://gist.github.com/kikuchan/3a0127b8eb19bcf49ad0a271d64b2ef0) 関数）を用いて次のようにすると良い感じになります: `{ light: { anchor: 'map', position: lightPosition(subsolar.lat, subsolar.lng) } }` （ただし無理やり感がすごいので、たぶんそのうち修正・変更されるでしょう）

しかし、一般的な2Dの地図上に表現しようとすると様々な理由で少し難しくなります。

でもどうしても夜の部分をいい感じに表現したかったので、前述のカスタムレイヤー機能を使って独自のレイヤー（**NightLayer** / **夜レイヤー**）を実装しました。

### そもそも「夜」とは

夕日が沈んだからと言って、すぐに完全な暗闇になるわけではもちろんありません。
この中間の時間帯は、黄昏時・かわたれ時・**薄明**などと呼ばれています。

…では、いつになったら「夜」と言えるのでしょうか？

実は、太陽が地平線に沈んだ後の、見えなくなった太陽の**中心**と地平線のなす角度（**俯角**）を用いて以下のように定義されています。

| 俯角 | 英語 | 日本語 |
|:-:|-|-|
| (日出/日没)[^sunset] ～ 6 | civil twilight | 市民薄明 |
| 6 ～ 12 | nautical twilight | 航海薄明 |
| 12 ～ 18 | astronomical twilight | 天文薄明 |
| 18 ～  | night | 夜 |

(参考: [薄明 - 国立天文台](https://eco.mtk.nao.ac.jp/koyomi/topics/html/topics1992.html))

これを図示するとこのようになります。

<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/f92a3ba6-ad30-aed5-05b0-5cd72c78612d.png" /> <img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/4f59c6a2-d6a8-5c78-7863-4f1715f5e2ed.png" />
(引用: https://en.wikipedia.org/wiki/Twilight)

よって定義上は俯角が18度以上になったところからが「夜」ということになります。

### 計算・実装方法

:::note warn
あくまでも簡易的な近似計算なので、精度も含めて如何なる保証もしません。
:::

まず[**太陽高度**](https://ja.wikipedia.org/wiki/%E5%A4%AA%E9%99%BD%E9%AB%98%E5%BA%A6)を考えます。太陽の中心と地平線のなす角度を、**地平線を0度**として（先ほどとは逆に）**見上げる方向をプラス**（仰角）とします。 頭上にあれば、90度です。

地図上のすべての点（ピクセル）それぞれに対応する地点（観測地点）での**太陽高度**を求め、**プラス**であれば（日中なので）**完全な透明**に、**マイナス**[^minus]であれば（薄明～夜なので）その度合いに応じた**黒**をいい感じの透明度で描画すれば、既存の地図に重畳しやすい**半透明の暗闇レイヤー**（夜レイヤー）ができあがります。

[^minus]: 厳密には、マイナスであっても太陽が完全に沈んでいなければ日中です。 しかしここでは簡単のため、太陽高度が0度になる（太陽が半分出ている）ところを市民薄明の開始として扱っています。

[^sunset]: 太陽の視半径や大気の屈折などを考慮に入れると、およそ 0.833 度です。

#### 太陽の位置の計算

太陽高度を求めるために、まず、太陽の位置から求めます。

ここでは、太陽の位置を表現するために[**太陽直下点**](https://ja.wikipedia.org/wiki/%E5%A4%AA%E9%99%BD%E7%9B%B4%E4%B8%8B%E7%82%B9)（[**Subsolar Point**](https://en.wikipedia.org/wiki/Subsolar_point)）を使います。太陽直下点とは太陽が**完全に真上**に来る**地球上の座標**のことで、緯度経度で表します。

計算式には https://en.wikipedia.org/wiki/Equation_of_time#Alternative_calculation を使います。 この近似方法は、従来の複雑な近似手法に比べてシンプルな上に、精度も若干良いそうです。

```ts
function getSubsolarPoint(date: Date) {
    // based on https://en.wikipedia.org/wiki/Equation_of_time#Alternative_calculation
    const D = (date.getTime() - Date.UTC(date.getUTCFullYear(), 0, 0)) / 86400000;
    const n = (2 * Math.PI) / 365.24;
    
    const e = radians(23.44); // Earth's axial tilt
    const E = 0.0167; // Earth's orbital eccentricity
    
    const A = (D + 9) * n;
    const B = A + 2 * E * Math.sin((D - 3) * n);
    const C = (A - Math.atan2(Math.sin(B), Math.cos(B) * Math.cos(e))) / Math.PI;
    
    const EOT = 720 * (C - Math.trunc(C + 0.5));
    
    const UTC = date.getUTCHours() + date.getUTCMinutes() / 60 + date.getUTCSeconds() / 3600;
    
    const lng = -15 * (UTC - 12 + EOT / 60);
    const lat = degrees(Math.asin(Math.sin(-e) * Math.cos(B)));
    
    return {
        lng,
        lat,
    };
}
```

#### 太陽高度の計算

太陽直下点が確定した後は、その情報をもとにそれぞれの観測地点での太陽高度を求めます。 **それぞれの観測地点＝地図上の各ピクセル**なので、計算はフラグメントシェーダ内で行います。 ここがカスタムレイヤー（≒ WebGL）であることの利点です。

```glsl
float A = sin(observer.y) * sin(subsolar.y);
float B = cos(observer.y) * cos(subsolar.y) * cos(subsolar.x - observer.x);
float altitude = degrees(asin(A + B));
```

#### 暗闇の計算

太陽高度が求まったので、実際の各ピクセルの色を計算します。
透明度は、薄明の定義を参考に、**6度ごとに下絵の輝度が半分になる**ように設定します。

```glsl
float twilightLevel = -altitude / 6.0;
float darkness = (1. - clamp(pow(0.5, twilightLevel), 0., 1.));

gl_FragColor = vec4(u_color / 255., 1.0) * darkness;
```

これをレンダリングするとこのようになります。

<img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/6a45fd6c-39e9-8d7f-25e2-1d8f46f90754.png" /> <img width="350" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3961633/745ed581-60bc-6975-14d1-6cedd8ecb6c0.png" />

### 完成

以上を NightLayer としてクラス化し、パッケージ化したものがこちらです。[^openlayers]

[^openlayers]: PoC 止まりで少し古いですが、実は [OpenLayers版](https://github.com/openlayers/openlayers/pull/16028) もあります。

https://github.com/kikuchan/maplibre-gl-nightlayer

ついでに色々と機能[^features]を追加してあります。

[^features]: 冒頭右側にあるような薄明の段階表現や、色や透明度の指定などが可能です。


## おわりに

以上、駆け足でカスタムレイヤーの使い方と、夜レイヤーの実装について見てきました。

カスタムレイヤー機能を使うことで、お手軽に WebGL の独自レイヤーを重ねることができる、というのが伝わっていたら嬉しいです。

実は今回この記事を書くために頑張って [maplibre-gl-nightlayer](https://github.com/kikuchan/maplibre-gl-nightlayer) を v4 / v5 両対応にさせました。 そのため、ここで例示したコードより実際のものは少しだけ複雑になっています。 少々泥臭いですが、気になる方はぜひコードの方も覗いてみてください。


地球儀モードが目玉となる v5 の正式リリースが待ち遠しいですね！