TypeScript + Mapbox-GL-JS⑥ 軌跡の可視化 - Space Time Cube編

シリーズ一覧 (更新:2023/2/9)
[① 環境構築 ~ 地図の表示](https://qiita.com/onk/items/864bf80d66afda9fb57e)
[② 地図のスタイルの変更](https://qiita.com/onk/items/362876cbce29dc2fdaf0)
[③ 円のプロット・ラベルの表示](https://qiita.com/onk/items/ba0cb9e79fd481ea0576)
[④ 軌跡の可視化 - 基本編](https://qiita.com/onk/items/7c528bf4d7791049666b)
[⑤ 軌跡の可視化 - アニメーション編](https://qiita.com/onk/items/e37232f52ea4d2244856)
[⑥ 軌跡の可視化 - Space Time Cube編](https://qiita.com/onk/items/77687712aec730e40046)

# 概要
時系列位置データの可視化手法の1つとしてSpace Time Cube(STC)というものがある。
地図の平面に対し垂直な方向に時間軸をとった3D空間上に、(緯度,経度,時刻)できまる点をプロットし、これを時系列順に結んだ軌跡で表現する。
参考：https://icaci.org/files/documents/ICC_proceedings/ICC2003/Papers/255.pdf

この記事ではMapbox-Gl-JSとThree.jsを用いた3次元表現(STC)の作り方を紹介する。
![スクリーンショット 2023-02-09 12.49.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2902914/a4b7d1f1-71f4-07ca-914f-399ae461117b.png)



# Mapbox-gl-js + Three.js
以下を参考にしました。
- https://docs.mapbox.com/jp/mapbox-gl-js/example/add-3d-model/
- https://github.com/hobby-overflow/yurui-mapboxgl-2object
- https://zenn.dev/hobby_overflow/articles/c0e7a796998484

## Three.jsの導入
インストール
```sh
npm i -S three @types/three
```
インポート
```ts
import * as THREE from 'three';
```

## Three.jsの基本設定
```ts
// Three.js のセットアップ
const camera: THREE.Camera = new THREE.Camera();
const scene: THREE.Scene = new THREE.Scene();
let renderer: THREE.WebGLRenderer = new THREE.WebGLRenderer();

// 光源の設定
const directionalLight = new THREE.DirectionalLight(0xffffff);
directionalLight.position.set(0, -70, 100).normalize();
scene.add(directionalLight);

const directionalLight2 = new THREE.DirectionalLight(0xffffff);
directionalLight2.position.set(0, 70, 100).normalize();
scene.add(directionalLight2);
```

## Three.jsを用いてMapbox-gl-jsに3Dオブジェクトを追加する
Mapbox上の座標系(緯度,経度,高度)からWebGL上の座標系(x,y,z)への座標変換を行う必要がある。
```ts
// 3Dモデルの座標変換のための変数
const modelRotate: THREE.Vector3 = new THREE.Vector3(Math.PI / 2, 0, 0);
const modelTranslate: mapboxgl.MercatorCoordinate = mapboxgl.MercatorCoordinate.fromLngLat(map.getCenter(), 0);
const modelScale: number = modelTranslate.meterInMercatorCoordinateUnits();
```

つづいて、WebGLのレイヤーを作成する
```ts
const customLayer: mapboxgl.AnyLayer = {
    id: '3d-model',
    type: 'custom',
    renderingMode: '3d',
    onAdd: function (_map: mapboxgl.Map, gl: WebGLRenderingContext) {
        renderer = new THREE.WebGLRenderer({
            canvas: map.getCanvas(),
            context: gl,
            antialias: true
        });
        renderer.autoClear = false;
    },
    render: function (_gl: WebGLRenderingContext, matrix: number[]) {
        const rotationX = new THREE.Matrix4().makeRotationAxis(
            new THREE.Vector3(1, 0, 0),
            modelRotate.x
        );
        const rotationY = new THREE.Matrix4().makeRotationAxis(
            new THREE.Vector3(0, 1, 0),
            modelRotate.y
        );
        const rotationZ = new THREE.Matrix4().makeRotationAxis(
            new THREE.Vector3(0, 0, 1),
            modelRotate.z
        );

        const m = new THREE.Matrix4().fromArray(matrix);
        const l = new THREE.Matrix4()
            .makeTranslation(
                Number(modelTranslate.x),
                Number(modelTranslate.y),
                Number(modelTranslate.z)
            )
            .scale(
                new THREE.Vector3(
                    modelScale,
                    -modelScale,
                    modelScale
                )
            )
            .multiply(rotationX)
            .multiply(rotationY)
            .multiply(rotationZ);

        camera.projectionMatrix = m.multiply(l);
        renderer.resetState();
        renderer.render(scene, camera);
        map.triggerRepaint();
    }
}
```

そして、レイヤーを追加する。
```ts
map.on('load', () => {
    map.addLayer(customLayer, 'waterway-label');
});
```

さらに、MeshをSceneに追加する。
```ts
map.on('load', () => {
    map.addLayer(customLayer, 'waterway-label');

    const material: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
        color: new THREE.Color('#f1ff57'),
        transparent: true,
        opacity: 0.5
    });
    const geometry: THREE.BoxGeometry = new THREE.BoxGeometry(radian);
    const box = new THREE.Mesh(geometry, material);
    scene.add(box);
});
```

## 3Dオブジェクトの座標変換
ここまでは3D空間のグローバル座標をMapboxの空間と一致させただけであり、ここの3DオブジェクトをMapboxの座標系へ変換することはできていない。できることなら、任意の緯度経度の位置に3Dオブジェクトを表示したいと考えていると、[同様の趣旨の記事](https://zenn.dev/hobby_overflow/articles/c0e7a796998484)があったので、[この方のGithub](https://github.com/hobby-overflow/yurui-mapboxgl-2object/tree/main/src/map)を参考にさせていただいた。

```ts
// MapBoxの座標(緯度経度)をWebGLの3D空間上の座標に変換する関数
function getLoaction(lng: number, lat: number): THREE.Vector3 {
    const lngLat = new mapboxgl.LngLat(lng, lat);
    const dist = map.getCenter().distanceTo(lngLat);
    const bearing = getBearing(map.getCenter(), lngLat);

    const x: number = dist * Math.sin((Math.PI * 2 * bearing) / 360);
    const z: number = dist * Math.cos((Math.PI * 2 * bearing) / 360);

    return new THREE.Vector3(x, 0, -z);
}

function getBearing(pos1: mapboxgl.LngLat, pos2: mapboxgl.LngLat): number {
    const radian = (deg: number): number => deg * Math.PI / 180;
    const degree = (rad: number): number => rad * 180 / Math.PI;

    let dLng = pos2.lng - pos1.lng;

    const x = Math.cos(radian(pos2.lat)) * Math.sin(radian(dLng));
    const z = Math.cos(radian(pos1.lat)) * Math.sin(radian(pos2.lat)) - Math.sin(radian(pos1.lat)) * Math.cos(radian(pos2.lat)) * Math.cos(radian(dLng));

    let bearingRadian = Math.atan2(x, z);
    let bearingDegree = degree(bearingRadian);

    if (bearingDegree > 0.0) return bearingDegree;
    return bearingDegree + 360;
}
```

# Space Time Cubeの実装
基本的なアルゴリズムは[④ 軌跡の可視化 - 基本編](https://qiita.com/onk/items/7c528bf4d7791049666b)で紹介したものと同じである。

MapboxのGeoJsonソースによる記述をThree.jsによる描画処理におきかえた。
```ts
const firstDate: number = new Date((<MyData>data[0]).date).valueOf();
const lastDate: number = new Date((<MyData>data[data.length - 1]).date).valueOf();

const radian: number = 2000;
const trajMaterial: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
    color: new THREE.Color('#f1ff57'),
    transparent: true,
    opacity: 0.5
});

for (let i = 0; i < data.length - 1; i++) {
    const d0: MyData = <MyData>data[i];
    const d1: MyData = <MyData>data[i + 1];

    const head: THREE.Vector3 = getLoaction(d0, firstDate, lastDate);
    const tail: THREE.Vector3 = getLoaction(d1, firstDate, lastDate);

    // 線分の描画
    const path: THREE.LineCurve3 = new THREE.LineCurve3(head, tail);
    const trajGeometry: THREE.TubeGeometry = new THREE.TubeGeometry(path, 1, radian, 32, false);
    const line: THREE.Mesh = new THREE.Mesh(trajGeometry, trajMaterial);
    scene.add(line);
}
```
さらに時刻をY座標に変換するために前出のgetLocation関数を書き換えた。
```ts
function getLoaction(d: MyData, fristDate: number, lastDate: number): THREE.Vector3 {
    const lngLat = new mapboxgl.LngLat(d.lng, d.lat);
    const dist = map.getCenter().distanceTo(lngLat);
    const bearing = getBearing(map.getCenter(), lngLat);

    const x: number = dist * Math.sin((Math.PI * 2 * bearing) / 360);
    const z: number = dist * Math.cos((Math.PI * 2 * bearing) / 360);

    // 変更箇所
    const date: number = new Date(d.date).valueOf();
    const ratio: number = (date - fristDate) / (lastDate - fristDate);
    const y: number = 400000 * ratio;
    // 変更箇所 (ここまで)

    return new THREE.Vector3(x, y, -z);
}
```

ここまででも十分ではあるが、線分の接点が綺麗でないことや、3D表現ゆえに地理的な座標が正確に読み取れないなどの欠点が生じる。
![スクリーンショット 2023-02-09 13.38.49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2902914/d79f0816-c31d-a5cc-cfcb-e3b0bb308fe7.png)
![スクリーンショット 2023-02-09 13.39.03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2902914/26d31f00-0972-a1a7-956a-f0160e33e4dd.png)


## 接点の描画
線と線の接点に球体を描画することで滑らかにした。
```ts
const trajMaterial: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
    color: new THREE.Color('#f1ff57'),
    transparent: true,
    opacity: 0.5
});
const sphereGeometry: THREE.SphereGeometry = new THREE.SphereGeometry(radian);
const joint0 = new THREE.Mesh(sphereGeometry, trajMaterial);
joint0.position.set(head.x, head.y, head.z)
scene.add(joint0);
const joint1 = new THREE.Mesh(sphereGeometry, trajMaterial);
joint1.position.set(tail.x, tail.y, tail.z)
scene.add(joint1);
```

## 基線の描画
軌跡の各点から地図へ下ろした垂線を描画することで、3D表現による弊害を解消する。
```ts
const axisMaterial: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
    color: new THREE.Color('#eeeeee'),
    transparent: true,
    opacity: 0.5
});

const bottom: THREE.Vector3 = new THREE.Vector3(tail.x, 0, tail.z);
const axisPath: THREE.LineCurve3 = new THREE.LineCurve3(tail, bottom);
const axisGeometry: THREE.TubeGeometry = new THREE.TubeGeometry(axisPath, 1, radian / 2, 32, false);
const axis: THREE.Mesh = new THREE.Mesh(axisGeometry, axisMaterial);
scene.add(axis);
```

## コードの全容
```ts
import mapboxgl from 'mapbox-gl';
import * as THREE from 'three';

mapboxgl.accessToken = /* your token */;

interface MyData {
    date: string, lng: number, lat: number
}


window.addEventListener('load', () => {
    const map: mapboxgl.Map = new mapboxgl.Map({
        container: 'map',
        style: 'mapbox://styles/mapbox/dark-v10',
        center: [139.7670516, 39],
        zoom: 5,
        pitch: 70,
        antialias: true
    });

    // 3Dモデルの座標変換のための変数
    const modelRotate: THREE.Vector3 = new THREE.Vector3(Math.PI / 2, 0, 0);
    const modelTranslate: mapboxgl.MercatorCoordinate = mapboxgl.MercatorCoordinate.fromLngLat(map.getCenter(), 0);
    const modelScale: number = modelTranslate.meterInMercatorCoordinateUnits();

    // Three.js のセットアップ
    const camera: THREE.Camera = new THREE.Camera();
    const scene: THREE.Scene = new THREE.Scene();
    let renderer: THREE.WebGLRenderer = new THREE.WebGLRenderer();

    // 光源の設定
    const directionalLight = new THREE.DirectionalLight(0xffffff);
    directionalLight.position.set(0, -70, 100).normalize();
    scene.add(directionalLight);

    const directionalLight2 = new THREE.DirectionalLight(0xffffff);
    directionalLight2.position.set(0, 70, 100).normalize();
    scene.add(directionalLight2);

    // 3Dモデルのレイヤー
    const customLayer: mapboxgl.AnyLayer = {
        id: '3d-model',
        type: 'custom',
        renderingMode: '3d',
        onAdd: function (_map: mapboxgl.Map, gl: WebGLRenderingContext) {
            renderer = new THREE.WebGLRenderer({
                canvas: map.getCanvas(),
                context: gl,
                antialias: true
            });
            renderer.autoClear = false;
        },
        render: function (_gl: WebGLRenderingContext, matrix: number[]) {

            const rotationX = new THREE.Matrix4().makeRotationAxis(
                new THREE.Vector3(1, 0, 0),
                modelRotate.x
            );
            const rotationY = new THREE.Matrix4().makeRotationAxis(
                new THREE.Vector3(0, 1, 0),
                modelRotate.y
            );
            const rotationZ = new THREE.Matrix4().makeRotationAxis(
                new THREE.Vector3(0, 0, 1),
                modelRotate.z
            );

            const m = new THREE.Matrix4().fromArray(matrix);
            const l = new THREE.Matrix4()
                .makeTranslation(
                    Number(modelTranslate.x),
                    Number(modelTranslate.y),
                    Number(modelTranslate.z)
                )
                .scale(
                    new THREE.Vector3(
                        modelScale,
                        -modelScale,
                        modelScale
                    )
                )
                .multiply(rotationX)
                .multiply(rotationY)
                .multiply(rotationZ);

            camera.projectionMatrix = m.multiply(l);
            renderer.resetState();
            renderer.render(scene, camera);
            map.triggerRepaint();
        }
    }

    map.on('load', () => {
        map.addLayer(customLayer, 'waterway-label');

        fetch('./trajectory.json')
            .then(res => res.json())
            .then(data => {
                const firstDate: number = new Date((<MyData>data[0]).date).valueOf();
                const lastDate: number = new Date((<MyData>data[data.length - 1]).date).valueOf();

                const radian: number = 2000;
                const trajMaterial: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
                    color: new THREE.Color('#f1ff57'),
                    transparent: true,
                    opacity: 0.5
                });
                const axisMaterial: THREE.MeshBasicMaterial = new THREE.MeshBasicMaterial({
                    color: new THREE.Color('#eeeeee'),
                    transparent: true,
                    opacity: 0.5
                });

                for (let i = 0; i < data.length - 1; i++) {
                    const d0: MyData = <MyData>data[i];
                    const d1: MyData = <MyData>data[i + 1];

                    const head: THREE.Vector3 = getLoaction(d0, firstDate, lastDate);
                    const tail: THREE.Vector3 = getLoaction(d1, firstDate, lastDate);

                    // 線分の描画
                    const path: THREE.LineCurve3 = new THREE.LineCurve3(head, tail);
                    const trajGeometry: THREE.TubeGeometry = new THREE.TubeGeometry(path, 1, radian, 32, false);
                    const line: THREE.Mesh = new THREE.Mesh(trajGeometry, trajMaterial);
                    scene.add(line);

                    // 基線の描画
                    const bottom: THREE.Vector3 = new THREE.Vector3(tail.x, 0, tail.z);
                    const axisPath: THREE.LineCurve3 = new THREE.LineCurve3(tail, bottom);
                    const axisGeometry: THREE.TubeGeometry = new THREE.TubeGeometry(axisPath, 1, radian / 2, 32, false);
                    const axis: THREE.Mesh = new THREE.Mesh(axisGeometry, axisMaterial);
                    scene.add(axis);

                    // 線分の接点に丸を追加
                    const sphereGeometry: THREE.SphereGeometry = new THREE.SphereGeometry(radian);
                    const joint0 = new THREE.Mesh(sphereGeometry, trajMaterial);
                    joint0.position.set(head.x, head.y, head.z)
                    scene.add(joint0);
                    const joint1 = new THREE.Mesh(sphereGeometry, trajMaterial);
                    joint1.position.set(tail.x, tail.y, tail.z)
                    scene.add(joint1);
                }
            });
    });

    // MapBoxの座標(緯度経度)をWebGLの3D空間上の座標に変換する関数
    function getLoaction(d: MyData, fristDate: number, lastDate: number): THREE.Vector3 {
        const lngLat = new mapboxgl.LngLat(d.lng, d.lat);
        const dist = map.getCenter().distanceTo(lngLat);
        const bearing = getBearing(map.getCenter(), lngLat);

        const x: number = dist * Math.sin((Math.PI * 2 * bearing) / 360);
        const z: number = dist * Math.cos((Math.PI * 2 * bearing) / 360);

        const date: number = new Date(d.date).valueOf();
        const ratio: number = (date - fristDate) / (lastDate - fristDate);
        const y: number = 400000 * ratio;

        return new THREE.Vector3(x, y, -z);
    }

    function getBearing(pos1: mapboxgl.LngLat, pos2: mapboxgl.LngLat): number {
        const radian = (deg: number): number => deg * Math.PI / 180;
        const degree = (rad: number): number => rad * 180 / Math.PI;

        let dLng = pos2.lng - pos1.lng;

        const x = Math.cos(radian(pos2.lat)) * Math.sin(radian(dLng));
        const z = Math.cos(radian(pos1.lat)) * Math.sin(radian(pos2.lat)) - Math.sin(radian(pos1.lat)) * Math.cos(radian(pos2.lat)) * Math.cos(radian(dLng));

        let bearingRadian = Math.atan2(x, z);
        let bearingDegree = degree(bearingRadian);

        if (bearingDegree > 0.0) return bearingDegree;
        return bearingDegree + 360;
    }
});
```