Add topography to mesh

app/components/Form.tsx

@@ -166,6 +166,7 @@ export function Form() {
     if (!sceneView) return;
 
     sceneView.toggleLayerVisibility(name);
+    sceneView.dispatchChangeEvent();
   }
 
   function handleChangeTopography() {
@@ -259,9 +260,12 @@ export function Form() {
             <div className="flex flex-col gap-2">
               {sceneView?.model.children.map((child) => {
                 const key = `toggle-visibility-${child.name}`;
-                const color = `#${(
+                let color = "transparent";
-                  (child as Mesh).material as MeshStandardMaterial
+                if ((child as Mesh).material instanceof MeshStandardMaterial) {
-                ).color.getHexString()}`;
+                  color = `#${(
+                    (child as Mesh).material as MeshStandardMaterial
+                  ).color.getHexString()}`;
+                }
                 const visible = (child as Mesh).visible;
 
                 return (

app/three/SceneView.ts

@@ -1,10 +1,13 @@
 import {
+  Frustum,
   Group,
   Material,
+  Matrix4,
   Mesh,
   MeshBasicMaterial,
   MeshPhongMaterial,
   MeshStandardMaterial,
+  Object3D,
   PerspectiveCamera,
   Plane,
   Raycaster,
@@ -16,7 +19,12 @@ import {
 } from "three";
 import { buildMeshes } from "./utils/build-meshes";
 import { Extent, buildScene } from "./utils/build-scene";
-import { getMetadata, transform } from "./utils/utils";
+import {
+  getFrustumIntersections,
+  getMetadata,
+  tileBounds,
+  transform,
+} from "./utils/utils";
 import { MODEL_ID, SERVICE_URL } from "./config";
 import {
   Orientation,
@@ -32,11 +40,11 @@ import {
   OrbitControls,
 } from "three/examples/jsm/Addons.js";
 import {
+  LODFrustum,
   LODRaycast,
   MapPlaneNode,
   MapView,
   OpenStreetMapsProvider,
-  UnitsUtils,
 } from "geo-three";
 import { Data, createSVG } from "./utils/create-borehole-svg";
 import { TileData, updateTiles } from "./ShaderMaterial";
@@ -435,6 +443,10 @@ export class SceneView extends EventTarget {
       this._resetClippingBox();
     }
   }
+
+  dispatchChangeEvent() {
+    this._orbitControls.dispatchEvent({ type: "change" });
+  }
 }
 
 async function init(container: HTMLElement, modelId = MODEL_ID) {
@@ -483,20 +495,12 @@ async function init(container: HTMLElement, modelId = MODEL_ID) {
 
   // Create a map tiles provider object
   const provider = new OpenStreetMapsProvider();
-  //const heightProvider = new MapTilerProvider(
-  //  MAPTILER_API_KEY,
-  //  "tiles",
-  //  "terrain-rgb",
-  //  "png"
-  //);
 
   // Create the map view for OSM topography
-  const lod = new LODRaycast();
+  const lod = new LODFrustum();
 
   const map = new MapView(MapView.PLANAR, provider);
   map.lod = lod;
-  // const customNode = new CustomMapHeightNodeShader(undefined, map);
-  // map.setRoot(customNode);
   map.rotateX(Math.PI / 2);
 
   map.name = "topography";
@@ -505,39 +509,15 @@ async function init(container: HTMLElement, modelId = MODEL_ID) {
 
   controls.addEventListener("change", () => {
     const tiles: TileData[] = [];
-    function traverse(node: MapPlaneNode) {
+    camera.updateMatrix();
-      if (node.isMesh) {
-        const bounds = UnitsUtils.tileBounds(node.level, node.x, node.y);
 
-        const xmin = bounds[0];
+    const frustumPoints = getFrustumIntersections(camera);
-        const ymin = bounds[2];
-        const xmax = xmin + bounds[1];
-        const ymax = ymin + bounds[3];
-
-        if (
-          (extent.xmax >= xmin && extent.ymax >= ymin) ||
-          (extent.xmin <= xmax && extent.ymax >= ymin) ||
-          (extent.xmin <= xmax && extent.ymin <= ymax) ||
-          (extent.xmax >= xmin && extent.ymin <= ymax)
-        ) {
-          tiles.push({
-            xmin,
-            ymin,
-            xmax,
-            ymax,
-            texture: (node.material as MeshPhongMaterial).map,
-          });
-        }
-      }
-      for (const c of node.children) {
-        traverse(c as MapPlaneNode);
-      }
-    }
 
     map.lod.updateLOD(map, camera, renderer, scene);
-    traverse(map.root);
+    traverse(map.root, extent, tiles, frustumPoints);
+    tiles.sort((a, b) => b.zoom - a.zoom);
 
-    updateTiles(tiles.reverse());
+    updateTiles(tiles);
   });
 
   return {
@@ -549,3 +529,48 @@ async function init(container: HTMLElement, modelId = MODEL_ID) {
     renderer,
   };
 }
+
+function traverse(
+  node: MapPlaneNode,
+  extent: Extent,
+  tiles: TileData[],
+  frustumPoints: Vector3[]
+) {
+  const bounds = tileBounds(node.level, node.x, node.y);
+
+  const xmin = bounds[0];
+  const ymin = bounds[2];
+  const xmax = xmin + bounds[1];
+  const ymax = ymin + bounds[3];
+
+  if (
+    ((xmax >= extent.xmin && xmax <= extent.xmax) ||
+      (xmin >= extent.xmin && xmin <= extent.xmax)) &&
+    ((ymax >= extent.ymin && ymax <= extent.ymax) ||
+      (ymin >= extent.ymin && ymin <= extent.ymax))
+  ) {
+    if (
+      frustumPoints.length < 2 ||
+      (((xmax >= frustumPoints[0].x && xmax <= frustumPoints[1].x) ||
+        (xmin >= frustumPoints[0].x && xmin <= frustumPoints[1].x)) &&
+        ((ymax >= frustumPoints[0].y && ymax <= frustumPoints[1].y) ||
+          (ymin >= frustumPoints[0].y && ymin <= frustumPoints[1].y)))
+    ) {
+      const texture = (node.material as MeshPhongMaterial).map;
+      if (texture) {
+        tiles.push({
+          xmin,
+          ymin,
+          xmax,
+          ymax,
+          zoom: node.level,
+          texture,
+        });
+      }
+    }
+  }
+
+  for (const c of node.children) {
+    traverse(c as MapPlaneNode, extent, tiles, frustumPoints);
+  }
+}

app/three/ShaderMaterial.ts

@@ -1,14 +1,23 @@
-import { ShaderMaterial, Texture, Vector4 } from "three";
+import {
+  DataArrayTexture,
+  LinearFilter,
+  RGBAFormat,
+  ShaderChunk,
+  ShaderMaterial,
+  Texture,
+  Vector4,
+} from "three";
 
 export interface TileData {
   xmin: number;
   ymin: number;
   xmax: number;
   ymax: number;
-  texture: Texture | null;
+  zoom: number;
+  texture: Texture;
 }
 
-const maxTiles = 16;
+const maxTiles = 24;
 
 // Initialize empty texture slots
 const dummyTexture = new Texture();
@@ -18,25 +27,39 @@ dummyTexture.needsUpdate = true;
 // Create shader material
 export const shaderMaterial = new ShaderMaterial({
   uniforms: {
-    tileTextures: { value: Array(maxTiles).fill(dummyTexture) },
     tileBounds: { value: Array(maxTiles).fill(new Vector4(0, 0, 0, 0)) },
     tileCount: { value: 0 },
+    tiles: { value: null },
   },
-  vertexShader: `
+  vertexShader:
-        varying vec3 vWorldPosition;
+    ShaderChunk.common +
-        void main() {
+    "\n" +
-            vWorldPosition = (modelMatrix * vec4(position, 1.0)).xyz;
+    ShaderChunk.logdepthbuf_pars_vertex +
-            gl_Position = projectionMatrix * viewMatrix * vec4(vWorldPosition, 1.0);
+    `
-        }
-    `,
-  fragmentShader: `
-        uniform sampler2D tileTextures[${maxTiles}];
-        uniform vec4 tileBounds[${maxTiles}];
-        uniform int tileCount;
         varying vec3 vWorldPosition;
+        varying float fragDepth;
 
         void main() {
-            vec4 color = vec4(1.0, 1.0, 1.0, 1.0); // Default color
+            vWorldPosition = (modelMatrix * vec4(position, 1.0)).xyz;
+            gl_Position = projectionMatrix * modelViewMatrix * vec4(position, 1.0);
+            fragDepth = (gl_Position.z / gl_Position.w + 1.0) * 0.5;
+        
+    ` +
+    ShaderChunk.logdepthbuf_vertex +
+    `
+  }
+`,
+  fragmentShader:
+    ShaderChunk.logdepthbuf_pars_fragment +
+    `
+        uniform vec4 tileBounds[${maxTiles}];
+        uniform int tileCount;
+        uniform sampler2DArray tiles;
+        varying vec3 vWorldPosition;
+        varying float fragDepth;
+
+        void main() {
+            vec4 color = vec4(191.0/255.0, 209.0/255.0, 229.0/255.0, 1.0); // Default color
 
             for (int i = 0; i < ${maxTiles}; i++) {
                 if (i >= tileCount) break; // Only process available tiles
@@ -47,32 +70,20 @@ export const shaderMaterial = new ShaderMaterial({
                     vWorldPosition.y >= bounds.z && vWorldPosition.y <= bounds.w) {
                     
                     vec2 uv = (vWorldPosition.xy - bounds.xz) / (bounds.yw - bounds.xz);
-                    switch (i) {
+                    uv = vec2(uv.x, 1.0 - uv.y);
-                      case 0: color = texture2D(tileTextures[0], uv); break;
+                    color = texture2D(tiles, vec3(uv, i));
-                      case 1: color = texture2D(tileTextures[1], uv); break;
-                      case 2: color = texture2D(tileTextures[2], uv); break;
-                      case 3: color = texture2D(tileTextures[3], uv); break;
-                      case 4: color = texture2D(tileTextures[4], uv); break;
-                      case 5: color = texture2D(tileTextures[5], uv); break;
-                      case 6: color = texture2D(tileTextures[6], uv); break;
-                      case 7: color = texture2D(tileTextures[7], uv); break;
-                      case 8: color = texture2D(tileTextures[8], uv); break;
-                      case 9: color = texture2D(tileTextures[9], uv); break;
-                      case 10: color = texture2D(tileTextures[10], uv); break;
-                      case 11: color = texture2D(tileTextures[11], uv); break;
-                      case 12: color = texture2D(tileTextures[12], uv); break;
-                      case 13: color = texture2D(tileTextures[13], uv); break;
-                      case 14: color = texture2D(tileTextures[14], uv); break;
-                      case 15: color = texture2D(tileTextures[15], uv); break;
-                    }
 
-                  break; // Stop checking once we find the correct tile
+                    break; // Stop checking once we find the correct tile
                 }
               }
 
             gl_FragColor = color;
-        }
+            gl_FragDepth = fragDepth;
-    `,
+    ` +
+    ShaderChunk.logdepthbuf_fragment +
+    `
+  }
+`,
 });
 
 export function updateTiles(newTiles: TileData[]) {
@@ -92,7 +103,53 @@ export function updateTiles(newTiles: TileData[]) {
   }
 
   // Update shader uniforms
-  shaderMaterial.uniforms.tileTextures.value = textures;
   shaderMaterial.uniforms.tileBounds.value = bounds;
   shaderMaterial.uniforms.tileCount.value = newTiles.length;
+  shaderMaterial.uniforms.tiles.value = createDataArrayTexture(textures);
+}
+
+// Create a buffer with color data
+const width = 256;
+const height = 256;
+const size = width * height;
+function createDataArrayTexture(textures: Texture[]) {
+  const depth = textures.length;
+
+  const data = new Uint8Array(4 * size * depth);
+
+  for (let i = 0; i < depth; i++) {
+    const texture = textures[i];
+    const imageData = getImageData(texture);
+
+    if (imageData) {
+      data.set(imageData, i * size * 4);
+    }
+  }
+
+  // Use the buffer to create a DataArrayTexture
+  const texture = new DataArrayTexture(data, width, height, depth);
+  texture.format = RGBAFormat;
+  texture.generateMipmaps = false;
+  texture.magFilter = LinearFilter;
+  texture.minFilter = LinearFilter;
+  texture.needsUpdate = true;
+  return texture;
+}
+
+// Create a canvas and draw the image on it
+const canvas = new OffscreenCanvas(width, height);
+const ctx = canvas.getContext("2d");
+function getImageData(texture: Texture) {
+  const image = texture.source.data;
+
+  // Draw the image onto the canvas
+  if (ctx) {
+    ctx.drawImage(image, 0, 0);
+
+    // Get the pixel data from the canvas
+    const imageData = ctx.getImageData(0, 0, canvas.width, canvas.height);
+    return imageData.data;
+  } else {
+    return null;
+  }
 }

app/three/config.ts

@@ -5,5 +5,3 @@ export const SERVICE_URL =
 export const VERTICES_URL = "https://geusegdi01.geus.dk/geom3d/data/nodes/";
 export const TRIANGLE_INDICES_URL =
   "https://geusegdi01.geus.dk/geom3d/data/triangles/";
-
-export const MAPTILER_API_KEY = "1JkD1W8u5UM5Tjd8r3Wl ";

app/three/utils/build-meshes.ts

@@ -8,6 +8,7 @@ import {
 
 import { fetchVertices, fetchTriangleIndices, transform } from "./utils";
 import { TRIANGLE_INDICES_URL, VERTICES_URL } from "../config";
+import { shaderMaterial } from "../ShaderMaterial";
 
 interface MappedFeature {
   featuregeom_id: number;
@@ -60,14 +61,17 @@ async function buildMesh(layerData: MappedFeature) {
 
   const material = new MeshStandardMaterial({
     color: color,
-    metalness: 0.0,
+    metalness: 0.1,
-    roughness: 5.0,
+    roughness: 0.5,
     flatShading: true,
     side: DoubleSide,
     wireframe: false,
   });
 
-  const mesh = new Mesh(geometry, material);
+  const mesh = new Mesh(
+    geometry,
+    name === "Topography" ? shaderMaterial : material
+  );
   mesh.name = name;
   mesh.userData.layerId = geomId;
   mesh.castShadow = true;

app/three/utils/utils.ts

@@ -1,4 +1,4 @@
-import { Vector3 } from "three";
+import { PerspectiveCamera, Plane, Raycaster, Vector2, Vector3 } from "three";
 import { Extent } from "./build-scene";
 import { unpackEdges, unpackVertices } from "./decoders";
 import proj4 from "proj4";
@@ -60,12 +60,69 @@ export async function fetchVertices(pointUrl: string, geomId: string) {
   return unpackVertices(buffer);
 }
 
-// Transformation from EPSG 3034 to EPSG 3857
+// Transformation from EPSG 3034 to a modified EPSG 900913 using the mean radius to align with geo-three
 const SOURCE = "EPSG:3034";
-const PROJ_STRING =
+const DEST = "EPSG:900913";
-  "+proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs";
+proj4.defs(
-const DEST = "EPSG:3857";
+  SOURCE,
-proj4.defs(SOURCE, PROJ_STRING);
+  "+proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs"
+);
+proj4.defs(
+  DEST,
+  "+proj=merc +a=6371008 +b=6371008 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs"
+);
 export function transform(p: number[]) {
   return proj4(SOURCE, DEST, p);
 }
+
+const MAX_EXTENT = 20015111.9287618;
+function getTileSize(zoom: number): number {
+  const numTiles = Math.pow(2, zoom);
+  return (2 * MAX_EXTENT) / numTiles;
+}
+
+export function tileBounds(zoom: number, x: number, y: number): number[] {
+  const tileSize = getTileSize(zoom);
+  const minX = -MAX_EXTENT + x * tileSize;
+  const minY = MAX_EXTENT - (y + 1) * tileSize;
+  return [minX, tileSize, minY, tileSize];
+}
+
+const plane = new Plane(new Vector3(0, 0, 1), 0);
+
+const corners = [new Vector2(-1, -1), new Vector2(1, 1)];
+
+export const getFrustumIntersections = (camera: PerspectiveCamera) => {
+  const points = [];
+
+  const raycaster = new Raycaster();
+
+  for (const ndc of corners) {
+    // Convert NDC to world space ray
+    raycaster.setFromCamera(ndc, camera);
+
+    // Find intersection with the XY plane
+    const intersection = new Vector3();
+    const hit = raycaster.ray.intersectPlane(plane, intersection);
+    if (hit) {
+      points.push(intersection.clone());
+    }
+  }
+
+  if (points.length > 1) {
+    return [
+      new Vector3(
+        Math.min(points[0].x, points[1].x),
+        Math.min(points[0].y, points[1].y),
+        0
+      ),
+      new Vector3(
+        Math.max(points[0].x, points[1].x),
+        Math.max(points[0].y, points[1].y),
+        0
+      ),
+    ];
+  }
+
+  return points;
+};