diff --git a/.idea/workspace.xml b/.idea/workspace.xml index 50ffb7e..2e1a892 100644 --- a/.idea/workspace.xml +++ b/.idea/workspace.xml @@ -7,11 +7,17 @@ - - + + + + + + + + - - + @@ -137,6 +154,7 @@ - \ No newline at end of file diff --git a/Cargo.lock b/Cargo.lock index d7067fd..cf10be2 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1761,6 +1761,15 @@ dependencies = [ "num-traits", ] +[[package]] +name = "ordered-float" +version = "5.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2c1f9f56e534ac6a9b8a4600bdf0f530fb393b5f393e7b4d03489c3cf0c3f01" +dependencies = [ + "num-traits", +] + [[package]] name = "overload" version = "0.1.1" @@ -3689,6 +3698,7 @@ name = "wxbox-interchange" version = "0.1.0" dependencies = [ "nexrad-decode", + "ordered-float 5.0.0", "prost", "prost-build", "prost-types", @@ -3706,7 +3716,7 @@ dependencies = [ name = "wxbox-pal" version = "0.1.0" dependencies = [ - "ordered-float", + "ordered-float 4.6.0", "thiserror 1.0.69", ] diff --git a/client/src/lib/Map.svelte b/client/src/lib/Map.svelte index c4d7d45..074646f 100644 --- a/client/src/lib/Map.svelte +++ b/client/src/lib/Map.svelte @@ -15,6 +15,8 @@ import vertexSource from '$lib/map/vertex.glsl?raw'; import { degreesToRadians } from '@turf/turf'; import type { DigitalRadarData, Radial } from './generated_interop/digitalRadarData'; + import {forwardGeodesic, radToDeg} from "$lib/vincenty"; + import {degToRad} from "$lib/vincenty.js"; const BELOW_THRESHOLD = -9999.0; const RANGE_FOLDED = -9998.0; @@ -181,6 +183,7 @@ gl.attachShader(this.program, vertexShader); gl.attachShader(this.program, fragmentShader); gl.linkProgram(this.program); + gl.useProgram(this.program); this.aPos = gl.getAttribLocation(this.program, 'a_pos'); @@ -191,69 +194,69 @@ 111.320 * cos(latitude) km = 1 deg */ - gl.useProgram(this.program); + const lat = 38.976111; - gl.uniform1f(gl.getUniformLocation(this.program, 'radarLat'), lat); const long = -77.4875; + + gl.uniform1f(gl.getUniformLocation(this.program, 'radarLat'), lat); gl.uniform1f(gl.getUniformLocation(this.program, 'radarLng'), long); + const vertexData: number[] = []; + const triangleAzimuthLookup: number[] = []; + + const radarCoordinate = maplibregl.MercatorCoordinate.fromLngLat({ + lng: long, + lat: lat + }); + + for (let i =0; i < drd.radials.length; i++) { + const radial = drd.radials[i]; + + if (radial.product && radial.product.data && radial.product.data.data) { + const angle_a = radial.azimuthAngleDegrees; + const angle_b = angle_a + 2*radial.azimuthSpacingDegrees; + + const start_range = radial.product.data.startRange; // m + const sample_interval = radial.product.data.sampleInterval; // m + const number_of_samples = radial.product.data.data.length; + const range = sample_interval * number_of_samples + start_range; // m + // add an extra sample for good measure + const padded_range = range + sample_interval; // m + // calculate the two points + const [pointALat, pointALong] = forwardGeodesic(degToRad(angle_a), degToRad(lat), degToRad(long), padded_range); + const pointA = maplibregl.MercatorCoordinate.fromLngLat({ + lng: radToDeg(pointALong), + lat: radToDeg(pointALat) + }); + const [pointBLat, pointBLong] = forwardGeodesic(degToRad(angle_b), degToRad(lat), degToRad(long), padded_range); + const pointB = maplibregl.MercatorCoordinate.fromLngLat({ + lng: radToDeg(pointBLong), + lat: radToDeg(pointBLat) + }); + vertexData.push(radarCoordinate.x); + vertexData.push(radarCoordinate.y); + vertexData.push(pointA.x); + vertexData.push(pointA.y); + vertexData.push(pointB.x); + vertexData.push(pointB.y); + triangleAzimuthLookup.push(angle_a); + } + } + const radar_range_maximum = 560; // ish km const degrees_of_lat = radar_range_maximum / 110.574; const degrees_of_long = radar_range_maximum / (111.32 * Math.cos(degreesToRadians(lat))); - // A-C - // B-D - const a = maplibregl.MercatorCoordinate.fromLngLat({ - lng: long - degrees_of_long, - lat: lat + degrees_of_lat - }); - const b = maplibregl.MercatorCoordinate.fromLngLat({ - lng: long - degrees_of_long, - lat: lat - degrees_of_lat - }); - const c = maplibregl.MercatorCoordinate.fromLngLat({ - lng: long + degrees_of_long, - lat: lat + degrees_of_lat - }); - const d = maplibregl.MercatorCoordinate.fromLngLat({ - lng: long + degrees_of_long, - lat: lat - degrees_of_lat - }); - this.buffer = gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, this.buffer); gl.bufferData( gl.ARRAY_BUFFER, - new Float32Array([a.x, a.y, b.x, b.y, c.x, c.y, b.x, b.y, d.x, d.y, c.x, c.y]), + new Float32Array(vertexData), gl.STATIC_DRAW ); - - gl.useProgram(this.program); - - gl.uniform1i(gl.getUniformLocation(this.program, 'azimuthCount'), drd.radials.length); - - gl.uniform1fv( - gl.getUniformLocation(this.program, 'azimuthAngles'), - new Float32Array( - drd.radials.map((u) => { - return u.azimuthAngleDegrees; - }) - ) - ); - gl.uniform1f( - gl.getUniformLocation(this.program, 'azimuthSpacing'), - drd.radials[0].azimuthSpacingDegrees - ); - gl.uniform1f( - gl.getUniformLocation(this.program, 'startRange'), - drd.radials[0].product?.data?.startRange - ); - console.log(drd.radials[0].product?.data?.startRange); - gl.uniform1f( - gl.getUniformLocation(this.program, 'sampleInterval'), - drd.radials[0].product?.data?.sampleInterval - ); + this.vertexCount = vertexData.length/2; + gl.uniform1fv(gl.getUniformLocation(this.program, 'triangleAzimuthLookup'), new Float32Array(triangleAzimuthLookup)); function scaleMomentData(radial: Radial, u: number): number { if (!radial.product || !radial.product.data) return BELOW_THRESHOLD; @@ -269,31 +272,34 @@ } } } - const rdata: number[] = []; - for (let i = 0; i < drd.radials.length; i++) { - const radial = drd.radials[drd.radials.length - 1 - i]; - for (const tdata of radial.product.data.data) { - rdata.push(scaleMomentData(radial, tdata)); + const raw_data = []; + let validRadials = 0; + + for (const radial of drd.radials) { + if (radial.product && radial.product.data && radial.product.data.data) { + for (const datapoint of radial.product.data.data) { + raw_data.push(scaleMomentData(radial, datapoint)); + } + validRadials++; } } - const data = new Float32Array(rdata); - gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 4); - - gl.activeTexture(gl.TEXTURE0 + 4); + const data = new Float32Array(raw_data); + gl.activeTexture(gl.TEXTURE0 + 10); this.texture = gl.createTexture(); gl.bindTexture(gl.TEXTURE_2D, this.texture); - 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); + gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32F, 1832, validRadials, 0, gl.RED, gl.FLOAT, data); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); 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.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); + gl.uniform1i(gl.getUniformLocation(this.program, 'data'), 10); + }, render(gl, args) { - gl.useProgram(this.program); - - gl.activeTexture(gl.TEXTURE0 + 4); + gl.activeTexture(gl.TEXTURE0 + 10); gl.bindTexture(gl.TEXTURE_2D, this.texture); + gl.useProgram(this.program); gl.uniformMatrix4fv( gl.getUniformLocation(this.program, 'u_matrix'), @@ -306,10 +312,10 @@ gl.vertexAttribPointer(this.aPos, 2, gl.FLOAT, false, 0, 0); gl.enable(gl.BLEND); gl.blendFunc(gl.SRC_ALPHA, gl.ONE_MINUS_SRC_ALPHA); - gl.drawArrays(gl.TRIANGLES, 0, 6); + gl.drawArrays(gl.TRIANGLES, 0, this.vertexCount); } }; - map.addLayer(dataLayer); + map.addLayer(dataLayer, 'alerts'); } }); map.on('error', (e) => { diff --git a/client/src/lib/map/fragment.glsl b/client/src/lib/map/fragment.glsl index 4f3ef87..4825006 100644 --- a/client/src/lib/map/fragment.glsl +++ b/client/src/lib/map/fragment.glsl @@ -2,135 +2,82 @@ precision highp float; -const float PI = 3.141592653589793238462643; -const float PI_OVER_TWO = PI / 2.0; -const float TWO_PI = PI + PI; -const float NEGATIVE_PI_OVER_TWO = -PI_OVER_TWO; -const float NEGATIVE_TWO_PI = -TWO_PI; +const float PI = 3.14159265358979323846264338327950288; -out highp vec4 fragColor; -uniform mat4 u_matrix; +layout(location = 0) out vec4 fragColor; -void xyToLngLat(in float x, in float y, out float lat, out float lng) { - lng = x * 360.0 - 180.0; - float y2 = 180.0 - y * 360.0; - lat = 360.0 / PI * atan(exp(y2 * PI / 180.0)) - 90.0; -} +in vec2 raw_pos; +in vec4 vertexColor; +flat in int azimuthNumber; +flat in float azimuthAngle; -in vec4 raw_pos; +uniform sampler2D data; -uniform int azimuthCount; -uniform float[720] azimuthAngles; -uniform float azimuthSpacing; -uniform sampler2D scaledData; -uniform float startRange; -uniform float sampleInterval; uniform float radarLat; uniform float radarLng; -struct LocateRadialResult { - int radialIndex; - float radialDistance; - bool didFindRadial; +struct Wgs84Coordinates { + float lat; + float lng; }; +Wgs84Coordinates rawPosToWgs84(vec2 raw_pos) { + float lng = raw_pos.x * 360.0 - 180.0; + float y2 = 180.0 - raw_pos.y * 360.0; + float lat = 360.0 / PI * atan(exp(y2 * PI / 180.0)) - 90.0; + return Wgs84Coordinates(lat,lng); +} -LocateRadialResult locateRadial(float forAzimuth) { - int closestRadial; - float bestDistance = 1.0; - bool foundAnything = false; +float haversineDistance(Wgs84Coordinates from, Wgs84Coordinates to, float radius) { + return 2.0 * radius * asin(sqrt(pow(2.0, sin((to.lat-from.lat)/2.0)) + cos(from.lat) * cos(to.lat) * pow(2.0, sin((to.lng - from.lng)/2.0)))); +} - for (int i = 0; i < azimuthCount; i++) { - float angle = azimuthAngles[i]; - float this_dist = abs(angle - forAzimuth); - if (this_dist > azimuthSpacing) { - continue; - } - if (foundAnything) { - if (this_dist < bestDistance) { - closestRadial = i; - bestDistance = this_dist; - } - } else { - closestRadial = i; - bestDistance = this_dist; - foundAnything = true; - } - } +const float R = 6371000.0; // m - return LocateRadialResult(closestRadial, bestDistance, foundAnything); +float random (vec2 st) { + return fract(sin(dot(st.xy, + vec2(12.9898,78.233)))* + 43758.5453123); } void main() { - float lat; - float lng; - xyToLngLat(raw_pos.x, raw_pos.y, lat, lng); + Wgs84Coordinates thisPixelCoordinates = rawPosToWgs84(raw_pos); + Wgs84Coordinates radarCoordinates = Wgs84Coordinates(radarLat,radarLng); - if (azimuthCount == 0) { - fragColor = vec4(0.0, 0.0, 0.0, 0.0); + float phi1 = radians(radarCoordinates.lat); + float phi2 = radians(thisPixelCoordinates.lat); + float deltaPhi = phi2 - phi1; + float deltaLambda = radians(thisPixelCoordinates.lng - radarCoordinates.lng); + 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; + + float distance_meters = d; + + float start_distance = 2125.0; + float interval = 250.0; + int number_of_gates = 1832; + + if (distance_meters < start_distance) { + fragColor = vec4(0,0,0,0); return; } - - float R = 6371000.0; // meters - float phi1 = radians(radarLat); - float phi2 = radians(lat); - float lambda2 = radians(lng); - float lambda1 = radians(radarLng); - float deltaLambda = radians(lng - radarLng); - float d_m = acos(sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2) * cos(deltaLambda)) * R; - - float distance_km = d_m / 1000.0; - float first_gate_distance_km = startRange / 1000.0; - float gate_spacing_km = sampleInterval / 1000.0; - int gate = int(round((distance_km - first_gate_distance_km) / gate_spacing_km)); - - float y = sin(lambda2 - lambda1) * cos(phi2); - float x = cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(lambda2 - lambda1); - float theta = atan(y, x); - float azimuth = mod((theta * 180.0 / PI + 360.0), 360.0); - - if (distance_km < 10.0) { - fragColor = vec4(0.0,1.0,0.0,1.0); - return; - } - - if (distance_km < first_gate_distance_km || gate < 0) { - fragColor = vec4(0.0, 0.0, 0.0, 0.0); - return; - } - - if (gate > 1832) { - fragColor = vec4(0.0, 0.0, 0.0, 0.0); - return; - } - - fragColor = vec4(float(gate) / 1832.0, 0.0, 0.0, 1.0); - //fragColor = vec4(float(gate) / 1832.0, 0.0, azimuth / 360.0, 1.0); - return; - - LocateRadialResult maybeRadial = locateRadial(azimuth); - if (!maybeRadial.didFindRadial) { - fragColor = vec4(0.0, 0.0, 0.0, 0.0); + if (distance_meters > (start_distance + interval * float(number_of_gates))) { + fragColor = vec4(0,0,0,0); return; } + int selector = int(mod(float(azimuthNumber), 2.0)); - - - 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; + if (selector == 0) { + fragColor = vec4(float(azimuthAngle) / 360.0, 0, 0, 1); + } else if (selector == 1) { + fragColor = vec4(0,float(azimuthAngle) / 360.0, 0,1); + } else { + fragColor = vec4(0,0,float(azimuthAngle)/360.0,1); } - //float rawValue = (float(value.r) - float(value.g)) / float(value.b); + int gate_number = int(floor((distance_meters - 2125.0) / 250.0)); + + float rawValue = texelFetch(data, ivec2(gate_number, azimuthNumber), 0).r; if (rawValue > 80.0) { fragColor = vec4(128.0 / 255.0, 128.0 / 255.0, 128.0 / 255.0, 1.0); @@ -151,4 +98,4 @@ void main() { } else { fragColor = vec4(0.0, 0.0, 0.0, 0.0); } -} +} \ No newline at end of file diff --git a/client/src/lib/map/fragment_old.glsl b/client/src/lib/map/fragment_old.glsl new file mode 100644 index 0000000..b4c169d --- /dev/null +++ b/client/src/lib/map/fragment_old.glsl @@ -0,0 +1,154 @@ +#version 300 es + +precision highp float; + +const float PI = 3.141592653589793238462643; +const float PI_OVER_TWO = PI / 2.0; +const float TWO_PI = PI + PI; +const float NEGATIVE_PI_OVER_TWO = -PI_OVER_TWO; +const float NEGATIVE_TWO_PI = -TWO_PI; + +out highp vec4 fragColor; +uniform mat4 u_matrix; + +void xyToLngLat(in float x, in float y, out float lat, out float lng) { + lng = x * 360.0 - 180.0; + float y2 = 180.0 - y * 360.0; + lat = 360.0 / PI * atan(exp(y2 * PI / 180.0)) - 90.0; +} + +in vec4 raw_pos; + +uniform int azimuthCount; +uniform float[720] azimuthAngles; +uniform float azimuthSpacing; +uniform sampler2D scaledData; +uniform float startRange; +uniform float sampleInterval; +uniform float radarLat; +uniform float radarLng; + +struct LocateRadialResult { + int radialIndex; + float radialDistance; + bool didFindRadial; +}; + +LocateRadialResult locateRadial(float forAzimuth) { + int closestRadial; + float bestDistance = 1.0; + bool foundAnything = false; + + for (int i = 0; i < azimuthCount; i++) { + float angle = azimuthAngles[i]; + float this_dist = abs(angle - forAzimuth); + if (this_dist > azimuthSpacing) { + continue; + } + if (foundAnything) { + if (this_dist < bestDistance) { + closestRadial = i; + bestDistance = this_dist; + } + } else { + closestRadial = i; + bestDistance = this_dist; + foundAnything = true; + } + } + + return LocateRadialResult(closestRadial, bestDistance, foundAnything); +} + +void main() { + float lat; + float lng; + xyToLngLat(raw_pos.x, raw_pos.y, lat, lng); + + if (azimuthCount == 0) { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + return; + } + + float R = 6371000.0; // meters + float phi1 = radians(radarLat); + float phi2 = radians(lat); + float lambda2 = radians(lng); + float lambda1 = radians(radarLng); + float deltaLambda = radians(lng - radarLng); + float d_m = acos(sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2) * cos(deltaLambda)) * R; + + float distance_km = d_m / 1000.0; + float first_gate_distance_km = startRange / 1000.0; + float gate_spacing_km = sampleInterval / 1000.0; + int gate = int(round((distance_km - first_gate_distance_km) / gate_spacing_km)); + + float y = sin(lambda2 - lambda1) * cos(phi2); + float x = cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(lambda2 - lambda1); + float theta = atan(y, x); + float azimuth = mod((theta * 180.0 / PI + 360.0), 360.0); + + if (distance_km - 10.0 < 0.01) { + fragColor = vec4(0.0,1.0,0.0,mix(0.0, 1.0,distance_km-10.0)); + return; + } + + if (distance_km < first_gate_distance_km || gate < 0) { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + return; + } + + if (gate > 1832) { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + return; + } + + fragColor = vec4(float(gate) / 1832.0, 0.0, 0.0, 1.0); + //fragColor = vec4(float(gate) / 1832.0, 0.0, azimuth / 360.0, 1.0); + return; + + LocateRadialResult maybeRadial = locateRadial(azimuth); + if (!maybeRadial.didFindRadial) { + 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); + } else if (rawValue > 70.0) { + fragColor = vec4(1.0, 1.0, 1.0, 1.0); + } else if (rawValue > 60.0) { + fragColor = vec4(1.0, 0.0, 1.0, 1.0); + } else if (rawValue > 50.0) { + fragColor = vec4(1.0, 0.0, 0.0, 1.0); + } else if (rawValue > 40.0) { + fragColor = vec4(1.0, 1.0, 0.0, 1.0); + } else if (rawValue > 30.0) { + fragColor = vec4(0.0, 1.0, 0.0, 1.0); + } else if (rawValue > 20.0) { + fragColor = vec4(64.0 / 255.0, 128.0 / 255.0, 1.0, 1.0); + } else if (rawValue > 10.0) { + fragColor = vec4(164.0 / 255.0, 164.0 / 255.0, 1.0, 1.0); + } else { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + } +} diff --git a/client/src/lib/map/vertex.glsl b/client/src/lib/map/vertex.glsl index a8584cf..a2bff5f 100644 --- a/client/src/lib/map/vertex.glsl +++ b/client/src/lib/map/vertex.glsl @@ -5,9 +5,25 @@ precision highp float; uniform mat4 u_matrix; in vec2 a_pos; -out vec4 raw_pos; +out vec2 raw_pos; +out vec4 vertexColor; +flat out int azimuthNumber; +flat out float azimuthAngle; + +uniform float[720] triangleAzimuthLookup; + +float rand(vec2 co){ + return fract(sin(dot(co, vec2(12.9898, 78.233))) * 43758.5453); +} void main() { - raw_pos = vec4(a_pos, 0.0, 1.0); + raw_pos = a_pos; gl_Position = u_matrix * vec4(a_pos, 0.0, 1.0); + azimuthAngle = triangleAzimuthLookup[gl_VertexID / 3]; + azimuthNumber = gl_VertexID / 3; + + float randColor = rand(vec2(triangleAzimuthLookup[gl_VertexID/3] / 720.0,triangleAzimuthLookup[gl_VertexID/3] / 720.0)); + float randColor2 = rand(vec2(randColor, randColor)); + float randColor3 = rand(vec2(randColor, randColor)); + vertexColor = vec4(randColor,randColor2,randColor3,1.0); } \ No newline at end of file diff --git a/client/src/lib/vincenty.ts b/client/src/lib/vincenty.ts new file mode 100644 index 0000000..cd1b117 --- /dev/null +++ b/client/src/lib/vincenty.ts @@ -0,0 +1,48 @@ +const a = 6378137.0; +const b = 6356752.314245; +const f = (a-b) / a; + +export function degToRad(deg: number): number { + return deg * (Math.PI / 180); +} +export function radToDeg(radian: number): number { + return radian * (180 / Math.PI); +} + + +export function forwardGeodesic(α1: number, φ1: number, λ1: number, s: number): [number, number, number] { + const sinα1 = Math.sin(α1); + const cosα1 = Math.cos(α1); + + const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; + const σ1 = Math.atan2(tanU1, cosα1); // σ1 = angular distance on the sphere from the equator to P1 + const sinα = cosU1 * sinα1; // α = azimuth of the geodesic at the equator + const cosSqα = 1 - sinα*sinα; + const uSq = cosSqα * (a*a - b*b) / (b*b); + const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); + const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); + + let σ = s / (b*A), sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere + let cos2σₘ = null; // σₘ = angular distance on the sphere from the equator to the midpoint of the line + + let σʹ = null; + do { + cos2σₘ = Math.cos(2*σ1 + σ); + sinσ = Math.sin(σ); + cosσ = Math.cos(σ); + const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ))); + σʹ = σ; + σ = s / (b*A) + Δσ; + } while (Math.abs(σ-σʹ) > 1e-12); + + const x = sinU1*sinσ - cosU1*cosσ*cosα1; + const φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + x*x)); + const λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1); + const C = f/16*cosSqα*(4+f*(4-3*cosSqα)); + const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ))); + const λ2 = λ1 + L; + + const α2 = Math.atan2(sinα, -x); // final bearing + + return [φ2, λ2, α2]; +} \ No newline at end of file diff --git a/crates/interchange/Cargo.toml b/crates/interchange/Cargo.toml index d560143..a6bd709 100644 --- a/crates/interchange/Cargo.toml +++ b/crates/interchange/Cargo.toml @@ -8,11 +8,12 @@ prost = "0.13" prost-types = "0.13" wxbox-ar2 = { path = "../ar2", optional = true } +ordered-float = { version = "5", optional = true } nexrad-decode = { path = "../nexrad-decode", optional = true } [features] default = [] -ar2 = ["dep:wxbox-ar2", "dep:nexrad-decode"] +ar2 = ["dep:wxbox-ar2", "dep:nexrad-decode", "dep:ordered-float" ] [build-dependencies] prost-build = "0.13" \ No newline at end of file diff --git a/crates/interchange/src/ar2.rs b/crates/interchange/src/ar2.rs index daae479..aa54b05 100644 --- a/crates/interchange/src/ar2.rs +++ b/crates/interchange/src/ar2.rs @@ -1,3 +1,4 @@ +use ordered_float::OrderedFloat; use crate::cif::cif_container::MessageType; use crate::{AsCif, cif}; use nexrad_decode::messages::digital_radar_data::RadialStatus; @@ -16,13 +17,15 @@ impl AsCif for Scan { digital_radar_data.vcp_number = self.coverage_pattern_number as i32; // find the elevation - let maybe_elevation = self + let mut maybe_elevation = self .sweeps .iter() .find(|u| u.elevation_number == params.requested_elevation); if let Some(elevation) = maybe_elevation { // parse out the radials let mut radials_we_can_use: Vec<(&Radial, &MomentData)> = vec![]; + let mut elevation = elevation.clone(); + elevation.radials.sort_by_key(|u| OrderedFloat(u.azimuth_angle_degrees)); for radial in &elevation.radials { match params.requested_product.as_str() { "REF" if radial.reflectivity.is_some() => { diff --git a/crates/tiler/Cargo.toml b/crates/tiler/Cargo.toml index d0e166a..4ea5ca3 100644 --- a/crates/tiler/Cargo.toml +++ b/crates/tiler/Cargo.toml @@ -46,4 +46,4 @@ wxbox-common = { path = "../common" } chrono = "0.4" quick-xml = { version = "0.37", features = ["serialize"] } -wxbox-interchange = { path = "../interchange", features = ["ar2"] } \ No newline at end of file +wxbox-interchange = { path = "../interchange", features = ["ar2"] }