Merge pull request #10 from open-pv/feature-4-diffuse-radiation-from-…

…file

Feature 4 diffuse radiation and elevation

.github/workflows/tests.yml

@@ -2,7 +2,7 @@ name: Unit tests
 
 on:
   pull_request:
-    branches: [ main ]
+    branches: [main]
   workflow_dispatch:
 
 jobs:
@@ -15,15 +15,15 @@ jobs:
         node-version: [18.x, 20.x]
 
     steps:
-    - uses: actions/checkout@v4
+      - uses: actions/checkout@v4
 
-    - name: Use Node.js ${{ matrix.node-version }}
-      uses: actions/setup-node@v4
-      with:
-        node-version: ${{ matrix.node-version }}
+      - name: Use Node.js ${{ matrix.node-version }}
+        uses: actions/setup-node@v4
+        with:
+          node-version: ${{ matrix.node-version }}
 
-    - name: Install dependencies
-      run: yarn install
+      - name: Install dependencies
+        run: yarn install
 
-    - name: Run Vitest
-      run: yarn test
+      - name: Run Vitest
+        run: yarn test

package.json

@@ -3,7 +3,6 @@
   "version": "0.0.1a1",
   "description": "Simulating Shadows for PV Potential Analysis on 3D Data on the GPU.",
   "main": "./dist/index.js",
-  "module": "./dist/index.mjs",
   "type": "module",
   "types": "./dist/index.d.ts",
   "files": [

src/elevation.ts

@@ -0,0 +1,81 @@
+import { CartesianPoint, SphericalPoint } from './utils';
+
+export function fillMissingAltitudes(maxAngles: SphericalPoint[]): void {
+  // First copy the maxAngles to a newAngles list, so that changes
+  // in the list do not affect the algorithm
+  let newAngles = maxAngles.map((angle) => ({ ...angle }));
+  for (let i = 0; i < newAngles.length; i++) {
+    if (newAngles[i].altitude != -Infinity) {
+      continue;
+    }
+    let distance = 1;
+    while (true) {
+      let prevIndex = (i - distance + newAngles.length) % newAngles.length;
+      let nextIndex = (i + distance) % newAngles.length;
+
+      if (maxAngles[nextIndex].altitude !== -Infinity) {
+        newAngles[i].altitude = maxAngles[nextIndex].altitude;
+        break;
+      } else if (maxAngles[prevIndex].altitude !== -Infinity) {
+        newAngles[i].altitude = maxAngles[prevIndex].altitude;
+        break;
+      } else distance++;
+    }
+  }
+  // Overwrite the maxAngles to make changes in this vector global
+  for (let i = 0; i < maxAngles.length; i++) {
+    maxAngles[i] = newAngles[i];
+  }
+}
+
+/**
+ * Returns the vector from start to end in the Horizontal coordinate system
+ * @param start
+ * @param end
+ * @returns
+ */
+export function calculateSphericalCoordinates(start: CartesianPoint, end: CartesianPoint): SphericalPoint {
+  const dx = end.x - start.x;
+  const dy = end.y - start.y;
+  const dz = end.z - start.z;
+  if (dx == 0 && dy == 0) {
+    return { radius: 1, azimuth: 0, altitude: 0 };
+  }
+
+  const r = Math.sqrt(dx * dx + dy * dy + dz * dz);
+  const altitude = Math.asin(dz / r);
+
+  let azimuth = (2 * Math.PI - Math.atan2(dy, dx)) % (2 * Math.PI);
+
+  return { radius: 1, azimuth, altitude };
+}
+
+/**
+ * Calculates the maximum heights visible from an observer in a set of directions.
+ * Returns a list of spherical points of length numDirections.
+ * @param elevation list of points with x,y,z component
+ * @param observer Point of interest for which the elevation angles are calculated.
+ * @param numDirections Number of steps for the azimuth angle.
+ * @returns
+ */
+export function getMaxElevationAngles(
+  elevation: CartesianPoint[],
+  observer: CartesianPoint,
+  numDirections: number,
+): SphericalPoint[] {
+  let maxAngles: SphericalPoint[] = Array.from({ length: numDirections }, (_, index) => ({
+    radius: 1,
+    azimuth: index * ((2 * Math.PI) / numDirections),
+    altitude: -Infinity,
+  }));
+
+  for (let point of elevation) {
+    const { azimuth, altitude } = calculateSphericalCoordinates(observer, point);
+    const closestIndex = Math.round(azimuth / ((2 * Math.PI) / numDirections)) % numDirections;
+    if (altitude > maxAngles[closestIndex].altitude) {
+      maxAngles[closestIndex].altitude = altitude;
+    }
+  }
+  fillMissingAltitudes(maxAngles);
+  return maxAngles;
+}

src/index.ts

@@ -2,37 +2,47 @@ import * as THREE from 'three';
 import { BufferAttribute, BufferGeometry, TypedArray } from 'three';
 import * as BufferGeometryUtils from 'three/examples/jsm/utils/BufferGeometryUtils.js';
 import { viridis } from './colormaps';
-import { getRandomSunVectors } from './sun';
+import * as elevation from './elevation';
+import * as sun from './sun';
 import * as triangleUtils from './triangleUtils.js';
-import { ArrayType, Triangle } from './triangleUtils.js';
+import { CartesianPoint, Point, SphericalPoint, SunVector, isValidUrl } from './utils';
 
 // @ts-ignore
 import { rayTracingWebGL } from './rayTracingWebGL.js';
 
 /**
  * This class holds all information about the scene that is simulated.
- * A scene is typically equipped with the following attributes:
+ * A ShadingScene is typically equipped with the following attributes:
  * * A pair of coordinates to locate the scene
  * * Simulation geometries, where the PV potential is calculated
  * * Shading geometries, where no PV potential is calculated but which are
  *   responsible for shading
  */
-export default class Scene {
+export default class ShadingScene {
   simulationGeometries: Array<BufferGeometry>;
   shadingGeometries: Array<BufferGeometry>;
+  elevationRaster: Array<CartesianPoint>;
+  elevationRasterMidpoint: CartesianPoint;
   latitude: number;
   longitude: number;
+  elevationAzimuthDivisions: number;
 
   /**
    *
    * @param latitude Latitude of the midpoint of the scene.
    * @param longitude Longitude of the midpoint of the scene.
    */
   constructor(latitude: number, longitude: number) {
+    if (latitude === undefined || longitude === undefined) {
+      throw new Error('Latitude and Longitude must be defined');
+    }
     this.simulationGeometries = [];
     this.shadingGeometries = [];
+    this.elevationRaster = [];
+    this.elevationRasterMidpoint = { x: 0, y: 0, z: 0 };
     this.latitude = latitude;
     this.longitude = longitude;
+    this.elevationAzimuthDivisions = 60;
   }
 
   /**
@@ -56,6 +66,19 @@ export default class Scene {
   addShadingGeometry(geometry: BufferGeometry) {
     this.shadingGeometries.push(geometry);
   }
+  /**
+   * IMPORTANT: Make sure that the DEM and the building mesh are in the same units, for example 1 step in
+   * DEM coordinates should be equal to 1 step in the SimulationGeometry coordinates.
+   * @param raster List of Points with x,y,z coordinates, representing a digital elevation model (DEM)
+   * @param midpoint The point of the observer, ie the center of the building
+   * @param azimuthDivisions Number of divisions of the azimuth Angle, i.e. the list of the azimuth
+   * angle will be [0, ..., 2Pi] where the list has a lenght of azimuthDivisions
+   */
+  addElevationRaster(raster: CartesianPoint[], midpoint: CartesianPoint, azimuthDivisions: number) {
+    this.elevationAzimuthDivisions = azimuthDivisions;
+    this.elevationRaster = raster;
+    this.elevationRasterMidpoint = midpoint;
+  }
 
   /** @ignore */
   refineMesh(mesh: BufferGeometry, maxLength: number): BufferGeometry {
@@ -94,14 +117,17 @@ export default class Scene {
    * @param numberSimulations Number of random sun positions that are used to calculate the PV yield
    * @returns
    */
+
   async calculate(
     numberSimulations: number = 80,
+    irradianceUrl: string | undefined,
     progressCallback: (progress: number, total: number) => void = (progress, total) =>
       console.log(`Progress: ${progress}/${total}%`),
   ) {
     console.log('Simulation package was called to calculate');
     let simulationGeometry = BufferGeometryUtils.mergeGeometries(this.simulationGeometries);
     let shadingGeometry = BufferGeometryUtils.mergeGeometries(this.shadingGeometries);
+
     // TODO: This breaks everything, why?
     simulationGeometry = this.refineMesh(simulationGeometry, 0.5); // TODO: make configurable
 
@@ -134,7 +160,16 @@ export default class Scene {
       }
     }
     // Compute unique intensities
-    const intensities = await this.rayTrace(midpointsArray, normalsArray, meshArray, numberSimulations, progressCallback);
+    console.log('Calling this.rayTrace');
+
+    const intensities = await this.rayTrace(
+      midpointsArray,
+      normalsArray,
+      meshArray,
+      numberSimulations,
+      irradianceUrl,
+      progressCallback,
+    );
 
     if (intensities === null) {
       throw new Error('Error raytracing in WebGL.');
@@ -153,10 +188,10 @@ export default class Scene {
       intensities[i] /= numberSimulations;
     }
 
-    return this.show(simulationGeometry, intensities);
+    return this.createMesh(simulationGeometry, intensities);
   }
   /** @ignore */
-  show(subdividedGeometry: BufferGeometry, intensities: Float32Array) {
+  createMesh(subdividedGeometry: BufferGeometry, intensities: Float32Array): THREE.Mesh {
     const Npoints = subdividedGeometry.attributes.position.array.length / 9;
     var newColors = new Float32Array(Npoints * 9);
     for (var i = 0; i < Npoints; i++) {
@@ -187,6 +222,8 @@ export default class Scene {
    * @param midpoints midpoints of triangles for which to calculate intensities
    * @param normals normals for each midpoint
    * @param meshArray array of vertices for the shading mesh
+   * @param numberSimulations number of random sun positions that are used for the simulation. Either numberSimulations or irradianceUrl need to be given.
+   * @param diffuseIrradianceUrl url where a 2D json of irradiance values lies. To generate such a json, visit https://github.com/open-pv/irradiance
    * @return
    * @memberof Scene
    */
@@ -195,9 +232,35 @@ export default class Scene {
     normals: TypedArray,
     meshArray: Float32Array,
     numberSimulations: number,
+    diffuseIrradianceUrl: string | undefined,
     progressCallback: (progress: number, total: number) => void,
   ) {
-    let sunDirections = getRandomSunVectors(numberSimulations, this.latitude, this.longitude);
-    return rayTracingWebGL(midpoints, normals, meshArray, sunDirections, progressCallback);
+    let directIrradiance: SunVector[] = [];
+    let diffuseIrradiance: SunVector[] = [];
+    let shadingElevationAngles: SphericalPoint[] = [];
+
+    if (typeof diffuseIrradianceUrl === 'string' && isValidUrl(diffuseIrradianceUrl)) {
+      const diffuseIrradianceSpherical = await sun.fetchIrradiance(diffuseIrradianceUrl, this.latitude, this.longitude);
+      diffuseIrradiance = sun.convertSpericalToEuclidian(diffuseIrradianceSpherical);
+    } else if (typeof diffuseIrradianceUrl != 'undefined') {
+      throw new Error('The given url for diffuse Irradiance is not valid.');
+    }
+    console.log('Calling getRandomSunVectors');
+    directIrradiance = sun.getRandomSunVectors(numberSimulations, this.latitude, this.longitude);
+    console.log(directIrradiance);
+    if (this.elevationRaster.length > 0) {
+      shadingElevationAngles = elevation.getMaxElevationAngles(
+        this.elevationRaster,
+        this.elevationRasterMidpoint,
+        this.elevationAzimuthDivisions,
+      );
+      sun.shadeIrradianceFromElevation(directIrradiance, shadingElevationAngles);
+      if (diffuseIrradiance.length > 0) {
+        sun.shadeIrradianceFromElevation(diffuseIrradiance, shadingElevationAngles);
+      }
+    }
+    console.log('Calling rayTracingWebGL');
+    normals = normals.filter((_, index) => index % 9 < 3);
+    return rayTracingWebGL(midpoints, normals, meshArray, directIrradiance, diffuseIrradiance, progressCallback);
   }
 }

src/rayTracingWebGL.ts

@@ -1,4 +1,5 @@
 import { TypedArray } from 'three';
+import { Point, SunVector } from './utils';
 
 function addToArray(ar1: Float32Array, ar2: Float32Array) {
   for (var i = 0; i < ar1.length; i++) {
@@ -10,7 +11,8 @@ export function rayTracingWebGL(
   pointsArray: TypedArray,
   normals: TypedArray,
   trianglesArray: TypedArray,
-  sunDirections: Float32Array,
+  directRadiance: SunVector[],
+  diffuseRadiance: SunVector[],
   progressCallback: (progress: number, total: number) => void,
 ): Float32Array | null {
   const N_TRIANGLES = trianglesArray.length / 9;
@@ -179,11 +181,20 @@ export function rayTracingWebGL(
 
   var colorCodedArray = null;
   var isShadowedArray = null;
-  for (var i = 0; i < sunDirections.length; i += 3) {
-    progressCallback(i/3, sunDirections.length/3);
+
+  for (var i = 0; i < directRadiance.length; i += 1) {
+    if (directRadiance[i].isShadedByElevation) {
+      continue;
+    }
+    progressCallback(i, directRadiance.length);
+
     // TODO: Iterate over sunDirection
     let sunDirectionUniformLocation = gl.getUniformLocation(program, 'u_sun_direction');
-    gl.uniform3fv(sunDirectionUniformLocation, [sunDirections[i], sunDirections[i + 1], sunDirections[i + 2]]);
+    gl.uniform3fv(sunDirectionUniformLocation, [
+      directRadiance[i].vector.cartesian.x,
+      directRadiance[i].vector.cartesian.y,
+      directRadiance[i].vector.cartesian.z,
+    ]);
 
     drawArraysWithTransformFeedback(gl, tf, gl.POINTS, N_POINTS);