Merge pull request #222 from open-pv/feature/geotiff-terrain

Feature/geotiff terrain

package.json

@@ -10,6 +10,7 @@
     "@emotion/styled": "^11.13.0",
     "@openpv/simshady": "^0.0.3",
     "@react-three/drei": "^9.111.1",
+    "geotiff": "^2.1.3",
     "i18next": "^23.2.11",
     "i18next-http-backend": "^2.2.1",
     "maplibre-gl": "^4.5.2",

src/components/ThreeViewer/Footer.js

@@ -22,6 +22,21 @@ function Footer({ federalState, frontendState }) {
               dl-de/by-2-0
             </a>
             )
+
+          |
+
+            Geländemodell:&nbsp;
+            <a href="https://sonny.4lima.de" target="_blank">
+              &copy;&nbsp;Sonny
+            </a>
+            &nbsp;(
+            <a href="https://creativecommons.org/licenses/by/4.0/deed.en" target="_blank">
+              CC-BY-4.0
+            </a>
+            ), erstellt aus 
+            <a href="https://drive.google.com/file/d/1rgGA22Ha42ulQORK9Pfp4JPpPAIKFx6Q/view" target="_blank">
+              verschiedenen Quellen 
+            </a>
           </p>
         )}
         {federalState && (

src/components/ThreeViewer/Terrain.js

@@ -1,26 +1,69 @@
 import React, { useEffect, useState } from "react"
 import * as THREE from "three"
-import { tile2meters } from "../../simulation/download"
-import { coordinatesXY15 } from "../../simulation/location"
+import { tile2meters, mercator2meters } from "../../simulation/download"
+import { coordinatesWebMercator, coordinatesXY15, xyzBounds } from "../../simulation/location"
+import * as geotiff from 'geotiff';
+
+class ElevationManager {
+  static instancePromise = null;
+  static state = 'uninitialized';
+
+  tiff = null;
+  image = null;
+  width = -1;
+  height = -1;
+  boundingBox = [0, 0, 0, 0];
+
+  static async toPoint3D(x, y) {
+    const me = await this.getInstance(); 
+    const px = (x - me.tiepoint[3]) / me.pixelScale[0]
+    const py = (y - me.tiepoint[4]) / -me.pixelScale[1]
+
+    // Retrieve pixel data
+    const rasterData = await me.image.readRasters({ window: [
+      Math.floor(px), Math.floor(py),
+      Math.ceil(px)+1, Math.ceil(py)+1] });
+
+    // bilinear interpolation
+    const tx = px % 1;
+    const ty = py % 1;
+    const qx = 1 - tx;
+    const qy = 1 - ty;
+    const z = qx * qy * rasterData[0][0] +
+              qx * ty * rasterData[0][2] +
+              tx * qy * rasterData[0][1] +
+              tx * ty * rasterData[0][3];
+    return [x, y, z];
+  }
+
+  static async getInstance() {
+    if (!this.instancePromise) {
+      this.instancePromise = this.init();
+    }
+    return this.instancePromise;
+  }
+
+  static async init() {
+    this.instance = new ElevationManager();
+    let me = this.instance
+    me.tiff = await geotiff.fromUrl('https://maps.heidler.info/sonny_dtm_20.tif');
+    me.image = await this.instance.tiff.getImage();
+    me.pixelScale = me.image.fileDirectory.ModelPixelScale;
+    me.tiepoint = me.image.fileDirectory.ModelTiepoint;
+    return me
+  }
+}
+
 
 /** Load an OSM map tile and return it as a THREE Mesh
  */
 const TerrainTile = (props) => {
   const zoom = props.zoom
   const tx = props.x
   const ty = props.y
+  const divisions = props.divisions;
 
   const url = `https://sgx.geodatenzentrum.de/wmts_basemapde/tile/1.0.0/de_basemapde_web_raster_farbe/default/GLOBAL_WEBMERCATOR/${zoom}/${ty}/${tx}.png`
-  const shift = zoom - 12
-  const scale = 1 << (zoom - 15)
-  if (zoom < 12) {
-    console.error("DEM is broken for zoom < 12!")
-  }
-
-  // const dem_url = `https://maps.heidler.info/dem-tiles-12/12/${tx >> shift}/${ty >> shift}.png`;
-  const dem_url = `https://web3d.basemap.de/maplibre/dgm5-rgb/12/${
-    tx >> shift
-  }/${ty >> shift}.png`
 
   let [geometry, setGeometry] = useState(null)
   let [material, setMaterial] = useState(null)
@@ -33,63 +76,48 @@ const TerrainTile = (props) => {
   useEffect(() => {
     async function fetchData() {
       const mapFuture = new THREE.TextureLoader().loadAsync(url)
-      const demFuture = new THREE.TextureLoader().loadAsync(dem_url)
-
-      // DEM Processinkeyg
-      const canvas = document.createElement("canvas")
-      const context = canvas.getContext("2d")
-      const dem = await demFuture
-      canvas.width = dem.image.width
-      canvas.height = dem.image.height
-      context.drawImage(dem.image, 0, 0, canvas.width, canvas.height)
-
-      function sampleDEM(fraction_x, fraction_y) {
-        // Ensure x and y are within bounds
-        if (
-          fraction_x >= 0 &&
-          fraction_x <= 1 &&
-          fraction_y >= 0 &&
-          fraction_y <= 1
-        ) {
-          const x0 = tx - ((tx >> shift) << shift)
-          const y0 = ty - ((ty >> shift) << shift)
-          const s = 1 << shift
-          const x = Math.round(((fraction_x + x0) / s) * (canvas.width - 1))
-          const y = Math.round(((fraction_y + y0) / s) * (canvas.height - 1))
-          // Get image data at the specific (x, y) location
-          const pixelData = context.getImageData(x, y, 1, 1).data
-          const [r, g, b, _] = pixelData
-          const height = -10000 + (r * 256 * 256 + g * 256 + b) * 0.1
-          return height
+
+      // Size of the world map in meters
+      const [x0, y0, x1, y1] = xyzBounds(tx, ty, zoom);
+      let points = [];
+      let uvs = [];
+      let indices = [];
+      let i = 0;
+
+      const row = divisions+1;
+      for (let ty = 0; ty <= divisions; ty++) {
+        for (let tx = 0; tx <= divisions; tx++) {
+          const x = x0 + tx / divisions * (x1 - x0);
+          const y = y0 + ty / divisions * (y1 - y0);
+          points.push(ElevationManager.toPoint3D(x, y));
+          // UV mapping for the texture
+          uvs = uvs.concat([tx / divisions, 1.0 - ty / divisions]);
+          // Triangle indices
+          if(tx > 0 && ty > 0) {
+            indices = indices.concat([
+              i-row-1, i-1, i-row, // 1st triangle
+              i-row, i-1, i // 2nd triangle
+            ]);
+          }
+          i += 1;
         }
       }
 
-      // TODO: Subdivide
-      const corners = [
-        [0, 0],
-        [1, 0],
-        [0, 1],
-        [1, 1],
-      ]
-      const vertices = corners.flatMap(([x, y]) => [
-        // [[tx, ty], [tx+1, ty], [tx, ty+1], [tx+1, ty+1]];
-        tile2meters() * ((tx + x) / scale - coordinatesXY15[0]),
-        -tile2meters() * ((ty + y) / scale - coordinatesXY15[1]),
-        sampleDEM(x, y),
+      const vertices = (await Promise.all(points)).flatMap(([x, y, z]) => [
+        mercator2meters() * (x - coordinatesWebMercator[0]),
+        mercator2meters() * (y - coordinatesWebMercator[1]),
+        z
       ])
-
       const vertexBuffer = new Float32Array(vertices)
-      // UV mapping for the texture
-      const uvs = new Float32Array([0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0])
-      // Triangle indices
-      const indices = new Uint32Array([0, 2, 1, 1, 2, 3])
+      const uvBuffer = new Float32Array(uvs)
+      const indexBuffer = new Uint32Array(indices);
       const geometry = new THREE.BufferGeometry()
       geometry.setAttribute(
         "position",
         new THREE.BufferAttribute(vertexBuffer, 3)
       )
-      geometry.setAttribute("uv", new THREE.BufferAttribute(uvs, 2))
-      geometry.setIndex(new THREE.BufferAttribute(indices, 1))
+      geometry.setAttribute("uv", new THREE.BufferAttribute(uvBuffer, 2))
+      geometry.setIndex(new THREE.BufferAttribute(indexBuffer, 1))
 
       setGeometry(geometry)
       const map = await mapFuture;
@@ -110,28 +138,47 @@ const TerrainTile = (props) => {
 
 const Terrain = ({visible}) => {
   const [x, y] = coordinatesXY15
-  let tiles = []
+  const [tiles, setTiles] = useState([]); // State to manage tiles
   const tx = Math.floor(x * 16)
   const ty = Math.floor(y * 16)
 
   let xys = []
   for (let dx = -11; dx <= 11; dx++) {
     for (let dy = -11; dy <= 11; dy++) {
-      xys.push({ dx: dx, dy: dy })
+      xys.push({ dx, dy, divisions: 2 })
     }
   }
 
   xys.sort((a, b) => a.dx * a.dx + a.dy * a.dy - (b.dx * b.dx + b.dy * b.dy))
-  for (let { dx, dy } of xys) {
-    const key = `${tx + dx}-${ty + dy}-${19}`
-    tiles.push(<TerrainTile key={key} x={tx + dx} y={ty + dy} zoom={19} />)
-  }
+  useEffect(() => {
+    let currentTiles = [];
+
+    // Function to load tiles progressively
+    const loadTiles = (index) => {
+      if (index < xys.length) {
+        const { dx, dy, divisions } = xys[index];
+        const key = `${tx + dx}-${ty + dy}-${19}`;
+        currentTiles.push(<TerrainTile key={key} x={tx + dx} y={ty + dy} divisions={divisions} zoom={19} />);
+
+        setTiles([...currentTiles]); // Update the state with the new set of tiles
+
+        // Schedule the next tile load
+        setTimeout(() => loadTiles(index + 1), 10); // Adjust the timeout for desired loading speed
+      }
+    };
+
+    loadTiles(0); // Start loading tiles
+
+    return () => {
+      setTiles([]); // Clean up on component unmount
+    };
+  }, [tx, ty]); // Dependency array to reset when the coordinates change
 
-  console.log(`Terrain visible:`)
-  console.log(visible);
-  return <group visible={visible}>
-    {tiles}
-  </group>
+  return (
+    <group visible={visible}>
+      {tiles}
+    </group>
+  );
 }
 
 export default Terrain

src/simulation/download.js

@@ -5,13 +5,16 @@ import { GLTFLoader } from "three/addons/loaders/GLTFLoader.js"
 import { attributions } from "../data/dataLicense"
 import {
   coordinatesLonLat,
-  coordinatesXY15,
   projectToWebMercator,
 } from "./location"
 
 export function tile2meters() {
+  return 1222.992452 * mercator2meters();
+}
+
+export function mercator2meters() {
   const lat = coordinatesLonLat[1]
-  return 1222.992452 * Math.cos((lat * Math.PI) / 180.0)
+  return Math.cos((lat * Math.PI) / 180.0)
 }
 
 const dracoLoader = new DRACOLoader()
@@ -116,83 +119,3 @@ async function downloadFile(download_spec) {
     return []
   }
 }
-
-/** Load an OSM map tile and return it as a THREE Mesh
- */
-export async function loadMapTile(tx, ty, zoom) {
-  // const url = `https://tile.openstreetmap.org/${zoom}/${tx}/${ty}.png`;
-  const url = `https://sgx.geodatenzentrum.de/wmts_basemapde/tile/1.0.0/de_basemapde_web_raster_farbe/default/GLOBAL_WEBMERCATOR/${zoom}/${ty}/${tx}.png`
-  const mapFuture = new THREE.TextureLoader().loadAsync(url)
-
-  if (zoom < 12) {
-    console.error("DEM is broken for zoom < 12!")
-  }
-
-  const shift = zoom - 12
-  // const dem_url = `https://maps.heidler.info/dem-tiles-12/12/${tx >> shift}/${ty >> shift}.png`;
-  const dem_url = `https://web3d.basemap.de/maplibre/dgm5-rgb/12/${
-    tx >> shift
-  }/${ty >> shift}.png`
-  const demFuture = new THREE.TextureLoader().loadAsync(dem_url)
-
-  // DEM Processing
-  const canvas = document.createElement("canvas")
-  const context = canvas.getContext("2d")
-  const dem = await demFuture
-  canvas.width = dem.image.width
-  canvas.height = dem.image.height
-  context.drawImage(dem.image, 0, 0, canvas.width, canvas.height)
-
-  function sampleDEM(fraction_x, fraction_y) {
-    // Ensure x and y are within bounds
-    if (
-      fraction_x >= 0 &&
-      fraction_x <= 1 &&
-      fraction_y >= 0 &&
-      fraction_y <= 1
-    ) {
-      const x0 = tx - ((tx >> shift) << shift)
-      const y0 = ty - ((ty >> shift) << shift)
-      const s = 1 << shift
-      const x = Math.round(((fraction_x + x0) / s) * (canvas.width - 1))
-      const y = Math.round(((fraction_y + y0) / s) * (canvas.height - 1))
-      // Get image data at the specific (x, y) location
-      const pixelData = context.getImageData(x, y, 1, 1).data
-      const [r, g, b, _] = pixelData
-      const height = -10000 + (r * 256 * 256 + g * 256 + b) * 0.1
-      return height
-    }
-  }
-
-  const scale = 1 << (zoom - 15)
-  // TODO: Subdivide
-  const corners = [
-    [0, 0],
-    [1, 0],
-    [0, 1],
-    [1, 1],
-  ]
-  const vertices = corners.flatMap(([x, y]) => [
-    // [[tx, ty], [tx+1, ty], [tx, ty+1], [tx+1, ty+1]];
-    tile2meters() * ((tx + x) / scale - coordinatesXY15[0]),
-    -tile2meters() * ((ty + y) / scale - coordinatesXY15[1]),
-    sampleDEM(x, y),
-  ])
-
-  const vertexBuffer = new Float32Array(vertices)
-  // UV mapping for the texture
-  const uvs = new Float32Array([0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0])
-  // Triangle indices
-  const indices = new Uint32Array([0, 2, 1, 1, 2, 3])
-  const geometry = new THREE.BufferGeometry()
-  geometry.setAttribute("position", new THREE.BufferAttribute(vertexBuffer, 3))
-  geometry.setAttribute("uv", new THREE.BufferAttribute(uvs, 2))
-  geometry.setIndex(new THREE.BufferAttribute(indices, 1))
-
-  const material = new THREE.MeshBasicMaterial({
-    map: await mapFuture,
-    side: THREE.DoubleSide,
-  })
-
-  return new THREE.Mesh(geometry, material)
-}