Compare commits
2 commits
2a5052f47f
...
3c94d3a32a
| Author | SHA1 | Date | |
|---|---|---|---|
| 3c94d3a32a | |||
| cee9d47dc2 |
3 changed files with 360 additions and 382 deletions
|
|
@ -2,18 +2,22 @@
|
|||
import maplibregl from 'maplibre-gl';
|
||||
import ToolbarProductSelector from '$lib/ToolbarProductSelector.svelte';
|
||||
import { Button } from '$lib/components/ui/button';
|
||||
import {mount, onMount} from 'svelte';
|
||||
import { mount, onMount } from 'svelte';
|
||||
import type { LayerList } from './layerList';
|
||||
import type { StationStatus } from './stationData';
|
||||
import {toast} from "svelte-sonner";
|
||||
import {Toggle} from "$lib/components/ui/toggle";
|
||||
import CloudAlertIcon from "@lucide/svelte/icons/cloud-alert";
|
||||
import {borderLUT, fillLUT} from "$lib/alertLayer";
|
||||
import AlertPopup from "$lib/AlertPopup.svelte";
|
||||
import {CifContainer} from "$lib/generated_interop/cifContainer";
|
||||
import fragmentSource from "$lib/map/fragment.glsl?raw";
|
||||
import vertexSource from "$lib/map/vertex.glsl?raw";
|
||||
import {degreesToRadians} from "@turf/turf";
|
||||
import { toast } from 'svelte-sonner';
|
||||
import { Toggle } from '$lib/components/ui/toggle';
|
||||
import CloudAlertIcon from '@lucide/svelte/icons/cloud-alert';
|
||||
import { borderLUT, fillLUT } from '$lib/alertLayer';
|
||||
import AlertPopup from '$lib/AlertPopup.svelte';
|
||||
import { CifContainer } from '$lib/generated_interop/cifContainer';
|
||||
import fragmentSource from '$lib/map/fragment.glsl?raw';
|
||||
import vertexSource from '$lib/map/vertex.glsl?raw';
|
||||
import { degreesToRadians } from '@turf/turf';
|
||||
import type { DigitalRadarData, Radial } from './generated_interop/digitalRadarData';
|
||||
|
||||
const BELOW_THRESHOLD = -9999.0;
|
||||
const RANGE_FOLDED = -9998.0;
|
||||
|
||||
interface Props {
|
||||
categories: LayerList;
|
||||
|
|
@ -35,8 +39,8 @@
|
|||
return;
|
||||
}
|
||||
|
||||
map.setLayoutProperty('alerts', 'visibility', showAlertLayer ? "visible" : "none");
|
||||
map.setLayoutProperty('alerts-outline', 'visibility', showAlertLayer ? "visible" : "none");
|
||||
map.setLayoutProperty('alerts', 'visibility', showAlertLayer ? 'visible' : 'none');
|
||||
map.setLayoutProperty('alerts-outline', 'visibility', showAlertLayer ? 'visible' : 'none');
|
||||
}
|
||||
|
||||
onMount(() => {
|
||||
|
|
@ -70,42 +74,47 @@
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
map.addSource('alerts', {
|
||||
type: 'geojson',
|
||||
data: 'https://api.weather.gov/alerts/active?status=actual'
|
||||
});
|
||||
map.addLayer({
|
||||
id: 'alerts',
|
||||
type: 'fill',
|
||||
source: 'alerts',
|
||||
paint: {
|
||||
// @ts-expect-error this type is too complicated
|
||||
'fill-color': fillLUT,
|
||||
}
|
||||
}, sym_layer);
|
||||
map.addLayer({
|
||||
id: 'alerts-outline',
|
||||
type: 'line',
|
||||
source: 'alerts',
|
||||
layout: {
|
||||
'line-join': 'round'
|
||||
map.addLayer(
|
||||
{
|
||||
id: 'alerts',
|
||||
type: 'fill',
|
||||
source: 'alerts',
|
||||
paint: {
|
||||
// @ts-expect-error this type is too complicated
|
||||
'fill-color': fillLUT
|
||||
}
|
||||
},
|
||||
paint: {
|
||||
// @ts-expect-error this type is too complicated
|
||||
'line-color': borderLUT,
|
||||
'line-width': 3
|
||||
}
|
||||
}, sym_layer);
|
||||
sym_layer
|
||||
);
|
||||
map.addLayer(
|
||||
{
|
||||
id: 'alerts-outline',
|
||||
type: 'line',
|
||||
source: 'alerts',
|
||||
layout: {
|
||||
'line-join': 'round'
|
||||
},
|
||||
paint: {
|
||||
// @ts-expect-error this type is too complicated
|
||||
'line-color': borderLUT,
|
||||
'line-width': 3
|
||||
}
|
||||
},
|
||||
sym_layer
|
||||
);
|
||||
let createPopup = (e, full: boolean) => {
|
||||
const randLetter = String.fromCharCode(65 + Math.floor(Math.random() * 26));
|
||||
const uniqid = randLetter + Date.now();
|
||||
|
||||
new maplibregl.Popup({ className: 'popup' })
|
||||
.setLngLat(e.lngLat)
|
||||
.setHTML(`<div id="${uniqid}"></div>`)
|
||||
.setLngLat(e.lngLat)
|
||||
.setHTML(`<div id="${uniqid}"></div>`)
|
||||
|
||||
.addTo(map);
|
||||
.addTo(map);
|
||||
mount(AlertPopup, {
|
||||
target: document.getElementById(uniqid)!,
|
||||
props: {
|
||||
|
|
@ -114,30 +123,31 @@
|
|||
}
|
||||
});
|
||||
};
|
||||
map.on('click', 'alerts', (e) => { createPopup(e, true) });
|
||||
map.on('click', 'alerts', (e) => {
|
||||
createPopup(e, true);
|
||||
});
|
||||
let currentFeatureCoordinates = undefined;
|
||||
const popup = new maplibregl.Popup({
|
||||
closeButton: false,
|
||||
closeOnClick: false
|
||||
});
|
||||
map.on('mousemove', 'alerts', (e) => {
|
||||
// Change the cursor style as a UI indicator.
|
||||
map.getCanvas().style.cursor = 'pointer';
|
||||
// Change the cursor style as a UI indicator.
|
||||
map.getCanvas().style.cursor = 'pointer';
|
||||
|
||||
// Populate the popup and set its coordinates
|
||||
// based on the feature found.
|
||||
const randLetter = String.fromCharCode(65 + Math.floor(Math.random() * 26));
|
||||
const uniqid = randLetter + Date.now();
|
||||
|
||||
popup.setLngLat(e.lngLat).setHTML(`<div id="${uniqid}"></div>`).addTo(map);
|
||||
mount(AlertPopup, {
|
||||
target: document.getElementById(uniqid)!,
|
||||
props: {
|
||||
showFull: false,
|
||||
...e.features[0].properties
|
||||
}
|
||||
});
|
||||
// Populate the popup and set its coordinates
|
||||
// based on the feature found.
|
||||
const randLetter = String.fromCharCode(65 + Math.floor(Math.random() * 26));
|
||||
const uniqid = randLetter + Date.now();
|
||||
|
||||
popup.setLngLat(e.lngLat).setHTML(`<div id="${uniqid}"></div>`).addTo(map);
|
||||
mount(AlertPopup, {
|
||||
target: document.getElementById(uniqid)!,
|
||||
props: {
|
||||
showFull: false,
|
||||
...e.features[0].properties
|
||||
}
|
||||
});
|
||||
});
|
||||
map.on('mouseenter', 'alerts', () => {
|
||||
map.getCanvas().style.cursor = 'pointer';
|
||||
|
|
@ -147,110 +157,163 @@
|
|||
popup.remove();
|
||||
});
|
||||
|
||||
let r = await fetch('http://localhost:3000/v2/nexrad/l2/KFSD/1/REF');
|
||||
let r = await fetch('http://localhost:3000/v2/nexrad/l2/KLZK/1/REF');
|
||||
const buf = await r.arrayBuffer();
|
||||
|
||||
const container = CifContainer.fromBinary(new Uint8Array(buf));
|
||||
console.log('wxrad http://localhost:3000/v2/nexrad/l2/KFSD/1/REF: ', container);
|
||||
console.log('wxrad http://localhost:3000/v2/nexrad/l2/KLZK/1/REF: ', container);
|
||||
if (container.messageType.oneofKind == 'digitalRadarData') {
|
||||
const drd: DigitalRadarData = container.messageType.digitalRadarData;
|
||||
|
||||
const dataLayer = {
|
||||
id: 'dataGl',
|
||||
type: 'custom',
|
||||
onAdd(map, gl) {
|
||||
const vertexShader = gl.createShader(gl.VERTEX_SHADER);
|
||||
gl.shaderSource(vertexShader, vertexSource);
|
||||
gl.compileShader(vertexShader);
|
||||
|
||||
const dataLayer = {
|
||||
id: 'dataGl',
|
||||
type: 'custom',
|
||||
onAdd(map, gl) {
|
||||
const vertexShader = gl.createShader(gl.VERTEX_SHADER);
|
||||
gl.shaderSource(vertexShader, vertexSource);
|
||||
gl.compileShader(vertexShader);
|
||||
const fragmentShader = gl.createShader(gl.FRAGMENT_SHADER);
|
||||
gl.shaderSource(fragmentShader, fragmentSource);
|
||||
gl.compileShader(fragmentShader);
|
||||
|
||||
const fragmentShader = gl.createShader(gl.FRAGMENT_SHADER);
|
||||
gl.shaderSource(fragmentShader, fragmentSource);
|
||||
gl.compileShader(fragmentShader);
|
||||
this.program = gl.createProgram();
|
||||
gl.attachShader(this.program, vertexShader);
|
||||
gl.attachShader(this.program, fragmentShader);
|
||||
gl.linkProgram(this.program);
|
||||
|
||||
this.program = gl.createProgram();
|
||||
gl.attachShader(this.program, vertexShader);
|
||||
gl.attachShader(this.program, fragmentShader);
|
||||
gl.linkProgram(this.program);
|
||||
this.aPos = gl.getAttribLocation(this.program, 'a_pos');
|
||||
|
||||
this.aPos = gl.getAttribLocation(this.program, 'a_pos');
|
||||
|
||||
/*
|
||||
/*
|
||||
110.574 km = 1 deg
|
||||
deg = 1/110.574 km
|
||||
|
||||
111.320 * cos(latitude) km = 1 deg
|
||||
|
||||
*/
|
||||
const lat = 43.58777778;
|
||||
const long = -96.72944444;
|
||||
gl.useProgram(this.program);
|
||||
const lat = 34.8365;
|
||||
gl.uniform1f(gl.getUniformLocation(this.program, 'radarLat'), lat);
|
||||
const long = -92.262194;
|
||||
gl.uniform1f(gl.getUniformLocation(this.program, 'radarLng'), long);
|
||||
|
||||
const radar_range_maximum = 560; // ish km
|
||||
const degrees_of_lat = radar_range_maximum / 110.574;
|
||||
const degrees_of_long = radar_range_maximum / (111.320 * Math.cos(degreesToRadians(lat)));
|
||||
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
|
||||
});
|
||||
// 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(
|
||||
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([a.x, a.y, b.x, b.y, c.x, c.y, b.x, b.y, d.x, d.y, c.x, c.y]),
|
||||
gl.STATIC_DRAW
|
||||
);
|
||||
},
|
||||
render(gl, args) {
|
||||
gl.useProgram(this.program);
|
||||
gl.uniformMatrix4fv(
|
||||
);
|
||||
|
||||
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
|
||||
);
|
||||
gl.uniform1f(
|
||||
gl.getUniformLocation(this.program, 'sampleInterval'),
|
||||
drd.radials[0].product?.data?.sampleInterval
|
||||
);
|
||||
|
||||
function scaleMomentData(radial: Radial, u: number): number {
|
||||
if (!radial.product || !radial.product.data) return BELOW_THRESHOLD;
|
||||
if (radial.product?.data?.scale == 0) {
|
||||
return u;
|
||||
} else {
|
||||
if (u == 0) {
|
||||
return BELOW_THRESHOLD;
|
||||
} else if (u == 1) {
|
||||
return RANGE_FOLDED;
|
||||
} else {
|
||||
return (u - radial.product.data?.offset) / radial.product?.data?.scale;
|
||||
}
|
||||
}
|
||||
}
|
||||
const rdata: number[] = [];
|
||||
for (const radial of drd.radials) {
|
||||
for (const tdata of radial.product.data.data) {
|
||||
rdata.push(scaleMomentData(radial, tdata));
|
||||
}
|
||||
}
|
||||
const data = new Float32Array(rdata);
|
||||
|
||||
gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 4);
|
||||
|
||||
gl.activeTexture(gl.TEXTURE0 + 4);
|
||||
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.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);
|
||||
},
|
||||
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,
|
||||
args.defaultProjectionData.mainMatrix
|
||||
);
|
||||
gl.bindBuffer(gl.ARRAY_BUFFER, this.buffer);
|
||||
gl.enableVertexAttribArray(this.aPos);
|
||||
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.bindBuffer(gl.ARRAY_BUFFER, this.buffer);
|
||||
gl.enableVertexAttribArray(this.aPos);
|
||||
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);
|
||||
}
|
||||
};
|
||||
map.addLayer(dataLayer);
|
||||
}
|
||||
map.addLayer(dataLayer);
|
||||
});
|
||||
map.on('error', (e) => {
|
||||
console.error(e);
|
||||
toast.error('Data loading failed!');
|
||||
});
|
||||
|
||||
|
||||
});
|
||||
</script>
|
||||
|
||||
|
|
@ -274,20 +337,17 @@
|
|||
<!-- center-center -->
|
||||
<div class="flex flex-1 grow flex-row justify-center">
|
||||
<ToolbarProductSelector
|
||||
bind:pickingSiteForCategory
|
||||
bind:selectedSite
|
||||
bind:selectedPrimaryLayer
|
||||
{categories}
|
||||
{map}
|
||||
{stations}
|
||||
bind:pickingSiteForCategory
|
||||
bind:selectedSite
|
||||
bind:selectedPrimaryLayer
|
||||
{categories}
|
||||
{map}
|
||||
{stations}
|
||||
/>
|
||||
</div>
|
||||
|
||||
<!-- center-center -->
|
||||
<div class="flex flex-1 grow flex-row justify-start">
|
||||
|
||||
</div>
|
||||
|
||||
<div class="flex flex-1 grow flex-row justify-start"></div>
|
||||
</div>
|
||||
|
||||
<div class="flex flex-1 grow flex-row justify-end gap-1"></div>
|
||||
|
|
|
|||
|
|
@ -19,242 +19,136 @@ 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;
|
||||
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;
|
||||
|
||||
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;
|
||||
struct LocateRadialResult {
|
||||
int radialIndex;
|
||||
float radialDistance;
|
||||
bool didFindRadial;
|
||||
};
|
||||
|
||||
const float tolerance = pow(10.0, -13.0);
|
||||
LocateRadialResult locateRadial(float forAzimuth) {
|
||||
int closestRadial;
|
||||
float bestDistance = 1.0;
|
||||
bool foundAnything = false;
|
||||
|
||||
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;
|
||||
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 {
|
||||
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);
|
||||
closestRadial = i;
|
||||
bestDistance = this_dist;
|
||||
foundAnything = true;
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
return LocateRadialResult(closestRadial, bestDistance, foundAnything);
|
||||
}
|
||||
|
||||
void main() {
|
||||
Ellipsoid WGS84 = fromAAndInverseF(6378137.0, 298.257223563);
|
||||
float lat;
|
||||
float lng;
|
||||
xyToLngLat(raw_pos.x, raw_pos.y, lat, lng);
|
||||
|
||||
float radarLat = 43.58777778;
|
||||
float radarLng = -96.72944444;
|
||||
if (azimuthCount == 0) {
|
||||
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
|
||||
return;
|
||||
}
|
||||
|
||||
float a1;
|
||||
float a2;
|
||||
float s;
|
||||
float R = 6371.0 * pow(10.0, 3.0); // meters
|
||||
float phi1 = radians(radarLat);
|
||||
float phi2 = radians(lat);
|
||||
float lambda1 = radians(radarLng);
|
||||
float lambda2 = radians(lng);
|
||||
|
||||
// should be 43.58777778
|
||||
// 88 > x > 89
|
||||
// >40.0
|
||||
float delta = 0.01;
|
||||
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
|
||||
|
||||
GlobalCoordinates radar = GlobalCoordinates(degreeAngle(radarLat), degreeAngle(radarLng));
|
||||
canonicalizeGlobalCoordinates(radar);
|
||||
GlobalCoordinates
|
||||
samplePoint = GlobalCoordinates(degreeAngle(lat), degreeAngle(lng));
|
||||
canonicalizeGlobalCoordinates(samplePoint);
|
||||
float d_m = d;
|
||||
|
||||
Geodedic vincentyResult = vincenty(WGS84, radar, samplePoint);
|
||||
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 = (theta * 180.0 / PI + 360.0); // degrees
|
||||
|
||||
//
|
||||
// abs(lat-43.5877) <delta && abs(lng+96.72944444)<delta
|
||||
if (vincentyResult.didGetGoodEstimate) {
|
||||
fragColor = vec4(0.0, 1.0, 0.0, 0.5);
|
||||
LocateRadialResult maybeRadial = locateRadial(azimuth);
|
||||
if (!maybeRadial.didFindRadial) {
|
||||
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
|
||||
return;
|
||||
}
|
||||
|
||||
float distance_km = d_m / 1000.0;
|
||||
float first_gate_distance_km = startRange / 1000.0;
|
||||
float gate_spacing_km = sampleInterval / 1000.0;
|
||||
|
||||
if (distance_km < first_gate_distance_km) {
|
||||
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
|
||||
return;
|
||||
}
|
||||
|
||||
int gate = int(round((distance_km - first_gate_distance_km) / gate_spacing_km));
|
||||
if (gate > 1832) {
|
||||
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
|
||||
return;
|
||||
}
|
||||
|
||||
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(1.0, 0.0, 0.0, 0.5);
|
||||
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,14 +1,13 @@
|
|||
use crate::cif::cif_container::MessageType;
|
||||
use crate::{AsCif, cif};
|
||||
use nexrad_decode::messages::digital_radar_data::RadialStatus;
|
||||
use wxbox_ar2::{MomentData, Radial, Scan};
|
||||
use crate::{cif, AsCif};
|
||||
use crate::cif::cif_container::MessageType;
|
||||
|
||||
pub struct Ar2AsCifParams {
|
||||
pub requested_elevation: u8,
|
||||
pub requested_product: String
|
||||
pub requested_product: String,
|
||||
}
|
||||
|
||||
|
||||
impl AsCif for Scan {
|
||||
type Params = Ar2AsCifParams;
|
||||
fn as_cif(&self, params: &Self::Params) -> cif::CifContainer {
|
||||
|
|
@ -17,7 +16,10 @@ impl AsCif for Scan {
|
|||
digital_radar_data.vcp_number = self.coverage_pattern_number as i32;
|
||||
|
||||
// find the elevation
|
||||
let maybe_elevation = self.sweeps.iter().find(|u| u.elevation_number == params.requested_elevation);
|
||||
let 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![];
|
||||
|
|
@ -25,30 +27,38 @@ impl AsCif for Scan {
|
|||
match params.requested_product.as_str() {
|
||||
"REF" if radial.reflectivity.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.reflectivity.as_ref().unwrap()));
|
||||
},
|
||||
}
|
||||
"VEL" if radial.velocity.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.velocity.as_ref().unwrap()));
|
||||
},
|
||||
}
|
||||
"SW" if radial.spectrum_width.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.spectrum_width.as_ref().unwrap()));
|
||||
},
|
||||
radials_we_can_use
|
||||
.push((&radial, &radial.spectrum_width.as_ref().unwrap()));
|
||||
}
|
||||
"ZDR" if radial.differential_reflectivity.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.differential_reflectivity.as_ref().unwrap()));
|
||||
},
|
||||
radials_we_can_use
|
||||
.push((&radial, &radial.differential_reflectivity.as_ref().unwrap()));
|
||||
}
|
||||
"PHI" if radial.differential_phase.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.differential_phase.as_ref().unwrap()));
|
||||
},
|
||||
radials_we_can_use
|
||||
.push((&radial, &radial.differential_phase.as_ref().unwrap()));
|
||||
}
|
||||
"RHO" if radial.correlation_coefficient.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.correlation_coefficient.as_ref().unwrap()));
|
||||
},
|
||||
radials_we_can_use
|
||||
.push((&radial, &radial.correlation_coefficient.as_ref().unwrap()));
|
||||
}
|
||||
"CFP" if radial.specific_differential_phase.is_some() => {
|
||||
radials_we_can_use.push((&radial, &radial.specific_differential_phase.as_ref().unwrap()));
|
||||
},
|
||||
radials_we_can_use.push((
|
||||
&radial,
|
||||
&radial.specific_differential_phase.as_ref().unwrap(),
|
||||
));
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
|
||||
digital_radar_data.radials = radials_we_can_use.iter()
|
||||
digital_radar_data.radials = radials_we_can_use
|
||||
.iter()
|
||||
.map(|u| {
|
||||
let mut radial = cif::digital_radar_data::Radial::default();
|
||||
radial.elevation_number = u.0.elevation_number as i32;
|
||||
|
|
@ -58,13 +68,26 @@ impl AsCif for Scan {
|
|||
radial.collection_timestamp = u.0.collection_timestamp as u64;
|
||||
radial.elevation_degrees = u.0.elevation_number_degrees;
|
||||
radial.radial_status = match u.0.radial_status {
|
||||
RadialStatus::ElevationStart => cif::digital_radar_data::RadialStatus::ElevationStart,
|
||||
RadialStatus::IntermediateRadialData => cif::digital_radar_data::RadialStatus::IntermediateRadialData,
|
||||
RadialStatus::ElevationEnd => cif::digital_radar_data::RadialStatus::ElevationEnd,
|
||||
RadialStatus::VolumeScanStart => cif::digital_radar_data::RadialStatus::VolumeScanStart,
|
||||
RadialStatus::VolumeScanEnd => cif::digital_radar_data::RadialStatus::VolumeScanEnd,
|
||||
RadialStatus::ElevationStartVCPFinal => cif::digital_radar_data::RadialStatus::ElevationStartVcpFinal,
|
||||
}.into();
|
||||
RadialStatus::ElevationStart => {
|
||||
cif::digital_radar_data::RadialStatus::ElevationStart
|
||||
}
|
||||
RadialStatus::IntermediateRadialData => {
|
||||
cif::digital_radar_data::RadialStatus::IntermediateRadialData
|
||||
}
|
||||
RadialStatus::ElevationEnd => {
|
||||
cif::digital_radar_data::RadialStatus::ElevationEnd
|
||||
}
|
||||
RadialStatus::VolumeScanStart => {
|
||||
cif::digital_radar_data::RadialStatus::VolumeScanStart
|
||||
}
|
||||
RadialStatus::VolumeScanEnd => {
|
||||
cif::digital_radar_data::RadialStatus::VolumeScanEnd
|
||||
}
|
||||
RadialStatus::ElevationStartVCPFinal => {
|
||||
cif::digital_radar_data::RadialStatus::ElevationStartVcpFinal
|
||||
}
|
||||
}
|
||||
.into();
|
||||
|
||||
let mut moment_meta = cif::digital_radar_data::MomentaryMeta::default();
|
||||
moment_meta.product_name = params.requested_product.clone();
|
||||
|
|
@ -78,9 +101,10 @@ impl AsCif for Scan {
|
|||
|
||||
let mut moment_product = cif::digital_radar_data::MomentaryProduct::default();
|
||||
moment_product.product_metadata = Some(moment_meta);
|
||||
moment_product.data = Some(moment_data);
|
||||
|
||||
radial.product = Some(moment_product);
|
||||
|
||||
|
||||
radial
|
||||
})
|
||||
.collect()
|
||||
|
|
@ -90,4 +114,4 @@ impl AsCif for Scan {
|
|||
container.message_type = Some(MessageType::DigitalRadarData(digital_radar_data));
|
||||
container
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue