at least something renders now

client/src/lib/Map.svelte

@@ -231,12 +231,8 @@
 
 						gl.useProgram(this.program);
 
-						console.log(drd.radials.length);
 						gl.uniform1i(gl.getUniformLocation(this.program, 'azimuthCount'), drd.radials.length);
-						console.log('done');
+
-						console.log(
-							gl.getUniform(this.program, gl.getUniformLocation(this.program, 'azimuthCount'))
-						);
 						gl.uniform1fv(
 							gl.getUniformLocation(this.program, 'azimuthAngles'),
 							new Float32Array(
@@ -272,49 +268,31 @@
 								}
 							}
 						}
-						const dataTexture = gl.createTexture();
+						const rdata: number[] = [];
-						gl.activeTexture(gl.TEXTURE0);
+						for (const radial of drd.radials) {
-						gl.bindTexture(gl.TEXTURE_2D, dataTexture);
+							for (const tdata of radial.product.data.data) {
-						const level = 0;
+								rdata.push(scaleMomentData(radial, tdata));
-						const internalFormat = gl.R32F;
-						const width = 720;
-						const height = 1832;
-						const border = 0;
-						const format = gl.RED;
-						const type = gl.FLOAT;
-						// interleave the data
-						const data: number[] = [];
-						for (let y = 0; y < height; y++) {
-							for (let x = 0; x < width; x++) {
-								const radial = drd.radials[x];
-								if (radial.product && radial.product.data) {
-									data.push(scaleMomentData(radial, radial.product.data.data[y]));
-								}
 							}
 						}
-						const typedData = new Float32Array(data);
+						const data = new Float32Array(rdata);
-						gl.texImage2D(
+
-							gl.TEXTURE_2D,
+						gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 4);
-							level,
+
-							internalFormat,
+						gl.activeTexture(gl.TEXTURE0 + 4);
-							width,
+						this.texture = gl.createTexture();
-							height,
+						gl.bindTexture(gl.TEXTURE_2D, this.texture);
-							border,
-							format,
-							type,
-							data
-						);
-						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
-						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
 						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE);
 						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);
-						const alignment = 1;
+						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
-						gl.pixelStorei(gl.UNPACK_ALIGNMENT, alignment);
+						gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
-
+						gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32F, 720, 1832, 0, gl.RED, gl.FLOAT, data);
-						gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 0);
 					},
 					render(gl, args) {
 						gl.useProgram(this.program);
+
+						gl.activeTexture(gl.TEXTURE0 + 4);
+						gl.bindTexture(gl.TEXTURE_2D, this.texture);
+
 						gl.uniformMatrix4fv(
 							gl.getUniformLocation(this.program, 'u_matrix'),
 							false,

client/src/lib/map/fragment.glsl

@@ -19,211 +19,6 @@ void xyToLngLat(in float x, in float y, out float lat, out float lng) {
 
 in vec4 raw_pos;
 
-struct Ellipsoid {
-    float semiMajorAxisMeters;
-    float semiMinorAxisMeters;
-    float flattening;
-    float inverseFlattening;
-};
-Ellipsoid fromAAndInverseF(float semiMajorAxisMeters, float inverseFlattening) {
-    float f = 1.0 / inverseFlattening;
-    float b = (1.0 - f) * semiMajorAxisMeters;
-
-    return Ellipsoid(semiMajorAxisMeters, b, f, inverseFlattening);
-}
-
-struct Angle {
-    float radians;
-};
-float angleAsDegrees(Angle angle) {
-    return degrees(angle.radians);
-}
-Angle radianAngle(float radians) {
-    return Angle(radians);
-}
-Angle degreeAngle(float degrees) {
-    return Angle(radians(degrees));
-}
-
-struct GlobalCoordinates {
-    Angle latitude;
-    Angle longitude;
-};
-void canonicalizeGlobalCoordinates(inout GlobalCoordinates coords) {
-    float latitudeRadians = coords.latitude.radians;
-    float longitudeRadians = coords.longitude.radians;
-
-    latitudeRadians = mod((latitudeRadians + PI), TWO_PI);
-    if (latitudeRadians < 0.0) {
-        latitudeRadians += TWO_PI;
-    }
-    latitudeRadians -= PI;
-
-    if (latitudeRadians > PI_OVER_TWO) {
-        latitudeRadians = PI - latitudeRadians;
-        longitudeRadians += PI;
-    } else if (latitudeRadians < NEGATIVE_PI_OVER_TWO) {
-        latitudeRadians = -PI - latitudeRadians;
-        longitudeRadians += PI;
-    }
-
-    longitudeRadians = mod((longitudeRadians + PI), TWO_PI);
-    if (longitudeRadians <= 0.0) {
-        longitudeRadians += TWO_PI;
-    }
-    longitudeRadians -= PI;
-
-    coords.latitude = radianAngle(latitudeRadians);
-    coords.longitude = radianAngle(longitudeRadians);
-}
-
-struct Geodedic {
-    float s;
-    Angle a1;
-    Angle a2;
-    bool didGetGoodEstimate;
-};
-
-const float tolerance = pow(10.0, -13.0);
-
-Geodedic vincenty(Ellipsoid ellipsoid, GlobalCoordinates start, GlobalCoordinates end) {
-    float a = ellipsoid.semiMajorAxisMeters;
-    float b = ellipsoid.semiMinorAxisMeters;
-    float f = ellipsoid.flattening;
-
-    float phi1 = start.latitude.radians;
-    float lambda1 = start.longitude.radians;
-    float phi2 = end.latitude.radians;
-    float lambda2 = end.longitude.radians;
-
-    float a2 = a * a;
-    float b2 = b * b;
-    float a2b2b2 = (a2 - b2) / b2;
-
-    float omega = lambda2 - lambda1;
-
-    float tanphi1 = tan(phi1);
-    float tanU1 = (1.0 - f) * tanphi1;
-    float U1 = atan(tanU1);
-    float sinU1 = sin(U1);
-    float cosU1 = cos(U1);
-
-    float tanphi2 = tan(phi2);
-    float tanU2 = (1.0 - f) * tanphi2;
-    float U2 = atan(tanU2);
-    float sinU2 = sin(U2);
-    float cosU2 = cos(U2);
-
-    float sinU1sinU2 = sinU1 * sinU2;
-    float cosU1sinU2 = cosU1 * sinU2;
-    float sinU1cosU2 = sinU1 * cosU2;
-    float cosU1cosU2 = cosU1 * cosU2;
-
-    // equ. 13
-    float lambda = omega;
-
-    // intermediates to compute 's'
-    float A = 0.0;
-    float B = 0.0;
-    float sigma = 0.0;
-    float deltasigma = 0.0;
-    float lambda0;
-    bool converged = false;
-
-    for (int i = 0; i < 20; i++) {
-        lambda0 = lambda;
-        float sinlambda = sin(lambda);
-        float coslambda = cos(lambda);
-        // equ. 14
-        float cosU1sinU2_sinU2cosU2coslambda = cosU1sinU2 - sinU1cosU2 * coslambda;
-        float sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + (cosU1sinU2_sinU2cosU2coslambda * cosU1sinU2_sinU2cosU2coslambda);
-        float sinsigma = sqrt(sin2sigma);
-
-        // equ. 15
-        float cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda);
-
-        // equ. 16
-        sigma = atan(sinsigma, cossigma);
-
-        // equ. 17 - careful, sin2sigma might be almost 0
-        float sinalpha = (sin2sigma == 0.0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma;
-        float alpha = asin(sinalpha);
-        float cosalpha = cos(alpha);
-        float cos2alpha = cosalpha * cosalpha;
-
-        // equ. 18 - careful, cos2alpha might be almost 0
-        float cos2sigmam = (cos2alpha == 0.0) ? 0.0 : cossigma - 2.0 * sinU1sinU2 / cos2alpha;
-        float u2 = cos2alpha * a2b2b2;
-        float cos2sigmam2 = cos2sigmam * cos2sigmam;
-
-        // equ. 3
-        A = 1.0 + u2 / 16384.0 * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
-        // equ. 4
-        B = u2 / 1024.0 * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
-        // equ. 6
-        deltasigma = B * sinsigma * (cos2sigmam + B / 4.0 * (cossigma * (-1.0 + 2.0 * cos2sigmam2) - B / 6.0 * cos2sigmam * (-3.0 + 4.0 * sin2sigma) * (-3.0 + 4.0 * cos2sigmam2)));
-        // equ. 10
-        float C = f / 16.0 * cos2alpha * (4.0 + f * (4.0 - 3.0 * cos2alpha));
-        // equ. 11 (modified)
-        lambda = omega + (1.0 - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1.0 + 2.0 * cos2sigmam2)));
-
-        if (i < 2) {
-            continue;
-        }
-
-        float change = abs((lambda - lambda0) / lambda);
-        if (change < tolerance) {
-            converged = true;
-            break;
-        }
-
-        // equ. 19
-        float s = b * A * (sigma - deltasigma);
-        Angle alpha1;
-        Angle alpha2;
-
-        bool didGetGoodEstimate = true;
-
-        // didn't converge? must be N/S
-        if (!converged) {
-            if (phi1 > phi2) {
-                alpha1 = degreeAngle(180.0);
-                alpha2 = degreeAngle(0.0);
-            } else if (phi1 < phi2) {
-                alpha1 = degreeAngle(0.0);
-                alpha2 = degreeAngle(180.0);
-            } else {
-                alpha1 = Angle(0.0 / 0.0); // NaN
-                alpha2 = Angle(0.0 / 0.0); // NaN
-                didGetGoodEstimate = false;
-            }
-        } else {
-            float radians;
-            // equ. 20
-            radians = atan(cosU2 * sin(lambda), (cosU1sinU2 - sinU1cosU2 * cos(lambda)));
-            if (radians < 0.0) {
-                radians += TWO_PI;
-            }
-            alpha1 = radianAngle(radians);
-
-            radians = atan(cosU1 * sin(lambda), (-sinU1cosU2 + cosU1sinU2 * cos(lambda))) + PI;
-            if (radians < 0.0) {
-                radians += TWO_PI;
-            }
-            alpha2 = radianAngle(radians);
-        }
-
-        if (alpha1.radians >= TWO_PI) {
-            alpha1 = radianAngle(alpha1.radians - TWO_PI);
-        }
-        if (alpha2.radians >= TWO_PI) {
-            alpha2 = radianAngle(alpha2.radians - TWO_PI);
-        }
-
-        return Geodedic(s, alpha1, alpha2, didGetGoodEstimate);
-    }
-}
-
 uniform int azimuthCount;
 uniform float[720] azimuthAngles;
 uniform float azimuthSpacing;
@@ -266,7 +61,6 @@ LocateRadialResult locateRadial(float forAzimuth) {
 }
 
 void main() {
-    Ellipsoid WGS84 = fromAAndInverseF(6378137.0, 298.257223563);
     float lat;
     float lng;
     xyToLngLat(raw_pos.x, raw_pos.y, lat, lng);
@@ -276,24 +70,24 @@ void main() {
         return;
     }
 
-    GlobalCoordinates radar = GlobalCoordinates(degreeAngle(radarLat), degreeAngle(radarLng));
+    float R = 6371.0 * pow(10.0, 3.0); // meters
-    canonicalizeGlobalCoordinates(radar);
+    float phi1 = radians(radarLat);
-    GlobalCoordinates
+    float phi2 = radians(lat);
-    samplePoint = GlobalCoordinates(degreeAngle(lat), degreeAngle(lng));
+    float lambda1 = radians(radarLng);
-    canonicalizeGlobalCoordinates(samplePoint);
+    float lambda2 = radians(lng);
 
-    Geodedic vincentyResult = vincenty(WGS84, radar, samplePoint);
+    float deltaPhi = radians(lat - radarLat);
+    float deltaLambda = radians(lng - radarLng);
+    float a = sin(deltaPhi / 2.0) * sin(deltaPhi / 2.0) + cos(phi1) * cos(phi2) * sin(deltaLambda / 2.0) * sin(deltaLambda / 2.0);
+    float c = 2.0 * atan(sqrt(a), sqrt(1.0 - a));
+    float d = R * c; // meters
 
-    if (!vincentyResult.didGetGoodEstimate) {
+    float d_m = d;
-        fragColor = vec4(0.0, 0.0, 0.0, 0.0);
-        return;
-    }
 
-    float d_m = vincentyResult.s;
+    float y = sin(lambda2 - lambda1) * cos(phi2);
-    float azimuth = degrees(vincentyResult.a1.radians);
+    float x = cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(lambda2 - lambda1);
-    if (azimuth < 0.0) {
+    float theta = atan(y, x);
-        azimuth += 360.0;
+    float azimuth = (theta * 180.0 / PI + 360.0); // degrees
-    }
 
     LocateRadialResult maybeRadial = locateRadial(azimuth);
     if (!maybeRadial.didFindRadial) {
@@ -306,7 +100,7 @@ void main() {
     float gate_spacing_km = sampleInterval / 1000.0;
 
     if (distance_km < first_gate_distance_km) {
-        fragColor = vec4(1.0, 0.0, 0.0, 1.0);
+        fragColor = vec4(0.0, 0.0, 0.0, 0.0);
         return;
     }
 
@@ -316,8 +110,27 @@ void main() {
         return;
     }
 
-    vec2 coords = vec2(float(maybeRadial.radialIndex), float(gate));
+    if (!maybeRadial.didFindRadial) {
-    float rawValue = texture(scaledData, coords).r;
+        fragColor = vec4(0.0, 0.0, 0.0, 0.0);
+        return;
+    }
+
+    ivec2 coords = ivec2(maybeRadial.radialIndex, gate);
+    float rawValue = texelFetch(scaledData, coords, 0).r;
+
+    const float BELOW_THRESHOLD = -9999.0;
+    const float RANGE_FOLDED = -9998.0;
+
+    if (rawValue == BELOW_THRESHOLD) {
+        fragColor = vec4(0.0, 0.0, 0.0, 0.0);
+        return;
+    }
+    if (rawValue == RANGE_FOLDED) {
+        fragColor = vec4(1.0, 0.0, 1.0, 1.0);
+        return;
+    }
+
+    //float rawValue = (float(value.r) - float(value.g)) / float(value.b);
 
     if (rawValue > 80.0) {
         fragColor = vec4(128.0 / 255.0, 128.0 / 255.0, 128.0 / 255.0, 1.0);