31:T545,| | Description | Cell | Time (ms) | |----|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|----------|---------------| | 1. | The LAZ point cloud is parsed using [@loaders.gl/las](https://loaders.gl/docs/modules/las/api-reference/las-loader) | data | ${loadingMS} | | 2. | X/Y coordinates of points are stored in a [spatial index (Flatbush)](https://github.com/mourner/flatbush) | index | ${indexingMS} | | 3. | A ${gridSize} * ${gridSize} grid is created. For each point in the grid, the spatial index is queried to find its closest neighbor. The elevation from that point is used in the final mesh. | terrain | ${terrainMS} | | 4. | a MARTINI mesh is created from the elevation grid | geometry | ${martiniMS} |32:T42d,bbox = { let minX = data.attributes.POSITION.value[0]; let minY = data.attributes.POSITION.value[1]; let minZ = data.attributes.POSITION.value[2]; let maxX = data.attributes.POSITION.value[0]; let maxY = data.attributes.POSITION.value[1]; let maxZ = data.attributes.POSITION.value[2]; for (let i = 0; i < data.attributes.POSITION.value.length; i += 3) { if (data.attributes.POSITION.value[i] < minX) minX = data.attributes.POSITION.value[i]; if (data.attributes.POSITION.value[i + 1] < minY) minY = data.attributes.POSITION.value[i + 1]; if (data.attributes.POSITION.value[i + 2] < minZ) minZ = data.attributes.POSITION.value[i + 2]; if (data.attributes.POSITION.value[i] > maxX) maxX = data.attributes.POSITION.value[i]; if (data.attributes.POSITION.value[i + 1] > maxY) maxY = data.attributes.POSITION.value[i + 1]; if (data.attributes.POSITION.value[i + 2] > maxZ) maxZ = data.attributes.POSITION.value[i + 2]; } return { minX, minY, minZ, maxX, maxY, maxZ }; }33:T71f,geometry = { const start = new Date(); const error = 0; const martini = new Martini(gridSize); const tile = martini.createTile(terrain); const mesh = tile.getMesh(error); const geometry = new THREE.BufferGeometry(); const vertices = new Float32Array((mesh.vertices.length / 2) * 3); const terrainExaggeration = 1; let index = 0; for (let i = 0; i < mesh.vertices.length / 2; i++) { let x = mesh.vertices[i * 2], y = mesh.vertices[i * 2 + 1]; vertices[index++] = x * aequidistance[0] + bbox.minX; vertices[index++] = y * aequidistance[1] + bbox.minY; vertices[index++] = terrain[x * gridSize + y] * terrainExaggeration; } geometry.setAttribute("position", new THREE.BufferAttribute(vertices, 3)); geometry.setIndex(new THREE.BufferAttribute(mesh.triangles, 1)); geometry.computeVertexNormals(); geometry.computeBoundingBox(); // calc uvs to correctly position texture let { min, max } = geometry.boundingBox; let offset = new THREE.Vector2(0 - min.x, 0 - min.y); let range = new THREE.Vector2(max.x - min.x, max.y - min.y); const position = geometry.attributes.position; const martiniEnd = new Date(); mutable martiniMS = martiniEnd - start; const uvs = []; for (let i = 0; i < position.count; i++) { const v3 = new THREE.Vector3().fromBufferAttribute(position, i); uvs.push((v3.x + offset.x) / range.x); uvs.push((v3.y + offset.y) / range.y); } geometry.setAttribute( "uv", new THREE.BufferAttribute(new Float32Array(uvs), 2) ); geometry.setAttribute( "uv2", new THREE.BufferAttribute(new Float32Array(uvs), 2) ); geometry.attributes.uv.needsUpdate = true; geometry.attributes.uv2.needsUpdate = true; const end = new Date(); console.log("Calculating UVs", end - martiniEnd); return geometry; }34:[{"id":"3811343421f8327b0965124fbb55a99a6683d0dba1448f1b100a216d6a90b23a3afa4dd7be0278a70a5cd6986377aca90b8eeddd200615486464a951b8e4db1e","href":"https://static.observableusercontent.com/files/3811343421f8327b0965124fbb55a99a6683d0dba1448f1b100a216d6a90b23a3afa4dd7be0278a70a5cd6986377aca90b8eeddd200615486464a951b8e4db1e","name":"nrw_grid@3.csv","type":"text/csv","lastModified":1684503519947.49,"size":1147524}] 2f:{"title":"Creating a 3D mesh from LIDAR point cloud data","cells":[{"id":0,"value":"# Creating a 3D mesh from LIDAR point cloud data","pinned":false,"hidden":false,"mode":"md"},{"id":721,"value":"viewof k = Inputs.range([7, 12], {\n label: \"LAZ interpolation (tile size)\",\n value: 8,\n step: 1,\n format: (v) => Math.pow(2, v)\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1225,"value":"viewof textureSize = Inputs.range([7, 12], {\n label: \"Texture size\",\n value: 11,\n step: 1,\n format: (v) => Math.pow(2, v) + 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1294,"value":"viewof poi = Inputs.select(pois, {\n label: \"Select one\",\n value: pois[1],\n format: (r) => r.name\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":589,"value":"renderer.domElement","pinned":false,"hidden":false,"mode":"ojs"},{"id":1148,"value":"## Overview","pinned":false,"hidden":false,"mode":"md"},{"id":1150,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":1236,"value":"md`LAZ Point cloud data and satellite imagery from [opengeodata.nrw.de](https://www.opengeodata.nrw.de/produkte/) \n*Datenlizenz Deutschland - Zero - Version 2.0*`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1107,"value":"md`## Data & Methods`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1190,"value":"lazUrl = `https://www.opengeodata.nrw.de/produkte/geobasis/hm/bdom50_las/bdom50_las/${curTile.tile}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1401,"value":"ab = fetch(lazUrl).then((res) => res.arrayBuffer())","pinned":true,"hidden":false,"mode":"ojs"},{"id":1413,"value":"ab.byteLength","pinned":true,"hidden":false,"mode":"ojs"},{"id":1200,"value":"textureWMSUrl = `https://www.wms.nrw.de/geobasis/wms_nw_dop?language=ger&width=${\n Math.pow(2, textureSize) + 1\n}&height=${Math.pow(2, textureSize) + 1}}&bbox=${bbox.minX},${bbox.minY},${\n bbox.maxX\n},${\n bbox.maxY\n}&crs=EPSG%3A25832&format=image%2Fpng&request=GetMap&service=WMS&transparent=true&version=1.3.0&layers=nw_dop_rgb`","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"data = {\n let start = new Date();\n const data = await load(ab, LASWorkerLoader, {\n las: {\n skip: 4\n }\n });\n let end = new Date();\n mutable loadingMS = end - start;\n return data;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":509,"value":"tileSize = Math.pow(2, k)","pinned":false,"hidden":false,"mode":"ojs"},{"id":511,"value":"gridSize = tileSize+1","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"height = 420","pinned":false,"hidden":false,"mode":"ojs"},{"id":624,"value":"aequidistance = {\n const aequidistantFactorX = (bbox.maxX - bbox.minX) / tileSize;\n const aequidistantFactorY = (bbox.maxY - bbox.minY) / tileSize;\n return [aequidistantFactorX, aequidistantFactorY];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":863,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":892,"value":"center = ({ x: calcCenter(\"X\"), y: calcCenter(\"Y\"), z: calcCenter(\"Z\") })","pinned":false,"hidden":false,"mode":"ojs"},{"id":889,"value":"function calcCenter(dim) {\n return (bbox[\"max\" + dim] - bbox[\"min\" + dim]) / 2 + bbox[\"min\" + dim];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":122,"value":"LASWorkerLoader = (\n await import(\"https://cdn.skypack.dev/@loaders.gl/las@3.3.3?min\")\n).LASWorkerLoader","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"load = (await import(\"https://cdn.skypack.dev/@loaders.gl/core@3.3.3?min\")).load","pinned":false,"hidden":false,"mode":"ojs"},{"id":975,"value":"loadTexture = {\n const loader = new THREE.TextureLoader();\n return (url) => new Promise((resolve) => loader.load(url, resolve));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1164,"value":"mutable loadingMS = 0;","pinned":false,"hidden":false,"mode":"ojs"},{"id":1173,"value":"mutable indexingMS = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":1176,"value":"mutable terrainMS = 0;","pinned":false,"hidden":false,"mode":"ojs"},{"id":1178,"value":"mutable martiniMS = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":1347,"value":"mutable textureMS = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":1261,"value":"curTile = getETRS89TileForPoint(poi.coordinates)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1256,"value":"function getETRS89TileForPoint(coordinates) {\n // reproject point from WGS84 to ETRS89\n const [x, y] = proj4(\n \"+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs\",\n coordinates\n );\n //return [x, y];\n // find tile in grid that contains point\n const tile = nrw_grid_etrs89.find(\n (t) => x >= t.x && x < t.x + 1000 && y >= t.y && y < t.y + 1000\n );\n return tile;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1244,"value":"nrw_grid_etrs89 = (\n await fetch(\n \"https://www.opengeodata.nrw.de/produkte/geobasis/hm/bdom50_las/bdom50_las/index.json\"\n ).then((r) => r.json())\n).datasets[0].files.map((t) => ({\n x: +t.name.split(\"bdom50_\")[1].split(\"_\")[0].split(\"32\")[1] * 1000,\n y: +t.name.split(\"bdom50_\")[1].split(\"_\")[1].split(\"_\")[0] * 1000,\n tile: t.name\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1290,"value":"pois = [\n {\n name: \"Eggeberg\",\n coordinates: [8.351943, 52.08824]\n },\n {\n name: \"Eggeberg 2\",\n coordinates: [8.3386301, 52.0912863]\n },\n {\n name: \"Kölner Dom\",\n coordinates: [6.9579167, 50.9413232]\n },\n {\n name: \"Externsteine\",\n coordinates: [8.9214436, 51.8689634]\n },\n {\n name: \"Zeche Zollverein\",\n coordinates: [7.0413369, 51.4877427]\n },\n {\n name: \"Drachenfels\",\n coordinates: [7.2056338, 50.6679878]\n }\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1110,"value":"md`## THREE.js & Flatbush`","pinned":false,"hidden":false,"mode":"ojs"},{"id":784,"value":"index = {\n const start = new Date();\n const points = data.attributes.POSITION.value;\n const index = new Flatbush(points.length / 3);\n\n for (let i = 0; i < points.length; i += 3) {\n index.add(points[i], points[i + 1], points[i], points[i + 1]);\n }\n // perform the indexing\n index.finish();\n const end = new Date();\n mutable indexingMS = end - start;\n return index;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":789,"value":"terrain = {\n const start = new Date();\n const arr = new Float32Array(gridSize * gridSize);\n let ti = 0;\n for (let x = 0; x < gridSize; ++x) {\n for (let y = 0; y < gridSize; ++y) {\n const xPos = bbox.minX + aequidistance[0] * x;\n const yPos = bbox.minY + aequidistance[1] * y;\n const i = index.neighbors(xPos, yPos, 1)[0];\n arr[ti++] = data.attributes.POSITION.value[i * 3 + 2];\n }\n }\n const end = new Date();\n mutable terrainMS = end - start;\n return arr;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":476,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":28,"value":"camera = {\n const camera = new THREE.PerspectiveCamera(\n 100,\n Math.min(width, 640) / height,\n 0.1,\n 100000\n );\n camera.up.set(0, 0, 1);\n\n camera.position.set(center.x, center.y, center.z + 10);\n\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":549,"value":"scene = {\n const scene = new THREE.Scene();\n scene.background = new THREE.Color(0xffffff);\n\n var light = new THREE.DirectionalLight(0xffffff, 1.5);\n light.position.setScalar(100);\n scene.add(light);\n scene.add(new THREE.AmbientLight(0xffffff, 0.5));\n\n // add geometry\n\n const material = new THREE.MeshNormalMaterial({\n flatShading: true,\n side: THREE.DoubleSide\n });\n\n let threeMesh = new THREE.Mesh(\n geometry,\n new THREE.MeshBasicMaterial({ map: map, side: THREE.DoubleSide })\n );\n scene.add(threeMesh);\n\n return scene;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1352,"value":"map = {\n const start = new Date();\n const map = await loadTexture(textureWMSUrl);\n const end = new Date();\n mutable textureMS = end - start;\n return map;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":545,"value":"renderer = {\n const renderer = new THREE.WebGLRenderer({ antialias: true });\n const controls = new THREE.OrbitControls(camera, renderer.domElement);\n\n controls.screenSpacePanning = true;\n\n controls.target.set(center.x, center.y, center.z);\n renderer.setSize(Math.min(width, 640), height);\n renderer.setPixelRatio(devicePixelRatio);\n controls.addEventListener(\"change\", () => renderer.render(scene, camera));\n renderer.render(scene, camera);\n return renderer;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1086,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1253,"value":"proj4 = (await import(\"https://cdn.skypack.dev/proj4@2.9.0?min\")).default","pinned":false,"hidden":false,"mode":"ojs"},{"id":505,"value":"Martini = require(\"@mapbox/martini@0.2.0/martini.min.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":779,"value":"Flatbush = (await import(\"https://unpkg.com/flatbush@4.1.0/index.js?module\"))\n .default","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"THREE = {\n const THREE = (window.THREE = await require(\"three@0.141\"));\n await require(\"three@0.141/examples/js/controls/OrbitControls.js\").catch(\n () => {}\n );\n return THREE;\n}","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}