Three.jsで月面クレーターを再現する

これは MIERUNE AdventCalendar 2023 11日目の記事です! 昨日は@soramiさんによる「アジアでオープンなジオの祭典: 「FOSS4G ASIA 2023」に参加&発表してきました！」でした。

はじめに

この記事では、NASAが公開している月面地形データ「SLDEM2015」を使用して、Three.js上に月面クレーターを表示させます。

DEMO
https://satoshi7190.github.io/threejs-lunar-crater-demo/

SLDEM2015 について

「SLDEM2015」は、月の地形を詳細に描写したデジタル標高モデル（DEM）です。このデータは、NASAのLunar Reconnaissance Orbiter(LRO）によって収集されたもので、月の地表の高さを高い解像度でマッピングしたものです。

このデータは、通常の数値標高モデル（DEM）と同じようにQGISで確認できます。　

ちなみに、このデータのCRSは「GCS_Moon_2000」というものになります。
これは「Geocentric Coordinate System for the Moon 2000」という意味で、月の地理的な位置や特徴を参照するための座標系です。地球の座標系と同様、月の座標系も緯度と経度で表されますが、月固有の基準点と基準面を使用します。

データを使用する際は下記の引用文を記載する必要があります

Barker, M. K., Mazarico, E., Neumann, G. A., Zuber, M. T., Haruyama, J., Smith, D. E. "A new lunar digital elevation model from the Lunar Orbiter Laser Altimeter and SELENE Terrain Camera," Icarus, Volume 273, p. 346-355. http://dx.doi.org/10.1016/j.icarus.2015.07.039

Three.jsで宇宙空間を作る

まずはシーンを作成し、宇宙空間を作ります。せっかくなのでシーンの背景にもNASAが公開しているこちらのexrデータを使用します。

こちらのデータを使用するときは下記のクレジット表記が必要とのことです。

NASA's Scientific Visualization Studio

EXRLoaderを使用して.exrファイルを読み込んで、シーンの背景にセットします。

例
https://threejs.org/examples/#webgl_materials_envmaps_exr

コード

import * as THREE from 'three';
import { MapControls } from 'three/examples/jsm/controls/MapControls.js';
import { EXRLoader } from 'three/examples/jsm/loaders/EXRLoader.js';

// 画面サイズ
const sizes = {
    width: window.innerWidth,
    height: window.innerHeight,
};

// キャンバスを作成
const canvas = document.createElement('canvas');
document.body.appendChild(canvas);

// シーンを作成
const scene = new THREE.Scene();

// .exrファイルを読み込み、背景に設定
new EXRLoader().load('./starmap_random_2020_4k.exr', (texture) => {
    texture.mapping = THREE.EquirectangularReflectionMapping;
    scene.background = texture;
});

// カメラ
const camera = new THREE.PerspectiveCamera(75, sizes.width / sizes.height, 0.1, 100000);
camera.position.set(900, 1100, 1800);

// コントロール
const controls = new MapControls(camera, canvas);
controls.enableDamping = true;
controls.maxDistance = 4000;

// レンダー
const renderer = new THREE.WebGLRenderer({
    canvas: canvas,
    alpha: true,
});
renderer.setSize(sizes.width, sizes.height);
renderer.setPixelRatio(Math.min(window.devicePixelRatio, 2));

// レンダリング
const animate = () => {
    requestAnimationFrame(animate);
    controls.update();
    renderer.render(scene, camera);
};
animate();

// 画面リサイズ時にキャンバスもリサイズ
window.addEventListener(
    'resize',
    () => {
        camera.aspect = window.innerWidth / window.innerHeight;
        camera.updateProjectionMatrix();
        renderer.setSize(window.innerWidth, window.innerHeight);
    },
    false,
);

宇宙空間が描画されました。

月面クレーターを表示する

次に、SLDEM2015を使用して、シーン内に月面クレーター（地形）を表示します。今回は、比較的大きなクレーターで、地球から見てもその存在が確認できるコペルニクスを表示することにします。

コペルニクスのおよその位置座標は以下になります。
緯度: 約9.7°N
経度: 約20.0°W

SLDEM2015のより改造の高いタイル分割された方のデータを使用します。

下記のファイル名の.JP2をダウンロードします。

• SLDEM2015_512_00N_30N_315_360

ファイル名には、データの解像度とその範囲が記載されており、このデータは512ピクセル毎度解像度で、赤道から30度北緯、および西経45度から西経0度（または東経315度から東経360度）の範囲のものであることがよみとれます。

QGISでダウンロードしたSLDEM2015を表示します。このラスター画像の中心近くにあるく大きなクレーターがコペルニクスです。

gdalツールを使用してコペルニクス部分のみを切り抜いたラスターを作成します。

データの用意ができたのでThree.js上に載せたいと思います。ラスターデータを読み込むためにgeotiff.jsを使用し、ラスターを読み込む関数を記述します。この関数でラスタの各ピクセルの値、ラスターの縦幅、横幅を返すようにします。

npm install geotiff

import { fromArrayBuffer } from 'geotiff';

// ラスターデータの読み込み
const loadRasterData = async (data) => {
    const response = await fetch(data);
    const arrayBuffer = await response.arrayBuffer();
    const tiff = await fromArrayBuffer(arrayBuffer);
    const image = await tiff.getImage();

    // ラスターデータを取得
    const rasters = await image.readRasters();
    const raster = rasters[0];

    // ラスターの高さと幅を取得
    const tiffWidth = image.getWidth();
    const tiffHeight = image.getHeight();

    return { raster, tiffWidth, tiffHeight };
};

そして、この関数を使用して、シーンにPlaneGeometryを作成し、ラスターの値を各メッシュの頂点に割り当てることでクレーターの地形をシーン内に再現します。今回はPlaneGeometryの頂点数が多くなるため、頂点を操作する処理は、シェーダーマテリアルでシェーダーを書き、GPU処理で行うことにします。

このときに、クレーターの高さ（深さ）のスケールを、作成するPlaneGeometryの幅と一致させる必要があります。各ピクセルの高さ情報はメートル単位になっているので、これに合わせていきます。

今回使用しているSLDEM2015は解像度が512ピクセル/度のものを使用しています。月の周囲距離を考慮して、この値を度単位からメートル単位に再計算します。
月の赤道周囲は約10,917キロメートルです。これを360度で割ると、1度あたり約30.32キロメートルになります。したがって、1ピクセルが表す距離は以下のように計算できます。

30.32km/度 ÷ 512ピクセル/度 = 0.05921875km

よって、このデータにおける1ピクセルが表す実際の距離は約59.22メートルです。
単純な計算をすれば、PlaneGeometryの１メッシュに対して59.22を乗算すればよいのですが、Three.js上のワールの座標の単位の場合はかなりの大きさになってしまうため、高さと幅を10分の１スケールで描画することにします。
PlaneGeometryを作成するときに、ラスターの幅に対して5.922を乗算したものをwidth、heightに割り当て、各頂点（segment）の高さはSLDEMの値に0.1を乗算したものを割り当てることにしました。

// uniform変数
const uniforms = {
    uMax: { value: 0 },
    uMin: { value: 0 },
};

// 地形の作成
const int = async (sldem) => {
    // ラスターデータの読み込み
    const sldemData = await loadRasterData(sldem);

    // マテリアルの作成
    const planeMaterial = new THREE.ShaderMaterial({
        vertexShader: vertexShader,
        fragmentShader: fragmentShader,
        uniforms: uniforms,
        transparent: true,
        side: THREE.DoubleSide,
    });

    // SLDEMの値の最小値を計算
    const minValue = sldemData.raster.reduce((min, value) => {
        return Math.min(min, value * 0.1);
    }, Infinity);

    // SLDEMの値の最大値を計算
    const maxValue = sldemData.raster.reduce((max, value) => Math.max(max, value * 0.1), -Infinity);

    // uniform変数にSLDEMの最大値と最小値を設定
    planeMaterial.uniforms.uMax.value = maxValue;
    planeMaterial.uniforms.uMin.value = minValue;

    // PlaneGeometryの作成。ラスターのピクセル数分の頂点を持つ
    const planegeometry = new THREE.PlaneGeometry(sldemData.tiffWidth * 5.922, sldemData.tiffHeight * 5.922, sldemData.tiffWidth - 1, sldemData.tiffHeight - 1);
    const demValues = new Float32Array(sldemData.raster.map((value) => value * 0.1));
    
    // 各頂点にSLDEMの値を渡す 
    planegeometry.setAttribute('demValues', new THREE.BufferAttribute(demValues, 1));

    // デフォルトは縦向きなのでをX軸を中心に90度に回転
    planegeometry.rotateX(-Math.PI / 2);
    planegeometry.computeVertexNormals();

    // メッシュの作成
    const planegmesh = new THREE.Mesh(planegeometry, planeMaterial);
    scene.add(planegmesh);
};

//ラスターの読み込み処理を実行
int('./SLDEM2015.tif').catch((error) => {
    console.error('Error loading raster data:', error);
});

頂点シェーダーとフラグメントシェーダーも書いていきます。

// 頂点シェーダー
const vertexShader = /* glsl */ `
    attribute float demValues;
    varying vec3 fragNormal;
    varying vec3 fragPosition;
    varying vec3 vPosition;

    void main() {
        vec3 pos = position;
        pos.y += demValues; // 各頂点をSLDEMの値分y方向に動かす
        vPosition = pos;
        gl_Position = projectionMatrix * modelViewMatrix * vec4(pos, 1.0);
    }
`;

// フラグメントシェーダー
const fragmentShader = /* glsl */ `
    varying vec3 fragPosition;
    varying vec3 vPosition;
    uniform float uMax; // SLDEMの最大値
    uniform float uMin; // SLDEMの最小値

    //法線ベクトルを求める関数
    vec3 getNormal ( vec3 position ) {
        vec3 dx = dFdx( position );
        vec3 dy = dFdy( position );
        return normalize( cross(dx, dy) );
    }

    void main() {
        // 陰影表示
        vec3 normal = getNormal(vPosition); // 法線ベクトルを取得
        vec3 lightPosition = vec3(500.0, 500.0, -1000.0); // 光源の位置
    
        // 光源からの方向ベクトルを計算
        vec3 lightDir = normalize(lightPosition - fragPosition); 
    
        // 拡散反射の計算
        float diff = max(dot(normal, lightDir), 0.0);
    
        // 拡散反射に基づく色の計算
        vec3 diffuseColor = vec3(1.0, 1.0, 1.0) * diff; 
    
        // 最終的な色をフラグメントの色として設定
        gl_FragColor = vec4(diffuseColor, 1.0); 
    }
`;

シーンにクレーターの地形が表示されました。

ちなみに、フラグメントシェーダーを書き換えることで、陰影表示以外の表現もできます。

標高値による色分け

const fragmentShader = /* glsl */ `
    varying vec3 fragPosition;
    varying vec3 vPosition;
    uniform float uMax; // SLDEMの最大値
    uniform float uMin; // SLDEMの最小値
    
    //法線ベクトルを求める関数
    vec3 getNormal ( vec3 position ) {
        vec3 dx = dFdx( position );
        vec3 dy = dFdy( position );
        return normalize( cross(dx, dy) );
    }
    
    void main() {
        // 色の設定
        vec3 lowColor = vec3(0.129,0.067,0.306); // 低い標高の色（濃紫）
        vec3 midColor = vec3(0.729,0.22,0.471);  // 中間の標高の色（ピンク）
        vec3 highColor = vec3(0.996,0.631,0.431); // 高い標高の色（オレンジ）
        
        // 標高を正規化
        float normalizedHeight = (vPosition.y - uMin) / (uMax - uMin);
        float intensity = clamp(normalizedHeight, 0.0, 1.0); // 標高の強度を0.0から1.0の間に制限
        
        vec3 color;
        // 標高に基づいて色を補間
        if (intensity < 0.5) {
            // 標高が中間より低い場合、低い標高の色と中間の標高の色の間で補間
            color = mix(lowColor, midColor, intensity * 2.0);
        } else {
            // 標高が中間より高い場合、中間の標高の色と高い標高の色の間で補間
            color = mix(midColor, highColor, (intensity - 0.5) * 2.0);
        }
        
        gl_FragColor = vec4(color, 1.0); // 計算された色でフラグメントの色を設定
    }
`;

傾斜量による色分け

const fragmentShader = /* glsl */ `
    varying vec3 fragPosition;
    varying vec3 vPosition;
    uniform float uMax; // SLDEMの最大値
    uniform float uMin; // SLDEMの最小値
    
    //法線ベクトルを求める関数
    vec3 getNormal ( vec3 position ) {
        vec3 dx = dFdx( position );
        vec3 dy = dFdy( position );
        return normalize( cross(dx, dy) );
    }
    
    void main() {
        // 色の設定
        vec3 highColor = vec3(0.847,0.949,0.878); // 傾斜が高い部分の色（浅い緑）
        vec3 midColor = vec3(0.522,0.851,0.694);  // 中間の傾斜の色（中間の緑）
        vec3 lowColor = vec3(0.204,0.518,0.647);  // 傾斜が低い部分の色（濃い青）
        
        vec3 normal = getNormal(vPosition); // 法線ベクトルを取得
        float intensity = abs(normal.y); // 傾斜量を計算（Y成分の絶対値を使用）
        
        vec3 color;
        // 傾斜量に基づいて色を補間
        if (intensity < 0.5) {
            // 傾斜が少ない場合は、中間色と低傾斜色の間で補間
            color = mix(highColor, midColor, intensity * 2.0);
        } else {
            // 傾斜が多い場合は、高傾斜色と中間色の間で補間
            color = mix(midColor, lowColor, (intensity - 0.5) * 2.0);
        }
        
        gl_FragColor = vec4(color, 1.0); // 計算された色でフラグメントの色を設定
    }
`;

おわりに

このやり方は通常の数値標高モデル（DEM）も可能ですが、単にDEMを表示するのは面白くないと思ったので、地球上の地形データではなく、月の地形データを使ってみました。

明日は@chizutodesignさんによる記事です！お楽しみに！！