【数字孪生平台】使用 Three.js 绘制月球陨石坑

在本文中，我们将使用 NASA 发布的月球地形数据“SLDEM2015”在 Three.js 上绘制月球陨石坑。更多精彩内容尽在数字孪生平台。

关于SLDEM2015

“SLDEM2015”是详细描绘月球地形的数字高程模型（DEM）。这些数据由 NASA 的月球勘测轨道器 (LRO)收集，并提供了月球表面高度的高分辨率地图。

这些数据可以在 QGIS 中查看，就像任何普通的数字高程模型 (DEM) 一样。　

这个数据的CRS是“GCS_Moon_2000”，也就是“2000年月球地心坐标系”，是参考月球地理位置和特征的坐标系。与地球坐标系一样，月球坐标系也以纬度和经度表示，但使用月球自己的参考点和平面。

使用 Three.js 创建外太空

首先，创建场景并创建外太空。这里我使用 NASA 发布的 EXR 数据作为场景背景。

使用EXRLoader加载.exr文件并将其设置为场景的背景。可见threejs官方示例。

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 的数据，该数据在此可下载。下载SLDEM2015_512_00N_30N_315_360.JP2 文件。文件名描述了数据的分辨率和范围，分辨率为从赤道到北纬30度、西经45度到西经0度（或东经315度到东经360度）每度512像素。

用QGIS显示下载的SLDEM2015，该图像中心附近的大陨石坑就是哥白尼陨石坑。

使用 gdal 工具创建仅剪切哥白尼部分的栅格：

现在数据已经准备好了，接下来我们用 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像素/度。月球赤道的周长约为10917公里。将其除以 360 度，每度大约 30.32 公里。因此，一个像素所代表的距离可以计算如下：
30.32 km/度 ÷ 512 像素/度 = 0.05921875 km
因此，该数据中一个像素代表的实际距离约为59.22米。

我们可以将 PlaneGeometry 的 1 个网格大小乘以 59.22，但是对于 Three.js 上的坐标单位来说，它会相当大，所以高度和宽度我们以 1/10 比例来绘制。创建PlaneGeometry时，将栅格宽度乘以5.922分配给宽度和高度，并将每个顶点（段）的高度分配为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));

    // PlaneGeometry默认为纵向，因此绕 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; // 将每个顶点沿 y 方向移动 SLDEM 的值
        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 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); // 使用计算的颜色设置片元颜色
    }
`;