diff --git a/src/lib/data/ZarrDataManager.ts b/src/lib/data/ZarrDataManager.ts index 1dbc5df..2dcd833 100644 --- a/src/lib/data/ZarrDataManager.ts +++ b/src/lib/data/ZarrDataManager.ts @@ -1,6 +1,8 @@ import QuickLRU from "quick-lru"; import * as zarr from "zarrita"; +import "./codecs.ts"; + import { createIcechunkStore, isIcechunkStorePath, diff --git a/src/lib/data/codecs.ts b/src/lib/data/codecs.ts index 41621ce..62d8a8e 100644 --- a/src/lib/data/codecs.ts +++ b/src/lib/data/codecs.ts @@ -2,5 +2,7 @@ import { registerPCodec } from "@eeholmes/zarrita-pcodec"; import { registry } from "zarrita"; import { Fletcher32Codec } from "./fletcher32.ts"; +import { GribscanRawGribCodec } from "./gribscan.ts"; registry.set("numcodecs.fletcher32", async () => Fletcher32Codec); +registry.set("numcodecs.gribscan.rawgrib", async () => GribscanRawGribCodec); registerPCodec(registry); diff --git a/src/lib/data/gribscan.ts b/src/lib/data/gribscan.ts new file mode 100644 index 0000000..7a07197 --- /dev/null +++ b/src/lib/data/gribscan.ts @@ -0,0 +1,1168 @@ +import type { DataType } from "zarrita"; + +const Endian = { + BIG: "big", + LITTLE: "little", +} as const; + +type TEndian = (typeof Endian)[keyof typeof Endian]; + +const GribEdition = { + ONE: 1, + TWO: 2, +} as const; + +const Grib2DataRepresentationTemplate = { + COMPLEX_PACKING: 2, + COMPLEX_PACKING_SPATIAL_DIFFERENCING: 3, + JPEG2000: 40, + PNG: 41, + CCSDS: 42, + SIMPLE_PACKING: 0, +} as const; + +const AecFlag = { + DATA_3BYTE: 1 << 1, + DATA_PREPROCESS: 1 << 3, + DATA_SIGNED: 1 << 0, + MSB: 1 << 2, + PAD_RSI: 1 << 5, + RESTRICTED: 1 << 4, +} as const; + +type TCodecMeta = { + codecs?: Array<{ + configuration?: Record; + name: string; + }>; + dataType?: DataType; +}; + +type TSection = { + length: number; + number: number; + offset: number; +}; + +type TGrib1Bitmap = { + count: number; + data: Uint8Array; +}; + +type TGrib1BinaryData = { + bitsPerValue: number; + dataBytes: Uint8Array; + packedCount: number; + referenceValue: number; + scale: number; +}; + +type TAecParams = { + bitsPerSample: number; + blockSize: number; + flags: number; + rsi: number; +}; + +type TDecodedArray = + | BigInt64Array + | BigUint64Array + | Float32Array + | Float64Array + | Int8Array + | Int16Array + | Int32Array + | Uint8Array + | Uint16Array + | Uint32Array; + +export class GribscanRawGribCodec { + readonly kind = "bytes_to_bytes"; + readonly #dataType: DataType; + readonly #endian: TEndian; + + constructor(dataType: DataType = "float64", endian: TEndian = Endian.LITTLE) { + this.#dataType = dataType; + this.#endian = endian; + } + + static fromConfig(_configuration: unknown, meta?: TCodecMeta) { + return new GribscanRawGribCodec( + meta?.dataType ?? "float64", + getBytesEndian(meta) + ); + } + + encode(arr: Uint8Array): Uint8Array { + return arr; + } + + decode(arr: Uint8Array): Uint8Array { + return valuesToBytes(decodeGribMessage(arr), this.#dataType, this.#endian); + } +} + +class BitReader { + #bytes: Uint8Array; + #position: number; + + constructor(bytes: Uint8Array, byteOffset = 0) { + this.#bytes = bytes; + this.#position = byteOffset * 8; + } + + alignToByte() { + const remainder = this.#position % 8; + if (remainder !== 0) { + this.#position += 8 - remainder; + } + } + + readUint(bits: number) { + if (bits === 0) { + return 0; + } + + let result = 0; + let remaining = bits; + while (remaining > 0) { + const byteIndex = this.#position >> 3; + if (byteIndex >= this.#bytes.byteLength) { + throw new Error("Unexpected end of GRIB bitstream"); + } + const bitOffset = this.#position & 7; + const available = 8 - bitOffset; + const take = Math.min(available, remaining); + const shift = available - take; + const mask = (1 << take) - 1; + + result = (result << take) | ((this.#bytes[byteIndex] >> shift) & mask); + this.#position += take; + remaining -= take; + } + + return result >>> 0; + } + + readBit() { + return this.readUint(1) === 1; + } +} + +function decodeGribMessage(buffer: ArrayBufferLike | Uint8Array) { + const buf = buffer instanceof Uint8Array ? buffer : new Uint8Array(buffer); + + if ( + buf[0] !== 0x47 || + buf[1] !== 0x52 || + buf[2] !== 0x49 || + buf[3] !== 0x42 + ) { + throw new Error("Not a GRIB message: missing 'GRIB' magic bytes"); + } + + const edition = buf[7]; + if (edition === GribEdition.ONE) { + return decodeGrib1(buf); + } + if (edition === GribEdition.TWO) { + return decodeGrib2(buf); + } + throw new Error(`Unsupported GRIB edition: ${edition}`); +} + +function decodeGrib1(buf: Uint8Array) { + const view = createDataView(buf); + const pdsStart = 8; + const pdsLength = readUint24(view, pdsStart); + const flagByte = view.getUint8(pdsStart + 7); + const gdsPresent = (flagByte & 0x80) !== 0; + const bmsPresent = (flagByte & 0x40) !== 0; + const decimalScaleFactor = readGrib1DecimalScaleFactor(view, pdsStart); + + let offset = pdsStart + pdsLength; + + if (gdsPresent) { + offset += readUint24(view, offset); + } + + const { bitmap, nextOffset } = readGrib1Bitmap(buf, view, offset, bmsPresent); + const binaryData = readGrib1BinaryData( + buf, + view, + nextOffset, + decimalScaleFactor + ); + const packed = unpackGrib1Values(binaryData); + + return bitmap ? applyGrib1Bitmap(bitmap, packed) : packed; +} + +function readGrib1Bitmap( + buf: Uint8Array, + view: DataView, + offset: number, + bmsPresent: boolean +) { + if (!bmsPresent) { + return { bitmap: undefined, nextOffset: offset }; + } + + const bmsLength = readUint24(view, offset); + const unusedBits = view.getUint8(offset + 3); + const bitmapByteCount = bmsLength - 6; + return { + bitmap: { + count: bitmapByteCount * 8 - unusedBits, + data: new Uint8Array( + buf.buffer, + buf.byteOffset + offset + 6, + bitmapByteCount + ), + }, + nextOffset: offset + bmsLength, + }; +} + +function readGrib1BinaryData( + buf: Uint8Array, + view: DataView, + bdsStart: number, + decimalScaleFactor: number +): TGrib1BinaryData { + const bdsLength = readUint24(view, bdsStart); + const flagBds = view.getUint8(bdsStart + 3); + const unusedBitsAtEnd = flagBds & 0x0f; + const binaryScaleFactor = signMagnitude( + view.getUint16(bdsStart + 4, false), + 16 + ); + const referenceValue = readIbmFloat32(view, bdsStart + 6); + const bitsPerValue = view.getUint8(bdsStart + 10); + const dataStart = bdsStart + 11; + const dataBytes = new Uint8Array( + buf.buffer, + buf.byteOffset + dataStart, + bdsLength - 11 + ); + const scale = + Math.pow(2, binaryScaleFactor) / Math.pow(10, decimalScaleFactor); + const packedCount = + bitsPerValue > 0 + ? Math.floor((dataBytes.length * 8 - unusedBitsAtEnd) / bitsPerValue) + : 0; + + return { bitsPerValue, dataBytes, packedCount, referenceValue, scale }; +} + +function unpackGrib1Values(binaryData: TGrib1BinaryData) { + const packed = new Float64Array(binaryData.packedCount); + const reader = new BitReader(binaryData.dataBytes); + for (let i = 0; i < binaryData.packedCount; i++) { + packed[i] = + binaryData.referenceValue + + reader.readUint(binaryData.bitsPerValue) * binaryData.scale; + } + return packed; +} + +function applyGrib1Bitmap(bitmap: TGrib1Bitmap, packed: Float64Array) { + const out = new Float64Array(bitmap.count).fill(NaN); + const bitmapReader = new BitReader(bitmap.data); + let packedIndex = 0; + for (let i = 0; i < bitmap.count; i++) { + if (bitmapReader.readUint(1) === 1) { + out[i] = packed[packedIndex++]; + } + } + return out; +} + +function readGrib1DecimalScaleFactor(view: DataView, pdsStart: number) { + return signMagnitude(view.getUint16(pdsStart + 26, false), 16); +} + +function decodeGrib2(buf: Uint8Array) { + const { sections, view } = parseGrib2Sections(buf); + const sectionsByNumber = new Map( + sections.map((section) => [section.number, section]) + ); + const section3 = sectionsByNumber.get(3); + const section5 = sectionsByNumber.get(5); + const section6 = sectionsByNumber.get(6); + const section7 = sectionsByNumber.get(7); + + if (!section3 || !section5 || !section6 || !section7) { + throw new Error("GRIB2 message missing required sections (3, 5, 6, or 7)"); + } + + const totalPoints = view.getUint32(section3.offset + 6, false); + const template = view.getUint16(section5.offset + 9, false); + switch (template) { + case Grib2DataRepresentationTemplate.SIMPLE_PACKING: + return decodeTemplate0(section5, section6, section7, view, totalPoints); + case Grib2DataRepresentationTemplate.COMPLEX_PACKING: + return decodeTemplate2(section5, section6, section7, view, totalPoints); + case Grib2DataRepresentationTemplate.COMPLEX_PACKING_SPATIAL_DIFFERENCING: + return decodeTemplate3(section5, section6, section7, view, totalPoints); + case Grib2DataRepresentationTemplate.JPEG2000: + throw new Error("GRIB2 template 40 (JPEG 2000) is not supported"); + case Grib2DataRepresentationTemplate.PNG: + throw new Error("GRIB2 template 41 (PNG) is not supported"); + case Grib2DataRepresentationTemplate.CCSDS: + return decodeTemplate42(section5, section6, section7, view, totalPoints); + default: + throw new Error( + `GRIB2 data representation template ${template} is not supported` + ); + } +} + +function decodeTemplate0( + section5: TSection, + section6: TSection, + section7: TSection, + view: DataView, + totalPoints: number +) { + const s5 = section5.offset; + const referenceValue = view.getFloat32(s5 + 11, false); + const binaryScaleFactor = view.getInt16(s5 + 15, false); + const decimalScaleFactor = view.getInt16(s5 + 17, false); + const bitsPerValue = view.getUint8(s5 + 19); + const valueCount = view.getUint32(s5 + 5, false); + const scale = + Math.pow(2, binaryScaleFactor) / Math.pow(10, decimalScaleFactor); + + return unpackSimple( + section6, + section7, + view, + referenceValue, + scale, + bitsPerValue, + valueCount, + totalPoints + ); +} + +function decodeTemplate2( + section5: TSection, + section6: TSection, + section7: TSection, + view: DataView, + totalPoints: number +) { + const params = readComplexPackingParams(section5, view); + const reader = createSection7Reader(section7, view); + const groupReference = readGroupReferences(reader, params); + reader.alignToByte(); + const groupWidth = readGroupWidths(reader, params); + reader.alignToByte(); + const groupLength = readGroupLengths(reader, params); + reader.alignToByte(); + + const out = new Float64Array(params.valueCount); + let outIndex = 0; + for (let group = 0; group < params.groupCount; group++) { + const width = groupWidth[group]; + const length = groupLength[group]; + for (let i = 0; i < length; i++) { + const value = width > 0 ? reader.readUint(width) : 0; + if ( + isComplexPackingMissingValue(value, width, params.missingManagement) + ) { + out[outIndex++] = NaN; + continue; + } + out[outIndex++] = + params.referenceValue + (groupReference[group] + value) * params.scale; + } + } + + return applyBitmap(section6, view, out, totalPoints); +} + +function decodeTemplate3( + section5: TSection, + section6: TSection, + section7: TSection, + view: DataView, + totalPoints: number +) { + const params = readComplexPackingParams(section5, view); + const orderSpatialDifferencing = view.getUint8(section5.offset + 47); + const extraDescriptorOctets = view.getUint8(section5.offset + 48); + const extraDescriptorBits = extraDescriptorOctets * 8; + const reader = createSection7Reader(section7, view); + const initialValues = Array.from({ length: orderSpatialDifferencing }, () => + signMagnitude(reader.readUint(extraDescriptorBits), extraDescriptorBits) + ); + const overallMinimum = signMagnitude( + reader.readUint(extraDescriptorBits), + extraDescriptorBits + ); + const groupReference = readGroupReferences(reader, params); + reader.alignToByte(); + const groupWidth = readGroupWidths(reader, params); + reader.alignToByte(); + const groupLength = readGroupLengths(reader, params); + reader.alignToByte(); + const values = unpackComplexIntegerValues( + reader, + params, + groupReference, + groupWidth, + groupLength + ); + + restoreSpatialDifferences( + values, + initialValues, + overallMinimum, + orderSpatialDifferencing + ); + + const out = new Float64Array(params.valueCount); + for (let i = 0; i < params.valueCount; i++) { + out[i] = params.referenceValue + values[i] * params.scale; + } + return applyBitmap(section6, view, out, totalPoints); +} + +function decodeTemplate42( + section5: TSection, + section6: TSection, + section7: TSection, + view: DataView, + totalPoints: number +) { + const s5 = section5.offset; + const binaryScaleFactor = view.getInt16(s5 + 15, false); + const decimalScaleFactor = view.getInt16(s5 + 17, false); + const params = { + bitsPerSample: view.getUint8(s5 + 19), + blockSize: view.getUint8(s5 + 22), + flags: aecFlagsFromGrib2(view.getUint8(s5 + 21)), + rsi: view.getUint16(s5 + 23, false), + }; + const sampleValues = decodeAecSamples( + createSection7Bytes(section7, view), + params, + view.getUint32(s5 + 5, false) + ); + const scale = + Math.pow(2, binaryScaleFactor) / Math.pow(10, decimalScaleFactor); + + return applyBitmap( + section6, + view, + scaleAecSamples(sampleValues, view.getFloat32(s5 + 11, false), scale), + totalPoints + ); +} + +function readComplexPackingParams(section5: TSection, view: DataView) { + const s5 = section5.offset; + const binaryScaleFactor = view.getInt16(s5 + 15, false); + const decimalScaleFactor = view.getInt16(s5 + 17, false); + return { + bitsGroupWidth: view.getUint8(s5 + 36), + bitsPerValue: view.getUint8(s5 + 19), + bitsScaledLength: view.getUint8(s5 + 46), + groupCount: view.getUint32(s5 + 31, false), + lastGroupLength: view.getUint32(s5 + 42, false), + lengthIncrement: view.getUint8(s5 + 41), + missingManagement: view.getUint8(s5 + 22), + referenceGroupLength: view.getUint32(s5 + 37, false), + referenceGroupWidth: view.getUint8(s5 + 35), + referenceValue: view.getFloat32(s5 + 11, false), + scale: Math.pow(2, binaryScaleFactor) / Math.pow(10, decimalScaleFactor), + valueCount: view.getUint32(s5 + 5, false), + }; +} + +function readGroupReferences( + reader: BitReader, + params: ReturnType +) { + const groupReference = new Uint32Array(params.groupCount); + for (let group = 0; group < params.groupCount; group++) { + groupReference[group] = reader.readUint(params.bitsPerValue); + } + return groupReference; +} + +function readGroupWidths( + reader: BitReader, + params: ReturnType +) { + const groupWidth = new Uint32Array(params.groupCount); + for (let group = 0; group < params.groupCount; group++) { + groupWidth[group] = + params.referenceGroupWidth + reader.readUint(params.bitsGroupWidth); + } + return groupWidth; +} + +function readGroupLengths( + reader: BitReader, + params: ReturnType +) { + const groupLength = new Uint32Array(params.groupCount); + for (let group = 0; group < params.groupCount - 1; group++) { + groupLength[group] = + params.referenceGroupLength + + reader.readUint(params.bitsScaledLength) * params.lengthIncrement; + } + groupLength[params.groupCount - 1] = params.lastGroupLength; + return groupLength; +} + +function unpackComplexIntegerValues( + reader: BitReader, + params: ReturnType, + groupReference: Uint32Array, + groupWidth: Uint32Array, + groupLength: Uint32Array +) { + const values = new Int32Array(params.valueCount); + let outIndex = 0; + for (let group = 0; group < params.groupCount; group++) { + const width = groupWidth[group]; + const length = groupLength[group]; + for (let i = 0; i < length; i++) { + const value = width > 0 ? reader.readUint(width) : 0; + values[outIndex++] = groupReference[group] + value; + } + } + return values; +} + +function restoreSpatialDifferences( + values: Int32Array, + initialValues: number[], + overallMinimum: number, + orderSpatialDifferencing: number +) { + if (orderSpatialDifferencing === 1) { + values[0] = initialValues[0]; + for (let i = 1; i < values.length; i++) { + values[i] += overallMinimum + values[i - 1]; + } + return; + } + if (orderSpatialDifferencing === 2) { + values[0] = initialValues[0]; + values[1] = initialValues[1]; + for (let i = 2; i < values.length; i++) { + values[i] += overallMinimum + 2 * values[i - 1] - values[i - 2]; + } + return; + } + throw new Error( + `GRIB2 spatial differencing order ${orderSpatialDifferencing} is not supported` + ); +} + +function unpackSimple( + section6: TSection, + section7: TSection, + view: DataView, + referenceValue: number, + scale: number, + bitsPerValue: number, + valueCount: number, + totalPoints: number +) { + const out = new Float64Array(valueCount); + if (bitsPerValue === 0) { + out.fill(referenceValue); + return applyBitmap(section6, view, out, totalPoints); + } + + const reader = createSection7Reader(section7, view); + for (let i = 0; i < valueCount; i++) { + out[i] = referenceValue + reader.readUint(bitsPerValue) * scale; + } + return applyBitmap(section6, view, out, totalPoints); +} + +function applyBitmap( + section6: TSection, + view: DataView, + packedValues: Float64Array, + totalPoints: number +) { + const bitmapIndicator = view.getUint8(section6.offset + 5); + if (bitmapIndicator === 255) { + return packedValues; + } + if (bitmapIndicator !== 0) { + throw new Error( + `GRIB2 bitmap indicator ${bitmapIndicator} is not supported` + ); + } + + const bitmapStart = section6.offset + 6; + const bitmapByteCount = section6.length - 6; + const out = new Float64Array(totalPoints).fill(NaN); + const reader = new BitReader( + new Uint8Array(view.buffer, view.byteOffset + bitmapStart, bitmapByteCount) + ); + let packedIndex = 0; + for (let i = 0; i < totalPoints; i++) { + if (reader.readUint(1) === 1) { + out[i] = packedValues[packedIndex++]; + } + } + return out; +} + +function parseGrib2Sections(buf: Uint8Array) { + const view = createDataView(buf); + const sections: TSection[] = [{ length: 16, number: 0, offset: 0 }]; + let position = 16; + + while (position < buf.byteLength - 4) { + if ( + buf[position] === 0x37 && + buf[position + 1] === 0x37 && + buf[position + 2] === 0x37 && + buf[position + 3] === 0x37 + ) { + sections.push({ length: 4, number: 8, offset: position }); + break; + } + + const sectionLength = view.getUint32(position, false); + sections.push({ + length: sectionLength, + number: view.getUint8(position + 4), + offset: position, + }); + position += sectionLength; + } + + return { sections, view }; +} + +function createSection7Reader(section7: TSection, view: DataView) { + return new BitReader(createSection7Bytes(section7, view)); +} + +function createSection7Bytes(section7: TSection, view: DataView) { + const dataStart = section7.offset + 5; + return new Uint8Array( + view.buffer, + view.byteOffset + dataStart, + section7.length - 5 + ); +} + +class AecDecoder { + readonly #idLength: number; + readonly #output: Float64Array; + readonly #params: TAecParams; + readonly #reader: BitReader; + #blockIndexWithinRsi = 0; + #predictor: number | undefined; + #sampleIndex = 0; + + constructor(input: Uint8Array, params: TAecParams, outputSamples: number) { + validateAecParams(params); + this.#idLength = getAecIdLength(params); + this.#output = new Float64Array(outputSamples); + this.#params = params; + this.#reader = new BitReader(input); + } + + decode() { + while (this.#sampleIndex < this.#output.length) { + this.#decodeBlock(); + } + return this.#output; + } + + #decodeBlock() { + if (this.#isRsiStart()) { + this.#predictor = undefined; + } + + const id = this.#reader.readUint(this.#idLength); + const maxId = Math.pow(2, this.#idLength) - 1; + if (id === 0) { + this.#decodeLowEntropyBlock(); + } else if (id === maxId) { + this.#decodeUncompressedBlock(); + } else { + this.#decodeSplitBlock(id - 1); + } + } + + #decodeLowEntropyBlock() { + const selector = this.#reader.readBit(); + const referenceSampleConsumed = this.#consumeRsiReferenceIfNeeded(); + const remaining = this.#remainingBlockSamples(referenceSampleConsumed); + + if (!selector) { + this.#decodeZeroRun(referenceSampleConsumed); + return; + } + + this.#decodeSecondExtension(remaining, referenceSampleConsumed); + this.#advanceBlockIndex(1); + } + + #decodeZeroRun(referenceSampleConsumed: boolean) { + const blockCount = this.#readZeroBlockCount(); + const sampleCount = + blockCount * this.#params.blockSize - Number(referenceSampleConsumed); + this.#emitRepeatedValue(0, sampleCount); + this.#advanceBlockIndex(blockCount); + } + + #decodeUncompressedBlock() { + const referenceSampleConsumed = this.#consumeRsiReferenceIfNeeded(); + const remaining = this.#remainingBlockSamples(referenceSampleConsumed); + for (let i = 0; i < remaining && this.#hasOutputSpace(); i++) { + this.#emitValue(this.#reader.readUint(this.#params.bitsPerSample)); + } + this.#advanceBlockIndex(1); + } + + #decodeSplitBlock(k: number) { + const referenceSampleConsumed = this.#consumeRsiReferenceIfNeeded(); + const sampleCount = Math.min( + this.#remainingBlockSamples(referenceSampleConsumed), + this.#output.length - this.#sampleIndex + ); + const values = this.#readSplitValues(k, sampleCount); + + for (const value of values) { + this.#emitValue(value); + } + this.#advanceBlockIndex(1); + } + + #readSplitValues(k: number, sampleCount: number) { + const values = new Uint32Array(sampleCount); + const multiplier = Math.pow(2, k); + for (let i = 0; i < sampleCount; i++) { + values[i] = readUnary(this.#reader) * multiplier; + } + for (let i = 0; k > 0 && i < sampleCount; i++) { + values[i] += this.#reader.readUint(k); + } + return values; + } + + #decodeSecondExtension(remaining: number, referenceSampleConsumed: boolean) { + let remainingInBlock = remaining; + let emitOddFirst = referenceSampleConsumed; + + while (remainingInBlock > 0 && this.#hasOutputSpace()) { + const [evenValue, oddValue] = readSecondExtensionPair(this.#reader); + if (emitOddFirst) { + this.#emitValue(oddValue); + remainingInBlock -= 1; + emitOddFirst = false; + continue; + } + this.#emitValue(evenValue); + remainingInBlock -= 1; + if (remainingInBlock > 0 && this.#hasOutputSpace()) { + this.#emitValue(oddValue); + remainingInBlock -= 1; + } + } + } + + #readZeroBlockCount() { + const runOptionStart = 5; + const runLength = readUnary(this.#reader); + let blockCount = runLength + 1; + + if (blockCount === runOptionStart) { + blockCount = Math.min( + this.#params.rsi - this.#blockIndexWithinRsi, + 64 - (this.#blockIndexWithinRsi % 64) + ); + } else if (blockCount > runOptionStart) { + blockCount -= 1; + } + + return blockCount; + } + + #consumeRsiReferenceIfNeeded() { + if (!this.#isRsiStart()) { + return false; + } + + const raw = this.#reader.readUint(this.#params.bitsPerSample); + const value = hasAecFlag(this.#params, AecFlag.DATA_SIGNED) + ? signExtend(raw, this.#params.bitsPerSample) + : raw; + this.#output[this.#sampleIndex++] = value; + this.#predictor = value; + return true; + } + + #emitRepeatedValue(value: number, count: number) { + for (let i = 0; i < count && this.#hasOutputSpace(); i++) { + this.#emitValue(value); + } + } + + #emitValue(value: number) { + if (!this.#hasOutputSpace()) { + return; + } + + const sampleValue = hasAecFlag(this.#params, AecFlag.DATA_PREPROCESS) + ? this.#restorePreprocessedValue(value) + : this.#readRawSampleValue(value); + this.#output[this.#sampleIndex++] = sampleValue; + } + + #restorePreprocessedValue(value: number) { + if (this.#predictor === undefined) { + throw new Error("CCSDS/AEC stream is missing a reference sample"); + } + + const next = inverseAecPreprocess(this.#predictor, value, this.#params); + this.#predictor = next; + return next; + } + + #readRawSampleValue(value: number) { + return hasAecFlag(this.#params, AecFlag.DATA_SIGNED) + ? signExtend(value, this.#params.bitsPerSample) + : value; + } + + #remainingBlockSamples(referenceSampleConsumed: boolean) { + return this.#params.blockSize - Number(referenceSampleConsumed); + } + + #advanceBlockIndex(blockCount: number) { + this.#blockIndexWithinRsi += blockCount; + if (!hasAecFlag(this.#params, AecFlag.DATA_PREPROCESS)) { + return; + } + + if (this.#blockIndexWithinRsi >= this.#params.rsi) { + this.#blockIndexWithinRsi %= this.#params.rsi; + if (hasAecFlag(this.#params, AecFlag.PAD_RSI)) { + this.#reader.alignToByte(); + } + } + } + + #isRsiStart() { + return ( + hasAecFlag(this.#params, AecFlag.DATA_PREPROCESS) && + this.#blockIndexWithinRsi === 0 + ); + } + + #hasOutputSpace() { + return this.#sampleIndex < this.#output.length; + } +} + +function decodeAecSamples( + input: Uint8Array, + params: TAecParams, + outputSamples: number +) { + return new AecDecoder(input, params, outputSamples).decode(); +} + +function readUnary(reader: BitReader) { + let count = 0; + while (!reader.readBit()) { + count += 1; + if (count > 1_000_000) { + throw new Error("CCSDS/AEC unary run is too long"); + } + } + return count; +} + +function readSecondExtensionPair(reader: BitReader): [number, number] { + const value = readUnary(reader); + if (value > 90) { + throw new Error("CCSDS/AEC second-extension symbol is too large"); + } + + let index = 0; + for (let sum = 0; sum <= 12; sum++) { + for (let k = 0; k <= sum; k++) { + if (index === value) { + return [sum - k, k]; + } + index += 1; + } + } + + return [0, 0]; +} + +function inverseAecPreprocess( + previousValue: number, + encodedDelta: number, + params: TAecParams +) { + const delta = + encodedDelta % 2 === 0 ? encodedDelta / 2 : -((encodedDelta + 1) / 2); + const halfDelta = Math.floor(encodedDelta / 2) + (encodedDelta % 2); + + return hasAecFlag(params, AecFlag.DATA_SIGNED) + ? inverseSignedAecPreprocess( + previousValue, + encodedDelta, + delta, + halfDelta, + params + ) + : inverseUnsignedAecPreprocess( + previousValue, + encodedDelta, + delta, + halfDelta, + params + ); +} + +function inverseSignedAecPreprocess( + previousValue: number, + encodedDelta: number, + delta: number, + halfDelta: number, + params: TAecParams +) { + const signedMax = Math.pow(2, params.bitsPerSample - 1) - 1; + if (previousValue < 0) { + return halfDelta <= signedMax + previousValue + 1 + ? previousValue + delta + : encodedDelta - signedMax - 1; + } + return halfDelta <= signedMax - previousValue + ? previousValue + delta + : signedMax - encodedDelta; +} + +function inverseUnsignedAecPreprocess( + previousValue: number, + encodedDelta: number, + delta: number, + halfDelta: number, + params: TAecParams +) { + const unsignedMax = Math.pow(2, params.bitsPerSample) - 1; + const midpoint = Math.pow(2, params.bitsPerSample - 1); + const highMask = Math.floor(previousValue / midpoint) % 2 === 1; + const threshold = highMask ? unsignedMax - previousValue : previousValue; + + if (halfDelta <= threshold) { + return previousValue + delta; + } + return highMask ? unsignedMax - encodedDelta : encodedDelta; +} + +function scaleAecSamples( + sampleValues: Float64Array, + referenceValue: number, + scale: number +) { + const out = new Float64Array(sampleValues.length); + for (let i = 0; i < sampleValues.length; i++) { + out[i] = referenceValue + sampleValues[i] * scale; + } + return out; +} + +function aecFlagsFromGrib2(ccsdsFlags: number) { + let flags = 0; + + if ((ccsdsFlags & (1 << 0)) !== 0) { + flags |= AecFlag.DATA_SIGNED; + } + if ((ccsdsFlags & (1 << 1)) !== 0) { + flags |= AecFlag.DATA_3BYTE; + } + if ((ccsdsFlags & (1 << 2)) !== 0) { + flags |= AecFlag.MSB; + } + if ((ccsdsFlags & (1 << 3)) !== 0) { + flags |= AecFlag.DATA_PREPROCESS; + } + if ((ccsdsFlags & (1 << 4)) !== 0) { + flags |= AecFlag.RESTRICTED; + } + if ((ccsdsFlags & (1 << 5)) !== 0) { + flags |= AecFlag.PAD_RSI; + } + + return flags; +} + +function validateAecParams(params: TAecParams) { + if (params.bitsPerSample < 1 || params.bitsPerSample > 32) { + throw new Error("CCSDS/AEC bits per sample must be between 1 and 32"); + } + if (![8, 16, 32, 64].includes(params.blockSize)) { + throw new Error("CCSDS/AEC block size must be one of 8, 16, 32, or 64"); + } + if (params.rsi <= 0) { + throw new Error("CCSDS/AEC reference sample interval must be positive"); + } +} + +function getAecIdLength(params: TAecParams) { + if (hasAecFlag(params, AecFlag.RESTRICTED) && params.bitsPerSample <= 4) { + return params.bitsPerSample <= 2 ? 1 : 2; + } + if (params.bitsPerSample > 16) { + return 5; + } + return params.bitsPerSample > 8 ? 4 : 3; +} + +function hasAecFlag(params: TAecParams, flag: number) { + return (params.flags & flag) !== 0; +} + +function valuesToBytes( + values: Float64Array, + dataType: DataType, + endian: TEndian +) { + const typedValues = castValues(values, dataType); + const bytes = new Uint8Array( + typedValues.buffer, + typedValues.byteOffset, + typedValues.byteLength + ); + + if (endian === Endian.LITTLE || bytesPerElement(typedValues) === 1) { + return bytes; + } + + const swapped = new Uint8Array(bytes); + byteswapInPlace(swapped, bytesPerElement(typedValues)); + return swapped; +} + +function castValues(values: Float64Array, dataType: DataType): TDecodedArray { + switch (dataType) { + case "int8": + return Int8Array.from(values); + case "int16": + return Int16Array.from(values); + case "int32": + return Int32Array.from(values); + case "int64": + return BigInt64Array.from(values, (value) => BigInt(Math.trunc(value))); + case "uint8": + return Uint8Array.from(values); + case "uint16": + return Uint16Array.from(values); + case "uint32": + return Uint32Array.from(values); + case "uint64": + return BigUint64Array.from(values, (value) => BigInt(Math.trunc(value))); + case "float32": + return Float32Array.from(values); + case "float64": + return values; + default: + throw new Error(`GRIB values cannot be decoded as ${dataType}`); + } +} + +function isComplexPackingMissingValue( + value: number, + width: number, + missingManagement: number +) { + if (missingManagement === 0 || width === 0) { + return false; + } + + const allOnes = Math.pow(2, width) - 1; + if (value === allOnes) { + return true; + } + + return missingManagement === 2 && width > 1 && value === allOnes - 1; +} + +function getBytesEndian(meta?: TCodecMeta): TEndian { + const bytesCodec = meta?.codecs?.find((codec) => codec.name === "bytes"); + return bytesCodec?.configuration?.endian === Endian.BIG + ? Endian.BIG + : Endian.LITTLE; +} + +function createDataView(buf: Uint8Array) { + return new DataView(buf.buffer, buf.byteOffset, buf.byteLength); +} + +function readUint24(view: DataView, offset: number) { + return ( + (view.getUint8(offset) << 16) | + (view.getUint8(offset + 1) << 8) | + view.getUint8(offset + 2) + ); +} + +function readIbmFloat32(view: DataView, offset: number) { + const b0 = view.getUint8(offset); + const b1 = view.getUint8(offset + 1); + const b2 = view.getUint8(offset + 2); + const b3 = view.getUint8(offset + 3); + const sign = b0 & 0x80 ? -1 : 1; + const exponent = (b0 & 0x7f) - 64; + const mantissa = ((b1 << 16) | (b2 << 8) | b3) / 0x1000000; + return sign * mantissa * Math.pow(16, exponent); +} + +function signMagnitude(raw: number, bits: number) { + if (bits === 0) { + return 0; + } + + const signBit = Math.pow(2, bits - 1); + if (raw >= signBit) { + return -(raw - signBit); + } + return raw; +} + +function signExtend(raw: number, bits: number) { + const signBit = Math.pow(2, bits - 1); + return raw >= signBit ? raw - Math.pow(2, bits) : raw; +} + +function bytesPerElement(values: TDecodedArray) { + return values.BYTES_PER_ELEMENT; +} + +function byteswapInPlace(bytes: Uint8Array, bytesPerElementValue: number) { + for ( + let offset = 0; + offset < bytes.byteLength; + offset += bytesPerElementValue + ) { + for (let i = 0; i < bytesPerElementValue / 2; i++) { + const left = offset + i; + const right = offset + bytesPerElementValue - 1 - i; + const tmp = bytes[left]; + bytes[left] = bytes[right]; + bytes[right] = tmp; + } + } +} diff --git a/src/lib/data/irregularGridHelpers.ts b/src/lib/data/irregularGridHelpers.ts index 7248560..b9dd625 100644 --- a/src/lib/data/irregularGridHelpers.ts +++ b/src/lib/data/irregularGridHelpers.ts @@ -46,9 +46,22 @@ export function reconcileCoordinates( ): { latitudes: Float32Array; longitudes: Float32Array } { const latitudes = latitudesVar.data as Float32Array; const longitudes = longitudesVar.data as Float32Array; - const latShape = latitudesVar.shape; - const lonShape = longitudesVar.shape; + return reconcileCoordinateArrays( + latitudes, + longitudes, + latitudesVar.shape, + longitudesVar.shape, + dataLength + ); +} +export function reconcileCoordinateArrays( + latitudes: Float32Array, + longitudes: Float32Array, + latitudeShape: number[], + longitudeShape: number[], + dataLength: number +): { latitudes: Float32Array; longitudes: Float32Array } { // Case 1: All arrays already have same length (typical irregular grid) if (latitudes.length === dataLength && longitudes.length === dataLength) { return { latitudes, longitudes }; @@ -56,9 +69,9 @@ export function reconcileCoordinates( // Case 2: 2D coordinate arrays - just ensure they're flattened and match data // (zarr.Chunk.data should already be flattened) - if (latShape.length === 2 && lonShape.length === 2) { - const latTotal = latShape[0] * latShape[1]; - const lonTotal = lonShape[0] * lonShape[1]; + if (latitudeShape.length === 2 && longitudeShape.length === 2) { + const latTotal = latitudeShape[0] * latitudeShape[1]; + const lonTotal = longitudeShape[0] * longitudeShape[1]; if (latTotal === dataLength && lonTotal === dataLength) { // Already flattened by zarr, should match return { latitudes, longitudes }; @@ -66,7 +79,7 @@ export function reconcileCoordinates( } // Case 3: 1D separate coordinates - expand via meshgrid - if (latShape.length === 1 && lonShape.length === 1) { + if (latitudeShape.length === 1 && longitudeShape.length === 1) { const product = latitudes.length * longitudes.length; if (product === dataLength) { return expandCoordinatesToMeshgrid(latitudes, longitudes, dataLength); diff --git a/src/lib/grids/curvilinear.worker.ts b/src/lib/grids/curvilinear.worker.ts new file mode 100644 index 0000000..308a825 --- /dev/null +++ b/src/lib/grids/curvilinear.worker.ts @@ -0,0 +1,113 @@ +/// + +import { + buildCurvilinearBatch, + buildCurvilinearHoverIndexData, + detectCurvilinearColumnPeriodicity, + detectCurvilinearLongitudeFlip, + getCurvilinearBatchCount, + type TCurvilinearGrid, +} from "./curvilinearCalculations.ts"; +import type { + TCurvilinearWorkerRequest, + TCurvilinearWorkerResponse, +} from "./curvilinearWorkerProtocol.ts"; + +import { GridGeometryWorkerMessageType } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import { + postGridGeometryBatch, + postGridGeometryHoverIndex, + postGridGeometryResponse, +} from "@/lib/grids/gridGeometryWorkerUtils.ts"; +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; + +function postResponse( + response: TCurvilinearWorkerResponse, + transfer: Transferable[] = [] +) { + postGridGeometryResponse(workerScope, response, transfer); +} + +function createGrid(request: TCurvilinearWorkerRequest): TCurvilinearGrid { + return { + latitudes: request.latitudes, + longitudes: request.longitudes, + data: request.data, + nj: request.nj, + ni: request.ni, + shouldFlipLongitude: detectCurvilinearLongitudeFlip( + request.longitudes, + request.latitudes, + request.missingValue, + request.fillValue, + request.nj, + request.ni + ), + isPeriodicI: detectCurvilinearColumnPeriodicity( + request.longitudes, + request.nj, + request.ni + ), + }; +} + +function postHoverIndex(requestId: number, grid: TCurvilinearGrid) { + postGridGeometryHoverIndex( + workerScope, + requestId, + buildCurvilinearHoverIndexData(grid) + ); +} + +function postBatches( + request: TCurvilinearWorkerRequest, + grid: TCurvilinearGrid, + totalBatches: number +) { + const projection = new ProjectionHelper( + request.projectionType, + request.projectionCenter + ); + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + const batch = buildCurvilinearBatch( + grid, + batchIndex, + request.batchSize, + projection + ); + postGridGeometryBatch(workerScope, request.requestId, batch); + } +} + +function buildGrid(request: TCurvilinearWorkerRequest) { + const grid = createGrid(request); + const totalBatches = getCurvilinearBatchCount(grid, request.batchSize); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { + totalBatches, + shouldFlipLongitude: grid.shouldFlipLongitude, + }, + }); + postHoverIndex(request.requestId, grid); + postBatches(request, grid, totalBatches); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.DONE, + }); +} + +workerScope.onmessage = (event: MessageEvent) => { + try { + buildGrid(event.data); + } catch (error) { + postResponse({ + requestId: event.data.requestId, + type: GridGeometryWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }); + } +}; diff --git a/src/lib/grids/curvilinearCalculations.ts b/src/lib/grids/curvilinearCalculations.ts new file mode 100644 index 0000000..69d74b9 --- /dev/null +++ b/src/lib/grids/curvilinearCalculations.ts @@ -0,0 +1,420 @@ +import { buildSerializedGeoSampleIndexData } from "./gridWorkerCalculations.ts"; +import type { + TGridGeometryBatch, + TSerializedGeoSampleIndexData, +} from "./gridWorkerTypes.ts"; + +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +type TCoordinate = { lat: number; lon: number }; +type TBoundaryFunction = ( + row: number, + fromColumn: number, + toColumn: number +) => TCoordinate; +type TMidpointFunction = (fromColumn: number, toColumn: number) => TCoordinate; + +export type TCurvilinearGrid = { + latitudes: Float32Array; + longitudes: Float32Array; + data: Float32Array; + nj: number; + ni: number; + shouldFlipLongitude: boolean; + isPeriodicI: boolean; +}; + +function isMissingCoordinate( + value: number, + missingValue: number, + fillValue: number +) { + return Number.isNaN(value) || value === missingValue || value === fillValue; +} + +export function detectCurvilinearLongitudeFlip( + longitudes: Float32Array, + latitudes: Float32Array, + missingValue: number, + fillValue: number, + nj: number, + ni: number +) { + for (let row = 0; row < nj - 1; row++) { + const columnLimit = Math.min(ni - 1, 10); + for (let column = 0; column < columnLimit; column++) { + const topLeft = row * ni + column; + const topRight = topLeft + 1; + const bottomLeft = (row + 1) * ni + column; + const bottomRight = bottomLeft + 1; + if ( + isMissingCoordinate(longitudes[topLeft], missingValue, fillValue) || + isMissingCoordinate(longitudes[topRight], missingValue, fillValue) || + isMissingCoordinate(longitudes[bottomLeft], missingValue, fillValue) || + isMissingCoordinate(longitudes[bottomRight], missingValue, fillValue) + ) { + continue; + } + const longitudeI = longitudes[topRight] - longitudes[topLeft]; + const latitudeI = latitudes[topRight] - latitudes[topLeft]; + const longitudeJ = longitudes[bottomLeft] - longitudes[topLeft]; + const latitudeJ = latitudes[bottomLeft] - latitudes[topLeft]; + return longitudeI * latitudeJ - latitudeI * longitudeJ < 0; + } + } + return false; +} + +export function detectCurvilinearColumnPeriodicity( + longitudes: Float32Array, + nj: number, + ni: number +) { + if (ni < 3) { + return false; + } + const sampleRow = Math.floor(nj / 2); + const sampleCount = Math.min(ni - 1, 10); + let spacingTotal = 0; + let validSpacingCount = 0; + for (let column = 0; column < sampleCount; column++) { + const longitude = longitudes[sampleRow * ni + column]; + const nextLongitude = longitudes[sampleRow * ni + column + 1]; + const spacing = Math.abs(((nextLongitude - longitude + 540) % 360) - 180); + if (spacing > 0) { + spacingTotal += spacing; + validSpacingCount++; + } + } + if (validSpacingCount === 0) { + return false; + } + const firstLongitude = longitudes[sampleRow * ni]; + const lastLongitude = longitudes[sampleRow * ni + ni - 1]; + const wrapGap = Math.abs( + ((firstLongitude - lastLongitude + 540) % 360) - 180 + ); + return wrapGap < (spacingTotal / validSpacingCount) * 4; +} + +function getNextColumn(column: number, ni: number, flipLongitude: boolean) { + if (flipLongitude) { + return column === 0 ? ni - 1 : column - 1; + } + return (column + 1) % ni; +} + +function getPreviousColumn(column: number, ni: number, flipLongitude: boolean) { + if (flipLongitude) { + return (column + 1) % ni; + } + return column === 0 ? ni - 1 : column - 1; +} + +function getCircularMeanLongitude( + first: number, + second: number, + third: number, + fourth: number +) { + const radians = Math.PI / 180; + const sine = + Math.sin(first * radians) + + Math.sin(second * radians) + + Math.sin(third * radians) + + Math.sin(fourth * radians); + const cosine = + Math.cos(first * radians) + + Math.cos(second * radians) + + Math.cos(third * radians) + + Math.cos(fourth * radians); + return Math.atan2(sine / 4, cosine / 4) / radians; +} + +function getCellCenter( + grid: TCurvilinearGrid, + topLeft: number, + topRight: number, + bottomLeft: number, + bottomRight: number +) { + return { + lat: + (grid.latitudes[topLeft] + + grid.latitudes[topRight] + + grid.latitudes[bottomLeft] + + grid.latitudes[bottomRight]) / + 4, + lon: getCircularMeanLongitude( + grid.longitudes[topLeft], + grid.longitudes[topRight], + grid.longitudes[bottomLeft], + grid.longitudes[bottomRight] + ), + }; +} + +function getBoundaryPoint( + grid: TCurvilinearGrid, + row: number, + fromColumn: number, + toColumn: number +) { + return getCellCenter( + grid, + row * grid.ni + fromColumn, + row * grid.ni + toColumn, + (row + 1) * grid.ni + fromColumn, + (row + 1) * grid.ni + toColumn + ); +} + +function getRowMidpoint( + grid: TCurvilinearGrid, + row: number, + fromColumn: number, + toColumn: number +) { + return getCellCenter( + grid, + row * grid.ni + fromColumn, + row * grid.ni + toColumn, + row * grid.ni + fromColumn, + row * grid.ni + toColumn + ); +} + +function mirrorBoundaryPoint(center: TCoordinate, boundary: TCoordinate) { + const longitudeOffset = ((boundary.lon - center.lon + 540) % 360) - 180; + return { + lat: 2 * center.lat - boundary.lat, + lon: center.lon - longitudeOffset, + }; +} + +function getPreviousBoundaryCorners( + row: number, + column: number, + nextColumn: number, + boundary: TBoundaryFunction, + midpoint: TMidpointFunction +) { + const forwardNext = boundary(row, column, nextColumn); + const forwardPrevious = mirrorBoundaryPoint( + midpoint(column, column), + forwardNext + ); + const backwardNext = + row === 0 + ? mirrorBoundaryPoint(midpoint(column, nextColumn), forwardNext) + : boundary(row - 1, column, nextColumn); + return { + forwardPrevious, + forwardNext, + backwardPrevious: mirrorBoundaryPoint( + midpoint(column, column), + backwardNext + ), + backwardNext, + }; +} + +function getNextBoundaryCorners( + row: number, + column: number, + previousColumn: number, + boundary: TBoundaryFunction, + midpoint: TMidpointFunction +) { + const forwardPrevious = boundary(row, previousColumn, column); + const forwardNext = mirrorBoundaryPoint( + midpoint(column, column), + forwardPrevious + ); + const backwardPrevious = + row === 0 + ? mirrorBoundaryPoint(midpoint(previousColumn, column), forwardPrevious) + : boundary(row - 1, previousColumn, column); + return { + forwardPrevious, + forwardNext, + backwardPrevious, + backwardNext: mirrorBoundaryPoint( + midpoint(column, column), + backwardPrevious + ), + }; +} + +function getPeriodicCorners( + row: number, + column: number, + previousColumn: number, + nextColumn: number, + boundary: TBoundaryFunction, + midpoint: TMidpointFunction +) { + const forwardPrevious = boundary(row, previousColumn, column); + const forwardNext = boundary(row, column, nextColumn); + return { + forwardPrevious, + forwardNext, + backwardPrevious: + row === 0 + ? mirrorBoundaryPoint(midpoint(previousColumn, column), forwardPrevious) + : boundary(row - 1, previousColumn, column), + backwardNext: + row === 0 + ? mirrorBoundaryPoint(midpoint(column, nextColumn), forwardNext) + : boundary(row - 1, column, nextColumn), + }; +} + +/* eslint-disable-next-line max-lines-per-function */ +function getCenteredCellCorners( + grid: TCurvilinearGrid, + row: number, + column: number +) { + const previousColumn = getPreviousColumn( + column, + grid.ni, + grid.shouldFlipLongitude + ); + const nextColumn = getNextColumn(column, grid.ni, grid.shouldFlipLongitude); + const boundary: TBoundaryFunction = (boundaryRow, fromColumn, toColumn) => + getBoundaryPoint(grid, boundaryRow, fromColumn, toColumn); + const midpoint: TMidpointFunction = (fromColumn, toColumn) => + getRowMidpoint(grid, row, fromColumn, toColumn); + const isFirst = grid.shouldFlipLongitude + ? column === grid.ni - 1 + : column === 0; + const isLast = grid.shouldFlipLongitude + ? column === 0 + : column === grid.ni - 1; + const corners = + !grid.isPeriodicI && isFirst + ? getPreviousBoundaryCorners(row, column, nextColumn, boundary, midpoint) + : !grid.isPeriodicI && isLast + ? getNextBoundaryCorners( + row, + column, + previousColumn, + boundary, + midpoint + ) + : getPeriodicCorners( + row, + column, + previousColumn, + nextColumn, + boundary, + midpoint + ); + return { + latitudes: [ + corners.backwardPrevious.lat, + corners.backwardNext.lat, + corners.forwardNext.lat, + corners.forwardPrevious.lat, + ], + longitudes: [ + corners.backwardPrevious.lon, + corners.backwardNext.lon, + corners.forwardNext.lon, + corners.forwardPrevious.lon, + ], + }; +} + +function fillCell( + grid: TCurvilinearGrid, + projection: ProjectionHelper, + row: number, + column: number, + cellIndex: number, + batch: TGridGeometryBatch +) { + const corners = getCenteredCellCorners(grid, row, column); + const invalid = + corners.latitudes.some((value) => !Number.isFinite(value)) || + corners.longitudes.some((value) => !Number.isFinite(value)); + if (invalid) { + batch.dataValues.fill(NaN, cellIndex * 4, cellIndex * 4 + 4); + } else { + for (let corner = 0; corner < 4; corner++) { + projection.projectLatLonToArrays( + corners.latitudes[corner], + corners.longitudes[corner], + batch.positionValues, + cellIndex * 12 + corner * 3, + batch.latLonValues, + cellIndex * 8 + corner * 2 + ); + } + const dataIndex = row * grid.ni + column; + batch.dataValues.fill( + grid.data[dataIndex], + cellIndex * 4, + cellIndex * 4 + 4 + ); + } + const vertex = cellIndex * 4; + batch.indices.set( + [vertex, vertex + 1, vertex + 2, vertex, vertex + 2, vertex + 3], + cellIndex * 6 + ); +} + +export function getCurvilinearBatchCount( + grid: TCurvilinearGrid, + batchSize: number +) { + return Math.ceil((grid.nj - 1) / batchSize); +} + +export function buildCurvilinearBatch( + grid: TCurvilinearGrid, + batchIndex: number, + batchSize: number, + projection: ProjectionHelper +): TGridGeometryBatch { + const rowStart = batchIndex * batchSize; + const rowEnd = Math.min(rowStart + batchSize, grid.nj - 1); + const cellCount = (rowEnd - rowStart) * grid.ni; + const batch: TGridGeometryBatch = { + batchIndex, + positionValues: new Float32Array(cellCount * 12), + dataValues: new Float32Array(cellCount * 4), + latLonValues: new Float32Array(cellCount * 8), + indices: new Uint32Array(cellCount * 6), + }; + let cellIndex = 0; + for (let row = rowStart; row < rowEnd; row++) { + for (let column = 0; column < grid.ni; column++) { + fillCell(grid, projection, row, column, cellIndex, batch); + cellIndex++; + } + } + return batch; +} + +export function buildCurvilinearHoverIndexData( + grid: TCurvilinearGrid +): TSerializedGeoSampleIndexData { + const sampleCount = (grid.nj - 1) * grid.ni; + const latitudes = new Float64Array(sampleCount); + const longitudes = new Float64Array(sampleCount); + const values = new Float32Array(sampleCount); + let sampleIndex = 0; + for (let row = 0; row < grid.nj - 1; row++) { + for (let column = 0; column < grid.ni; column++) { + const dataIndex = row * grid.ni + column; + latitudes[sampleIndex] = grid.latitudes[dataIndex]; + longitudes[sampleIndex] = grid.longitudes[dataIndex]; + values[sampleIndex] = grid.data[dataIndex]; + sampleIndex++; + } + } + return buildSerializedGeoSampleIndexData(latitudes, longitudes, values); +} diff --git a/src/lib/grids/curvilinearWorkerClient.ts b/src/lib/grids/curvilinearWorkerClient.ts new file mode 100644 index 0000000..c8d0711 --- /dev/null +++ b/src/lib/grids/curvilinearWorkerClient.ts @@ -0,0 +1,69 @@ +import type { + TCurvilinearWorkerMetadata, + TCurvilinearWorkerRequest, +} from "./curvilinearWorkerProtocol.ts"; +import { + copyGridWorkerArray, + createGridGeometryWorkerClient, +} from "./gridGeometryWorkerClient.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import type { TGridGeometryBatch } from "./gridWorkerTypes.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +type TCurvilinearBuildRequest = { + latitudes: Float32Array; + longitudes: Float32Array; + data: Float32Array; + nj: number; + ni: number; + batchSize: number; + missingValue: number; + fillValue: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +const client = createGridGeometryWorkerClient< + TCurvilinearWorkerRequest, + TCurvilinearWorkerMetadata +>( + () => + new Worker(new URL("./curvilinear.worker.ts", import.meta.url), { + type: "module", + }) +); + +export function buildCurvilinearGrid( + request: TCurvilinearBuildRequest, + callbacks: { + onMetadata: (metadata: TCurvilinearWorkerMetadata) => void; + onBatch: (batch: TGridGeometryBatch) => void; + } +) { + return client.build((requestId) => { + const latitudes = copyGridWorkerArray(request.latitudes); + const longitudes = copyGridWorkerArray(request.longitudes); + const data = copyGridWorkerArray(request.data); + return { + message: { + ...request, + requestId, + type: GridGeometryWorkerMessageType.BUILD, + latitudes, + longitudes, + data, + projectionCenter: { + lat: request.projectionCenter.lat, + lon: request.projectionCenter.lon, + }, + }, + transfer: [latitudes.buffer, longitudes.buffer, data.buffer], + }; + }, callbacks); +} + +export const terminateCurvilinearWorker = client.terminate; diff --git a/src/lib/grids/curvilinearWorkerProtocol.ts b/src/lib/grids/curvilinearWorkerProtocol.ts new file mode 100644 index 0000000..b630f58 --- /dev/null +++ b/src/lib/grids/curvilinearWorkerProtocol.ts @@ -0,0 +1,29 @@ +import type { TGridGeometryWorkerResponse } from "./gridGeometryWorkerProtocol.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +export type TCurvilinearWorkerRequest = { + requestId: number; + type: "build"; + latitudes: Float32Array; + longitudes: Float32Array; + data: Float32Array; + nj: number; + ni: number; + batchSize: number; + missingValue: number; + fillValue: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +export type TCurvilinearWorkerMetadata = { + totalBatches: number; + shouldFlipLongitude: boolean; +}; + +export type TCurvilinearWorkerResponse = + TGridGeometryWorkerResponse; diff --git a/src/lib/grids/gaussianReduced.worker.ts b/src/lib/grids/gaussianReduced.worker.ts new file mode 100644 index 0000000..d1c0b3d --- /dev/null +++ b/src/lib/grids/gaussianReduced.worker.ts @@ -0,0 +1,88 @@ +/// + +import { + buildGaussianReducedBatch, + buildGaussianReducedHoverIndexData, + buildGaussianReducedRows, + getGaussianReducedBatchCount, +} from "./gaussianReducedCalculations.ts"; +import { + type TGaussianReducedWorkerRequest, + type TGaussianReducedWorkerResponse, +} from "./gaussianReducedWorkerProtocol.ts"; + +import { GridGeometryWorkerMessageType } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import { + postGridGeometryBatch, + postGridGeometryHoverIndex, + postGridGeometryResponse, +} from "@/lib/grids/gridGeometryWorkerUtils.ts"; +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; + +function postResponse( + response: TGaussianReducedWorkerResponse, + transfer: Transferable[] = [] +) { + postGridGeometryResponse(workerScope, response, transfer); +} + +function postBatches( + request: TGaussianReducedWorkerRequest, + grid: ReturnType, + totalBatches: number +) { + const projection = new ProjectionHelper( + request.projectionType, + request.projectionCenter + ); + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + const batch = buildGaussianReducedBatch( + grid, + batchIndex, + request.batchSize, + request.epsilon, + projection + ); + postGridGeometryBatch(workerScope, request.requestId, batch); + } +} + +function buildGrid(request: TGaussianReducedWorkerRequest) { + const grid = buildGaussianReducedRows( + request.latitudes, + request.longitudes, + request.data + ); + const totalBatches = getGaussianReducedBatchCount(grid, request.batchSize); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { totalBatches }, + }); + postGridGeometryHoverIndex( + workerScope, + request.requestId, + buildGaussianReducedHoverIndexData(grid) + ); + postBatches(request, grid, totalBatches); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.DONE, + }); +} + +workerScope.onmessage = ( + event: MessageEvent +) => { + try { + buildGrid(event.data); + } catch (error) { + postResponse({ + requestId: event.data.requestId, + type: GridGeometryWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }); + } +}; diff --git a/src/lib/grids/gaussianReducedCalculations.ts b/src/lib/grids/gaussianReducedCalculations.ts new file mode 100644 index 0000000..3be2ca4 --- /dev/null +++ b/src/lib/grids/gaussianReducedCalculations.ts @@ -0,0 +1,195 @@ +import { buildSerializedGeoSampleIndexData } from "./gridWorkerCalculations.ts"; +import type { + TGridGeometryBatch, + TSerializedGeoSampleIndexData, +} from "./gridWorkerTypes.ts"; + +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +type TGaussianReducedCell = { + lon: number; + value: number; +}; + +export type TGaussianReducedRows = { + rows: Record; + uniqueLatitudes: number[]; +}; + +export function buildGaussianReducedRows( + latitudes: Float64Array, + longitudes: Float64Array, + data: Float32Array +): TGaussianReducedRows { + const rows: Record = {}; + for (let index = 0; index < latitudes.length; index++) { + const latitude = latitudes[index]; + if (!rows[latitude]) { + rows[latitude] = []; + } + rows[latitude].push({ lon: longitudes[index], value: data[index] }); + } + + const uniqueLatitudes = Object.keys(rows) + .map(Number) + .sort((first, second) => second - first); + return { rows, uniqueLatitudes }; +} + +function getCellData(row: TGaussianReducedCell[], index: number) { + const cell = row[index]; + const nextCell = row[(index + 1) % row.length]; + const longitudeWidth = (nextCell.lon - cell.lon + 360) % 360; + return { cell, longitudeWidth }; +} + +function countCells( + grid: TGaussianReducedRows, + latitudeStart: number, + latitudeEnd: number +) { + let totalCells = 0; + for (let index = latitudeStart; index < latitudeEnd; index++) { + totalCells += grid.rows[grid.uniqueLatitudes[index]].length; + } + return totalCells; +} + +export function getGaussianReducedBatchCount( + grid: TGaussianReducedRows, + batchSize: number +) { + return Math.ceil((grid.uniqueLatitudes.length - 1) / batchSize); +} + +function projectQuad( + projection: ProjectionHelper, + latitudeTop: number, + longitudeLeft: number, + latitudeBottom: number, + longitudeWidth: number, + epsilon: number, + positionValues: Float32Array, + latLonValues: Float32Array, + positionOffset: number, + latLonOffset: number +) { + const longitudeRight = longitudeLeft - longitudeWidth; + projection.projectLatLonToArrays( + latitudeTop, + longitudeLeft + epsilon, + positionValues, + positionOffset, + latLonValues, + latLonOffset + ); + projection.projectLatLonToArrays( + latitudeTop, + longitudeRight - epsilon, + positionValues, + positionOffset + 3, + latLonValues, + latLonOffset + 2 + ); + projection.projectLatLonToArrays( + latitudeBottom - epsilon, + longitudeRight - epsilon, + positionValues, + positionOffset + 6, + latLonValues, + latLonOffset + 4 + ); + projection.projectLatLonToArrays( + latitudeBottom - epsilon, + longitudeLeft + epsilon, + positionValues, + positionOffset + 9, + latLonValues, + latLonOffset + 6 + ); +} + +/* eslint-disable-next-line max-lines-per-function */ +export function buildGaussianReducedBatch( + grid: TGaussianReducedRows, + batchIndex: number, + batchSize: number, + epsilon: number, + projection: ProjectionHelper +): TGridGeometryBatch { + const latitudeStart = batchIndex * batchSize; + const latitudeEnd = Math.min( + latitudeStart + batchSize, + grid.uniqueLatitudes.length - 1 + ); + const totalCells = countCells(grid, latitudeStart, latitudeEnd); + const latLonValues = new Float32Array(totalCells * 8); + const positionValues = new Float32Array(totalCells * 12); + const dataValues = new Float32Array(totalCells * 4); + const indices = new Uint32Array(totalCells * 6); + let cellIndex = 0; + + for ( + let latitudeIndex = latitudeStart; + latitudeIndex < latitudeEnd; + latitudeIndex++ + ) { + const latitudeTop = grid.uniqueLatitudes[latitudeIndex]; + const latitudeBottom = grid.uniqueLatitudes[latitudeIndex + 1]; + const row = grid.rows[latitudeTop]; + + for (let index = 0; index < row.length; index++) { + const { cell, longitudeWidth } = getCellData(row, index); + projectQuad( + projection, + latitudeTop, + cell.lon, + latitudeBottom, + longitudeWidth, + epsilon, + positionValues, + latLonValues, + cellIndex * 12, + cellIndex * 8 + ); + dataValues.fill(cell.value, cellIndex * 4, cellIndex * 4 + 4); + const vertex = cellIndex * 4; + indices.set( + [vertex, vertex + 1, vertex + 2, vertex, vertex + 2, vertex + 3], + cellIndex * 6 + ); + cellIndex++; + } + } + + return { batchIndex, positionValues, dataValues, latLonValues, indices }; +} + +export function buildGaussianReducedHoverIndexData( + grid: TGaussianReducedRows +): TSerializedGeoSampleIndexData { + const rowCount = Math.max(0, grid.uniqueLatitudes.length - 1); + const sampleCount = countCells(grid, 0, rowCount); + const latitudes = new Float64Array(sampleCount); + const longitudes = new Float64Array(sampleCount); + const values = new Float32Array(sampleCount); + let sampleIndex = 0; + + for (let latitudeIndex = 0; latitudeIndex < rowCount; latitudeIndex++) { + const latitudeTop = grid.uniqueLatitudes[latitudeIndex]; + const latitudeBottom = grid.uniqueLatitudes[latitudeIndex + 1]; + const row = grid.rows[latitudeTop]; + for (let indexInRow = 0; indexInRow < row.length; indexInRow++) { + const { cell, longitudeWidth } = getCellData(row, indexInRow); + const latitude = (latitudeTop + latitudeBottom) / 2; + const longitude = ProjectionHelper.normalizeLongitude( + cell.lon - longitudeWidth / 2 + ); + latitudes[sampleIndex] = latitude; + longitudes[sampleIndex] = longitude; + values[sampleIndex] = cell.value; + sampleIndex++; + } + } + return buildSerializedGeoSampleIndexData(latitudes, longitudes, values); +} diff --git a/src/lib/grids/gaussianReducedWorkerClient.ts b/src/lib/grids/gaussianReducedWorkerClient.ts new file mode 100644 index 0000000..87c5246 --- /dev/null +++ b/src/lib/grids/gaussianReducedWorkerClient.ts @@ -0,0 +1,74 @@ +import type { + TGaussianReducedWorkerMetadata, + TGaussianReducedWorkerRequest, +} from "./gaussianReducedWorkerProtocol.ts"; +import { + copyGridWorkerArray, + createGridGeometryWorkerClient, +} from "./gridGeometryWorkerClient.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import type { TGridGeometryBatch } from "./gridWorkerTypes.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +type TGaussianReducedBuildRequest = { + latitudes: Float64Array; + longitudes: Float64Array; + data: Float32Array; + batchSize: number; + epsilon: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +const client = createGridGeometryWorkerClient< + TGaussianReducedWorkerRequest, + TGaussianReducedWorkerMetadata +>( + () => + new Worker(new URL("./gaussianReduced.worker.ts", import.meta.url), { + type: "module", + }) +); + +export function buildGaussianReducedGrid( + request: TGaussianReducedBuildRequest, + callbacks: { + onMetadata: (totalBatches: number) => void; + onBatch: (batch: TGridGeometryBatch) => void; + } +) { + return client + .build( + (requestId) => { + const latitudes = copyGridWorkerArray(request.latitudes); + const longitudes = copyGridWorkerArray(request.longitudes); + const data = copyGridWorkerArray(request.data); + return { + message: { + ...request, + requestId, + type: GridGeometryWorkerMessageType.BUILD, + latitudes, + longitudes, + data, + projectionCenter: { + lat: request.projectionCenter.lat, + lon: request.projectionCenter.lon, + }, + }, + transfer: [latitudes.buffer, longitudes.buffer, data.buffer], + }; + }, + { + onMetadata: ({ totalBatches }) => callbacks.onMetadata(totalBatches), + onBatch: callbacks.onBatch, + } + ) + .then(({ hoverIndexData }) => hoverIndexData); +} + +export const terminateGaussianReducedWorker = client.terminate; diff --git a/src/lib/grids/gaussianReducedWorkerProtocol.ts b/src/lib/grids/gaussianReducedWorkerProtocol.ts new file mode 100644 index 0000000..418a688 --- /dev/null +++ b/src/lib/grids/gaussianReducedWorkerProtocol.ts @@ -0,0 +1,24 @@ +import type { TGridGeometryWorkerResponse } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +export type TGaussianReducedWorkerRequest = { + requestId: number; + type: "build"; + latitudes: Float64Array; + longitudes: Float64Array; + data: Float32Array; + batchSize: number; + epsilon: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +export type TGaussianReducedWorkerMetadata = { + totalBatches: number; +}; + +export type TGaussianReducedWorkerResponse = + TGridGeometryWorkerResponse; diff --git a/src/lib/grids/gridData.worker.ts b/src/lib/grids/gridData.worker.ts new file mode 100644 index 0000000..af94fb7 --- /dev/null +++ b/src/lib/grids/gridData.worker.ts @@ -0,0 +1,42 @@ +/// + +import { ZarrDataManager } from "@/lib/data/ZarrDataManager.ts"; +import { + GridDataWorkerMessageType, + type TGridDataWorkerRequest, + type TGridDataWorkerResponse, +} from "@/lib/grids/gridDataWorkerProtocol.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; + +workerScope.onmessage = async (event: MessageEvent) => { + const { requestId, source, variable, format, selection } = event.data; + try { + const array = await ZarrDataManager.getVariableInfo( + source, + variable, + format + ); + const chunk = await ZarrDataManager.getVariableDataFromArray( + array, + selection + ); + const response: TGridDataWorkerResponse = { + requestId, + type: GridDataWorkerMessageType.RESULT, + data: chunk.data, + }; + const transfer = + ArrayBuffer.isView(chunk.data) && chunk.data.buffer instanceof ArrayBuffer + ? [chunk.data.buffer] + : []; + workerScope.postMessage(response, transfer); + } catch (error) { + const response: TGridDataWorkerResponse = { + requestId, + type: GridDataWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }; + workerScope.postMessage(response); + } +}; diff --git a/src/lib/grids/gridDataWorkerClient.ts b/src/lib/grids/gridDataWorkerClient.ts new file mode 100644 index 0000000..79c19f3 --- /dev/null +++ b/src/lib/grids/gridDataWorkerClient.ts @@ -0,0 +1,93 @@ +import type * as zarr from "zarrita"; + +import { + GridDataWorkerMessageType, + type TGridDataWorkerRequest, + type TGridDataWorkerResponse, +} from "./gridDataWorkerProtocol.ts"; + +import type { TDataSource, TZarrFormat } from "@/lib/types/GlobeTypes.ts"; + +export type TGridDataRequest = { + source: Pick; + variable: string; + format: TZarrFormat; + selection: (number | null | zarr.Slice)[]; +}; + +type TPendingRequest = { + resolve: (data: zarr.TypedArray) => void; + reject: (reason?: unknown) => void; +}; + +let worker: Worker | null = null; +let nextRequestId = 0; +const pendingRequests = new Map(); + +function rejectPendingRequests(error: Error) { + for (const pending of pendingRequests.values()) { + pending.reject(error); + } + pendingRequests.clear(); +} + +function handleWorkerMessage(message: TGridDataWorkerResponse) { + const pending = pendingRequests.get(message.requestId); + if (!pending) { + return; + } + pendingRequests.delete(message.requestId); + if (message.type === GridDataWorkerMessageType.ERROR) { + pending.reject(new Error(message.message)); + return; + } + pending.resolve(message.data); +} + +function getWorker() { + if (worker) { + return worker; + } + worker = new Worker(new URL("./gridData.worker.ts", import.meta.url), { + type: "module", + }); + worker.onmessage = (event: MessageEvent) => { + handleWorkerMessage(event.data); + }; + worker.onerror = (event) => { + rejectPendingRequests(new Error(event.message)); + worker?.terminate(); + worker = null; + }; + return worker; +} + +export function getGridVariableData(request: TGridDataRequest) { + const requestId = ++nextRequestId; + const message: TGridDataWorkerRequest = { + requestId, + type: GridDataWorkerMessageType.GET_DATA, + source: { + store: request.source.store, + dataset: request.source.dataset, + }, + variable: request.variable, + format: request.format, + selection: request.selection.map((selection) => + typeof selection === "object" && selection !== null + ? { ...selection } + : selection + ), + }; + + return new Promise>((resolve, reject) => { + pendingRequests.set(requestId, { resolve, reject }); + getWorker().postMessage(message); + }); +} + +export function terminateGridDataWorker() { + worker?.terminate(); + worker = null; + rejectPendingRequests(new Error("Grid data worker terminated.")); +} diff --git a/src/lib/grids/gridDataWorkerProtocol.ts b/src/lib/grids/gridDataWorkerProtocol.ts new file mode 100644 index 0000000..b93daa3 --- /dev/null +++ b/src/lib/grids/gridDataWorkerProtocol.ts @@ -0,0 +1,36 @@ +import type * as zarr from "zarrita"; + +import type { TDataSource, TZarrFormat } from "@/lib/types/GlobeTypes.ts"; + +export const GridDataWorkerMessageType = { + GET_DATA: "getData", + RESULT: "result", + ERROR: "error", +} as const; + +type TGridDataWorkerMessageType = + (typeof GridDataWorkerMessageType)[keyof typeof GridDataWorkerMessageType]; + +export type TGridDataWorkerRequest = { + requestId: number; + type: typeof GridDataWorkerMessageType.GET_DATA; + source: Pick; + variable: string; + format: TZarrFormat; + selection: (number | null | zarr.Slice)[]; +}; + +type TGridDataWorkerResponseBase = { + requestId: number; + type: TGridDataWorkerMessageType; +}; + +export type TGridDataWorkerResponse = + | (TGridDataWorkerResponseBase & { + type: typeof GridDataWorkerMessageType.RESULT; + data: zarr.TypedArray; + }) + | (TGridDataWorkerResponseBase & { + type: typeof GridDataWorkerMessageType.ERROR; + message: string; + }); diff --git a/src/lib/grids/gridGeometryWorkerClient.ts b/src/lib/grids/gridGeometryWorkerClient.ts new file mode 100644 index 0000000..1915053 --- /dev/null +++ b/src/lib/grids/gridGeometryWorkerClient.ts @@ -0,0 +1,151 @@ +import { + GridGeometryWorkerMessageType, + type TGridGeometryWorkerResponse, +} from "./gridGeometryWorkerProtocol.ts"; +import type { + TGridGeometryBatch, + TSerializedGeoSampleIndexData, + TGridWorkerBatch, +} from "./gridWorkerTypes.ts"; + +type TGridGeometryBuildCallbacks = { + onMetadata: (metadata: TMetadata) => void; + onBatch: (batch: TBatch) => void; +}; + +type TWorkerRequest = { + message: TRequest; + transfer: Transferable[]; +}; + +type TBuildResult = { + metadata: TMetadata; + hoverIndexData: TSerializedGeoSampleIndexData; +}; + +type TActiveBuild = { + requestId: number; + worker: Worker; + callbacks: TGridGeometryBuildCallbacks; + metadata: TMetadata | null; + hoverIndexData: TSerializedGeoSampleIndexData | null; + resolve: (value: TBuildResult) => void; + reject: (reason?: unknown) => void; +}; + +export function copyGridWorkerArray( + array: TArray +): TArray { + return array.slice() as TArray; +} + +/* eslint-disable-next-line max-lines-per-function */ +export function createGridGeometryWorkerClient< + TRequest, + TMetadata, + TBatch extends TGridWorkerBatch = TGridGeometryBatch, +>(createWorker: () => Worker, keepWorkerAlive = false) { + let activeBuild: TActiveBuild | null = null; + let nextRequestId = 0; + let worker: Worker | null = null; + + function terminate() { + if (activeBuild) { + activeBuild.reject(new Error("Grid geometry worker terminated.")); + } + activeBuild = null; + worker?.terminate(); + worker = null; + } + + function fail(error: Error) { + if (!activeBuild) { + return; + } + const { worker: activeWorker, reject } = activeBuild; + activeBuild = null; + activeWorker.terminate(); + worker = null; + reject(error); + } + + function finish() { + if (!activeBuild?.metadata || !activeBuild.hoverIndexData) { + fail(new Error("Grid geometry worker returned incomplete results.")); + return; + } + const { + worker: activeWorker, + metadata, + hoverIndexData, + resolve, + } = activeBuild; + activeBuild = null; + if (!keepWorkerAlive) { + activeWorker.terminate(); + worker = null; + } + resolve({ metadata, hoverIndexData }); + } + + function handleMessage( + message: TGridGeometryWorkerResponse + ) { + if (activeBuild?.requestId !== message.requestId) { + return; + } + if (message.type === GridGeometryWorkerMessageType.METADATA) { + activeBuild.metadata = message.metadata; + activeBuild.callbacks.onMetadata(message.metadata); + } else if (message.type === GridGeometryWorkerMessageType.BATCH) { + activeBuild.callbacks.onBatch(message.batch); + } else if (message.type === GridGeometryWorkerMessageType.HOVER_INDEX) { + activeBuild.hoverIndexData = message.hoverIndexData; + } else if (message.type === GridGeometryWorkerMessageType.ERROR) { + fail(new Error(message.message)); + } else { + finish(); + } + } + + function build( + createRequest: (requestId: number) => TWorkerRequest, + callbacks: TGridGeometryBuildCallbacks + ) { + if (activeBuild) { + terminate(); + } + const requestId = ++nextRequestId; + worker ??= createWorker(); + const activeWorker = worker; + return new Promise>((resolve, reject) => { + activeBuild = { + requestId, + worker: activeWorker, + callbacks, + metadata: null, + hoverIndexData: null, + resolve, + reject, + }; + activeWorker.onmessage = ( + event: MessageEvent> + ) => { + try { + handleMessage(event.data); + } catch (error) { + fail(error instanceof Error ? error : new Error(String(error))); + } + }; + activeWorker.onerror = (event) => { + if (activeBuild?.requestId === requestId) { + fail(new Error(event.message)); + } + }; + const request = createRequest(requestId); + activeWorker.postMessage(request.message, request.transfer); + }); + } + + return { build, terminate }; +} diff --git a/src/lib/grids/gridGeometryWorkerProtocol.ts b/src/lib/grids/gridGeometryWorkerProtocol.ts new file mode 100644 index 0000000..2185d92 --- /dev/null +++ b/src/lib/grids/gridGeometryWorkerProtocol.ts @@ -0,0 +1,46 @@ +import type { + TGridGeometryBatch, + TSerializedGeoSampleIndexData, + TGridWorkerBatch, +} from "./gridWorkerTypes.ts"; + +export const GridGeometryWorkerMessageType = { + BUILD: "build", + METADATA: "metadata", + BATCH: "batch", + HOVER_INDEX: "hoverIndex", + DONE: "done", + ERROR: "error", +} as const; + +type TGridGeometryWorkerMessageType = + (typeof GridGeometryWorkerMessageType)[keyof typeof GridGeometryWorkerMessageType]; + +type TGridGeometryWorkerResponseBase = { + requestId: number; + type: TGridGeometryWorkerMessageType; +}; + +export type TGridGeometryWorkerResponse< + TMetadata, + TBatch extends TGridWorkerBatch = TGridGeometryBatch, +> = + | (TGridGeometryWorkerResponseBase & { + type: typeof GridGeometryWorkerMessageType.METADATA; + metadata: TMetadata; + }) + | (TGridGeometryWorkerResponseBase & { + type: typeof GridGeometryWorkerMessageType.BATCH; + batch: TBatch; + }) + | (TGridGeometryWorkerResponseBase & { + type: typeof GridGeometryWorkerMessageType.HOVER_INDEX; + hoverIndexData: TSerializedGeoSampleIndexData; + }) + | (TGridGeometryWorkerResponseBase & { + type: typeof GridGeometryWorkerMessageType.DONE; + }) + | (TGridGeometryWorkerResponseBase & { + type: typeof GridGeometryWorkerMessageType.ERROR; + message: string; + }); diff --git a/src/lib/grids/gridGeometryWorkerUtils.ts b/src/lib/grids/gridGeometryWorkerUtils.ts new file mode 100644 index 0000000..fc09ba9 --- /dev/null +++ b/src/lib/grids/gridGeometryWorkerUtils.ts @@ -0,0 +1,115 @@ +import type { TGridGeometryWorkerResponse } from "./gridGeometryWorkerProtocol.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import type { + TGridDataValueBatch, + TGridGeometryBatch, + TGridPointBatch, + TGridPositionBatch, + TSerializedGeoSampleIndexData, + TGridWorkerBatch, +} from "./gridWorkerTypes.ts"; + +export function postGridGeometryResponse< + TMetadata, + TBatch extends TGridWorkerBatch = TGridGeometryBatch, +>( + workerScope: DedicatedWorkerGlobalScope, + response: TGridGeometryWorkerResponse, + transfer: Transferable[] = [] +) { + workerScope.postMessage(response, transfer); +} + +export function postGridGeometryHoverIndex( + workerScope: DedicatedWorkerGlobalScope, + requestId: number, + hoverIndexData: TSerializedGeoSampleIndexData +) { + postGridGeometryResponse( + workerScope, + { + requestId, + type: GridGeometryWorkerMessageType.HOVER_INDEX, + hoverIndexData, + }, + [ + hoverIndexData.indexData, + hoverIndexData.latitudes.buffer, + hoverIndexData.longitudes.buffer, + hoverIndexData.values.buffer, + ] + ); +} + +export function postGridGeometryBatch( + workerScope: DedicatedWorkerGlobalScope, + requestId: number, + batch: TGridGeometryBatch +) { + postGridGeometryResponse( + workerScope, + { + requestId, + type: GridGeometryWorkerMessageType.BATCH, + batch, + }, + [ + batch.positionValues.buffer, + batch.dataValues.buffer, + batch.latLonValues.buffer, + batch.indices.buffer, + ] + ); +} + +export function postGridPointBatch( + workerScope: DedicatedWorkerGlobalScope, + requestId: number, + batch: TGridPointBatch +) { + postGridGeometryResponse( + workerScope, + { + requestId, + type: GridGeometryWorkerMessageType.BATCH, + batch, + }, + [ + batch.positionValues.buffer, + batch.dataValues.buffer, + batch.latLonValues.buffer, + ] + ); +} + +export function postGridDataValueBatch( + workerScope: DedicatedWorkerGlobalScope, + requestId: number, + batch: TGridDataValueBatch +) { + postGridGeometryResponse( + workerScope, + { + requestId, + type: GridGeometryWorkerMessageType.BATCH, + batch, + }, + [batch.dataValues.buffer] + ); +} + +export function postGridPositionBatch( + workerScope: DedicatedWorkerGlobalScope, + requestId: number, + batch: TGridPositionBatch +) { + postGridGeometryResponse( + workerScope, + { + requestId, + type: GridGeometryWorkerMessageType.BATCH, + batch, + }, + [batch.positionValues.buffer, batch.latLonValues.buffer] + ); +} diff --git a/src/lib/grids/gridWorkerCalculations.ts b/src/lib/grids/gridWorkerCalculations.ts new file mode 100644 index 0000000..465013d --- /dev/null +++ b/src/lib/grids/gridWorkerCalculations.ts @@ -0,0 +1,53 @@ +import KDBush from "kdbush"; + +import type { TSerializedGeoSampleIndexData } from "./gridWorkerTypes.ts"; + +export function shouldFlipCartesianTriangle( + x0: number, + y0: number, + z0: number, + x1: number, + y1: number, + z1: number, + x2: number, + y2: number, + z2: number +) { + const edgeAx = x1 - x0; + const edgeAy = y1 - y0; + const edgeAz = z1 - z0; + const edgeBx = x2 - x0; + const edgeBy = y2 - y0; + const edgeBz = z2 - z0; + const normalX = edgeAy * edgeBz - edgeAz * edgeBy; + const normalY = edgeAz * edgeBx - edgeAx * edgeBz; + const normalZ = edgeAx * edgeBy - edgeAy * edgeBx; + const centerX = x0 + x1 + x2; + const centerY = y0 + y1 + y2; + const centerZ = z0 + z1 + z2; + return normalX * centerX + normalY * centerY + normalZ * centerZ < 0; +} + +export function buildSerializedGeoSampleIndexData( + latitudes: Float64Array, + longitudes: Float64Array, + values: Float32Array +): TSerializedGeoSampleIndexData { + const index = new KDBush(latitudes.length); + for (let sampleIndex = 0; sampleIndex < latitudes.length; sampleIndex++) { + latitudes[sampleIndex] = Math.max( + -90, + Math.min(90, latitudes[sampleIndex]) + ); + longitudes[sampleIndex] = + (((longitudes[sampleIndex] % 360) + 540) % 360) - 180; + index.add(longitudes[sampleIndex], latitudes[sampleIndex]); + } + index.finish(); + return { + indexData: index.data as ArrayBuffer, + latitudes, + longitudes, + values, + }; +} diff --git a/src/lib/grids/gridWorkerTypes.ts b/src/lib/grids/gridWorkerTypes.ts new file mode 100644 index 0000000..fafd9dc --- /dev/null +++ b/src/lib/grids/gridWorkerTypes.ts @@ -0,0 +1,29 @@ +export type TGridGeometryBatch = { + batchIndex: number; + positionValues: Float32Array; + dataValues: Float32Array; + latLonValues: Float32Array; + indices: Uint32Array; +}; + +export type TGridPointBatch = Omit; + +export type TGridPositionBatch = Omit; + +export type TGridDataValueBatch = { + batchIndex: number; + dataValues: Float32Array; +}; + +export type TGridWorkerBatch = + | TGridGeometryBatch + | TGridPointBatch + | TGridPositionBatch + | TGridDataValueBatch; + +export type TSerializedGeoSampleIndexData = { + indexData: ArrayBuffer; + latitudes: Float64Array; + longitudes: Float64Array; + values: Float32Array; +}; diff --git a/src/lib/grids/irregular.worker.ts b/src/lib/grids/irregular.worker.ts new file mode 100644 index 0000000..fad9331 --- /dev/null +++ b/src/lib/grids/irregular.worker.ts @@ -0,0 +1,78 @@ +/// + +import { + buildIrregularBatch, + buildIrregularGridData, + buildIrregularHoverIndexData, + getIrregularBatchCount, +} from "./irregularCalculations.ts"; +import type { + TIrregularWorkerRequest, + TIrregularWorkerResponse, +} from "./irregularWorkerProtocol.ts"; + +import { GridGeometryWorkerMessageType } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import { + postGridGeometryHoverIndex, + postGridGeometryResponse, + postGridPointBatch, +} from "@/lib/grids/gridGeometryWorkerUtils.ts"; +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; + +function postResponse( + response: TIrregularWorkerResponse, + transfer: Transferable[] = [] +) { + postGridGeometryResponse(workerScope, response, transfer); +} + +function buildGrid(request: TIrregularWorkerRequest) { + const projection = new ProjectionHelper( + request.projectionType, + request.projectionCenter + ); + const grid = buildIrregularGridData( + request.latitudes, + request.longitudes, + request.latitudeShape, + request.longitudeShape, + request.data, + projection + ); + const totalBatches = getIrregularBatchCount(grid, request.batchSize); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { totalBatches, estimatedSpacing: grid.estimatedSpacing }, + }); + postGridGeometryHoverIndex( + workerScope, + request.requestId, + buildIrregularHoverIndexData(grid) + ); + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + postGridPointBatch( + workerScope, + request.requestId, + buildIrregularBatch(grid, batchIndex, request.batchSize) + ); + } + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.DONE, + }); +} + +workerScope.onmessage = (event: MessageEvent) => { + try { + buildGrid(event.data); + } catch (error) { + postResponse({ + requestId: event.data.requestId, + type: GridGeometryWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }); + } +}; diff --git a/src/lib/grids/irregularCalculations.ts b/src/lib/grids/irregularCalculations.ts new file mode 100644 index 0000000..9f2dbf2 --- /dev/null +++ b/src/lib/grids/irregularCalculations.ts @@ -0,0 +1,129 @@ +import { buildSerializedGeoSampleIndexData } from "./gridWorkerCalculations.ts"; +import type { + TGridPointBatch, + TSerializedGeoSampleIndexData, +} from "./gridWorkerTypes.ts"; + +import { reconcileCoordinateArrays } from "@/lib/data/irregularGridHelpers.ts"; +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +export type TIrregularGridData = { + positions: Float32Array; + latLonValues: Float32Array; + data: Float32Array; + latitudes: Float32Array; + longitudes: Float32Array; + estimatedSpacing: number; +}; + +export function estimateIrregularAverageSpacing( + positions: Float32Array, + sampleSize = 5000 +) { + const pointCount = positions.length / 3; + const stride = Math.max(1, Math.floor(pointCount / sampleSize)); + let totalDistance = 0; + let totalWeight = 0; + + for (let pointIndex = 0; pointIndex < pointCount; pointIndex += stride) { + const offset = pointIndex * 3; + const x = positions[offset]; + const y = positions[offset + 1]; + const z = positions[offset + 2]; + const weight = Math.cos(Math.asin(z)); + if (weight < 0.1) { + continue; + } + let nearestDistance = Number.POSITIVE_INFINITY; + const start = Math.max(0, pointIndex - 50); + const end = Math.min(pointCount, pointIndex + 50); + for (let neighborIndex = start; neighborIndex < end; neighborIndex++) { + if (pointIndex === neighborIndex) { + continue; + } + const neighborOffset = neighborIndex * 3; + const deltaX = positions[neighborOffset] - x; + const deltaY = positions[neighborOffset + 1] - y; + const deltaZ = positions[neighborOffset + 2] - z; + const distance = Math.sqrt( + deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ + ); + if (distance > 0 && distance < nearestDistance) { + nearestDistance = distance; + } + } + if (nearestDistance !== Number.POSITIVE_INFINITY) { + totalDistance += nearestDistance * weight; + totalWeight += weight; + } + } + return totalWeight > 0 ? totalDistance / totalWeight : 0.01; +} + +export function buildIrregularGridData( + latitudes: Float32Array, + longitudes: Float32Array, + latitudeShape: number[], + longitudeShape: number[], + data: Float32Array, + projection: ProjectionHelper +): TIrregularGridData { + const coordinates = reconcileCoordinateArrays( + latitudes, + longitudes, + latitudeShape, + longitudeShape, + data.length + ); + const positions = new Float32Array(data.length * 3); + const latLonValues = new Float32Array(data.length * 2); + for (let index = 0; index < data.length; index++) { + projection.projectLatLonToArrays( + coordinates.latitudes[index], + coordinates.longitudes[index], + positions, + index * 3, + latLonValues, + index * 2 + ); + } + return { + positions, + latLonValues, + data, + ...coordinates, + estimatedSpacing: estimateIrregularAverageSpacing(positions), + }; +} + +export function getIrregularBatchCount( + grid: TIrregularGridData, + batchSize: number +) { + return Math.ceil(grid.data.length / batchSize); +} + +export function buildIrregularBatch( + grid: TIrregularGridData, + batchIndex: number, + batchSize: number +): TGridPointBatch { + const start = batchIndex * batchSize; + const end = Math.min(start + batchSize, grid.data.length); + return { + batchIndex, + positionValues: grid.positions.slice(start * 3, end * 3), + latLonValues: grid.latLonValues.slice(start * 2, end * 2), + dataValues: grid.data.slice(start, end), + }; +} + +export function buildIrregularHoverIndexData( + grid: TIrregularGridData +): TSerializedGeoSampleIndexData { + return buildSerializedGeoSampleIndexData( + Float64Array.from(grid.latitudes), + Float64Array.from(grid.longitudes), + grid.data.slice() + ); +} diff --git a/src/lib/grids/irregularDelaunay.worker.ts b/src/lib/grids/irregularDelaunay.worker.ts new file mode 100644 index 0000000..f4d11af --- /dev/null +++ b/src/lib/grids/irregularDelaunay.worker.ts @@ -0,0 +1,114 @@ +/// + +import { + buildIrregularDelaunayDataBatch, + buildIrregularDelaunayGeometryBatch, + buildIrregularDelaunayGrid, + buildIrregularDelaunayHoverIndexData, + getIrregularDelaunayBatchCount, + type TIrregularDelaunayGrid, +} from "./irregularDelaunayCalculations.ts"; +import type { + TIrregularDelaunayWorkerRequest, + TIrregularDelaunayWorkerResponse, +} from "./irregularDelaunayWorkerProtocol.ts"; + +import { GridGeometryWorkerMessageType } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import { + postGridDataValueBatch, + postGridGeometryHoverIndex, + postGridGeometryResponse, + postGridPointBatch, +} from "@/lib/grids/gridGeometryWorkerUtils.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; +let cachedGrid: TIrregularDelaunayGrid | null = null; +let cachedPointCount = 0; + +function postResponse( + response: TIrregularDelaunayWorkerResponse, + transfer: Transferable[] = [] +) { + postGridGeometryResponse(workerScope, response, transfer); +} + +function postBatches( + request: TIrregularDelaunayWorkerRequest, + grid: TIrregularDelaunayGrid, + totalBatches: number, + rebuildGeometry: boolean +) { + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + if (rebuildGeometry) { + postGridPointBatch( + workerScope, + request.requestId, + buildIrregularDelaunayGeometryBatch( + grid, + request.data, + batchIndex, + request.batchSize + ) + ); + } else { + postGridDataValueBatch( + workerScope, + request.requestId, + buildIrregularDelaunayDataBatch( + grid, + request.data, + batchIndex, + request.batchSize + ) + ); + } + } +} + +function buildGrid(request: TIrregularDelaunayWorkerRequest) { + const rebuildGeometry = + request.forceGeometryRebuild || + !cachedGrid || + cachedPointCount !== request.data.length; + if (rebuildGeometry) { + cachedGrid = buildIrregularDelaunayGrid( + request.latitudes, + request.longitudes, + request.latitudeShape, + request.longitudeShape, + request.data.length + ); + cachedPointCount = request.data.length; + } + const grid = cachedGrid!; + const totalBatches = getIrregularDelaunayBatchCount(grid, request.batchSize); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { totalBatches, rebuildGeometry }, + }); + postGridGeometryHoverIndex( + workerScope, + request.requestId, + buildIrregularDelaunayHoverIndexData(grid, request.data) + ); + postBatches(request, grid, totalBatches, rebuildGeometry); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.DONE, + }); +} + +workerScope.onmessage = ( + event: MessageEvent +) => { + try { + buildGrid(event.data); + } catch (error) { + postResponse({ + requestId: event.data.requestId, + type: GridGeometryWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }); + } +}; diff --git a/src/lib/grids/irregularDelaunayCalculations.ts b/src/lib/grids/irregularDelaunayCalculations.ts new file mode 100644 index 0000000..052561e --- /dev/null +++ b/src/lib/grids/irregularDelaunayCalculations.ts @@ -0,0 +1,313 @@ +import { Delaunay } from "d3-delaunay"; + +import { + buildSerializedGeoSampleIndexData, + shouldFlipCartesianTriangle, +} from "./gridWorkerCalculations.ts"; +import type { + TGridDataValueBatch, + TGridPointBatch, + TSerializedGeoSampleIndexData, +} from "./gridWorkerTypes.ts"; + +import { reconcileCoordinateArrays } from "@/lib/data/irregularGridHelpers.ts"; +import { + PROJECTION_TYPES, + ProjectionHelper, +} from "@/lib/projection/projectionUtils.ts"; + +type TExtendedPointSet = { + points: [number, number][]; + indexMap: number[]; +}; + +export type TIrregularDelaunayGrid = { + triangleIndices: Uint32Array; + positions: Float32Array; + latLonValues: Float32Array; + latitudes: Float32Array; + longitudes: Float32Array; +}; + +function createExtendedPointSet( + latitudes: Float32Array, + longitudes: Float32Array, + wrapThreshold: number +): TExtendedPointSet { + const points: [number, number][] = []; + const indexMap: number[] = []; + for (let index = 0; index < latitudes.length; index++) { + const longitude = longitudes[index]; + const latitude = latitudes[index]; + points.push([longitude, latitude]); + indexMap.push(index); + if (longitude < -wrapThreshold) { + points.push([longitude + 360, latitude]); + indexMap.push(index); + } + if (longitude > wrapThreshold) { + points.push([longitude - 360, latitude]); + indexMap.push(index); + } + } + return { points, indexMap }; +} + +function isValidTriangle( + points: [number, number][], + pointIndices: [number, number, number], + maxEdgeLength: number +) { + const [point0, point1, point2] = pointIndices.map( + (pointIndex) => points[pointIndex] + ); + if ( + (point0[0] < -180 || point0[0] > 180) && + (point1[0] < -180 || point1[0] > 180) && + (point2[0] < -180 || point2[0] > 180) + ) { + return false; + } + const edges = [ + [point0, point1], + [point1, point2], + [point2, point0], + ] as const; + for (const [[longitude0, latitude0], [longitude1, latitude1]] of edges) { + const longitudeDelta = longitude1 - longitude0; + const latitudeDelta = latitude1 - latitude0; + if ( + Math.sqrt( + longitudeDelta * longitudeDelta + latitudeDelta * latitudeDelta + ) > maxEdgeLength + ) { + return false; + } + } + return true; +} + +export function computeDelaunayTriangulation( + latitudes: Float32Array, + longitudes: Float32Array +): Uint32Array { + const { points, indexMap } = createExtendedPointSet( + latitudes, + longitudes, + 150 + ); + const triangles = Delaunay.from(points).triangles; + const validTriangles: number[] = []; + for ( + let triangleOffset = 0; + triangleOffset < triangles.length; + triangleOffset += 3 + ) { + const pointIndices: [number, number, number] = [ + triangles[triangleOffset], + triangles[triangleOffset + 1], + triangles[triangleOffset + 2], + ]; + if (isValidTriangle(points, pointIndices, 30)) { + validTriangles.push( + indexMap[pointIndices[0]], + indexMap[pointIndices[1]], + indexMap[pointIndices[2]] + ); + } + } + return new Uint32Array(validTriangles); +} + +export function computeTriangleAverage( + value0: number, + value1: number, + value2: number +) { + const valid0 = Number.isFinite(value0); + const valid1 = Number.isFinite(value1); + const valid2 = Number.isFinite(value2); + if (!valid0 && !valid1 && !valid2) { + return NaN; + } + let sum = 0; + let count = 0; + if (valid0) { + sum += value0; + count++; + } + if (valid1) { + sum += value1; + count++; + } + if (valid2) { + sum += value2; + count++; + } + return sum / count; +} + +function shouldFlipTriangle( + positions: Float32Array, + index0: number, + index1: number, + index2: number +) { + const offset0 = index0 * 3; + const offset1 = index1 * 3; + const offset2 = index2 * 3; + const x0 = positions[offset0]; + const y0 = positions[offset0 + 1]; + const z0 = positions[offset0 + 2]; + return shouldFlipCartesianTriangle( + x0, + y0, + z0, + positions[offset1], + positions[offset1 + 1], + positions[offset1 + 2], + positions[offset2], + positions[offset2 + 1], + positions[offset2 + 2] + ); +} + +export function buildIrregularDelaunayGrid( + latitudes: Float32Array, + longitudes: Float32Array, + latitudeShape: number[], + longitudeShape: number[], + dataLength: number +): TIrregularDelaunayGrid { + const coordinates = reconcileCoordinateArrays( + latitudes, + longitudes, + latitudeShape, + longitudeShape, + dataLength + ); + const positions = new Float32Array(dataLength * 3); + const latLonValues = new Float32Array(dataLength * 2); + const projection = new ProjectionHelper( + PROJECTION_TYPES.NEARSIDE_PERSPECTIVE, + { lat: 0, lon: 0 } + ); + for (let index = 0; index < dataLength; index++) { + projection.projectLatLonToArrays( + coordinates.latitudes[index], + coordinates.longitudes[index], + positions, + index * 3, + latLonValues, + index * 2 + ); + } + return { + ...coordinates, + positions, + latLonValues, + triangleIndices: computeDelaunayTriangulation( + coordinates.latitudes, + coordinates.longitudes + ), + }; +} + +export function getIrregularDelaunayBatchCount( + grid: TIrregularDelaunayGrid, + batchSize: number +) { + return Math.ceil(grid.triangleIndices.length / 3 / batchSize); +} + +function getTriangleAverage( + grid: TIrregularDelaunayGrid, + data: Float32Array, + triangleIndex: number +) { + const offset = triangleIndex * 3; + const index0 = grid.triangleIndices[offset]; + const index1 = grid.triangleIndices[offset + 1]; + const index2 = grid.triangleIndices[offset + 2]; + return computeTriangleAverage(data[index0], data[index1], data[index2]); +} + +export function buildIrregularDelaunayGeometryBatch( + grid: TIrregularDelaunayGrid, + data: Float32Array, + batchIndex: number, + batchSize: number +): TGridPointBatch { + const start = batchIndex * batchSize; + const end = Math.min(start + batchSize, grid.triangleIndices.length / 3); + const vertexCount = (end - start) * 3; + const positionValues = new Float32Array(vertexCount * 3); + const latLonValues = new Float32Array(vertexCount * 2); + const dataValues = new Float32Array(vertexCount); + for (let triangleIndex = start; triangleIndex < end; triangleIndex++) { + const triangleOffset = triangleIndex * 3; + const index0 = grid.triangleIndices[triangleOffset]; + let index1 = grid.triangleIndices[triangleOffset + 1]; + let index2 = grid.triangleIndices[triangleOffset + 2]; + const average = computeTriangleAverage( + data[index0], + data[index1], + data[index2] + ); + if (shouldFlipTriangle(grid.positions, index0, index1, index2)) { + [index1, index2] = [index2, index1]; + } + const batchTriangleIndex = triangleIndex - start; + for (let vertexIndex = 0; vertexIndex < 3; vertexIndex++) { + const sourceIndex = + vertexIndex === 0 ? index0 : vertexIndex === 1 ? index1 : index2; + const destinationIndex = batchTriangleIndex * 3 + vertexIndex; + const sourcePositionOffset = sourceIndex * 3; + const destinationPositionOffset = destinationIndex * 3; + positionValues[destinationPositionOffset] = + grid.positions[sourcePositionOffset]; + positionValues[destinationPositionOffset + 1] = + grid.positions[sourcePositionOffset + 1]; + positionValues[destinationPositionOffset + 2] = + grid.positions[sourcePositionOffset + 2]; + const sourceLatLonOffset = sourceIndex * 2; + const destinationLatLonOffset = destinationIndex * 2; + latLonValues[destinationLatLonOffset] = + grid.latLonValues[sourceLatLonOffset]; + latLonValues[destinationLatLonOffset + 1] = + grid.latLonValues[sourceLatLonOffset + 1]; + dataValues[destinationIndex] = average; + } + } + return { batchIndex, positionValues, latLonValues, dataValues }; +} + +export function buildIrregularDelaunayDataBatch( + grid: TIrregularDelaunayGrid, + data: Float32Array, + batchIndex: number, + batchSize: number +): TGridDataValueBatch { + const start = batchIndex * batchSize; + const end = Math.min(start + batchSize, grid.triangleIndices.length / 3); + const dataValues = new Float32Array((end - start) * 3); + for (let triangleIndex = start; triangleIndex < end; triangleIndex++) { + const average = getTriangleAverage(grid, data, triangleIndex); + const batchOffset = (triangleIndex - start) * 3; + dataValues[batchOffset] = average; + dataValues[batchOffset + 1] = average; + dataValues[batchOffset + 2] = average; + } + return { batchIndex, dataValues }; +} + +export function buildIrregularDelaunayHoverIndexData( + grid: TIrregularDelaunayGrid, + data: Float32Array +): TSerializedGeoSampleIndexData { + return buildSerializedGeoSampleIndexData( + Float64Array.from(grid.latitudes), + Float64Array.from(grid.longitudes), + data.slice() + ); +} diff --git a/src/lib/grids/irregularDelaunayWorkerClient.ts b/src/lib/grids/irregularDelaunayWorkerClient.ts new file mode 100644 index 0000000..6db6fea --- /dev/null +++ b/src/lib/grids/irregularDelaunayWorkerClient.ts @@ -0,0 +1,61 @@ +import { + copyGridWorkerArray, + createGridGeometryWorkerClient, +} from "./gridGeometryWorkerClient.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import type { + TIrregularDelaunayBatch, + TIrregularDelaunayWorkerMetadata, + TIrregularDelaunayWorkerRequest, +} from "./irregularDelaunayWorkerProtocol.ts"; + +type TIrregularDelaunayBuildRequest = { + latitudes: Float32Array; + longitudes: Float32Array; + latitudeShape: number[]; + longitudeShape: number[]; + data: Float32Array; + batchSize: number; + forceGeometryRebuild: boolean; +}; + +const client = createGridGeometryWorkerClient< + TIrregularDelaunayWorkerRequest, + TIrregularDelaunayWorkerMetadata, + TIrregularDelaunayBatch +>( + () => + new Worker(new URL("./irregularDelaunay.worker.ts", import.meta.url), { + type: "module", + }), + true +); + +export function buildIrregularDelaunayGrid( + request: TIrregularDelaunayBuildRequest, + callbacks: { + onMetadata: (metadata: TIrregularDelaunayWorkerMetadata) => void; + onBatch: (batch: TIrregularDelaunayBatch) => void; + } +) { + return client.build((requestId) => { + const latitudes = copyGridWorkerArray(request.latitudes); + const longitudes = copyGridWorkerArray(request.longitudes); + const data = copyGridWorkerArray(request.data); + return { + message: { + ...request, + requestId, + type: GridGeometryWorkerMessageType.BUILD, + latitudes, + longitudes, + data, + latitudeShape: [...request.latitudeShape], + longitudeShape: [...request.longitudeShape], + }, + transfer: [latitudes.buffer, longitudes.buffer, data.buffer], + }; + }, callbacks); +} + +export const terminateIrregularDelaunayWorker = client.terminate; diff --git a/src/lib/grids/irregularDelaunayWorkerProtocol.ts b/src/lib/grids/irregularDelaunayWorkerProtocol.ts new file mode 100644 index 0000000..f60691a --- /dev/null +++ b/src/lib/grids/irregularDelaunayWorkerProtocol.ts @@ -0,0 +1,29 @@ +import type { TGridGeometryWorkerResponse } from "./gridGeometryWorkerProtocol.ts"; +import type { + TGridDataValueBatch, + TGridPointBatch, +} from "./gridWorkerTypes.ts"; + +export type TIrregularDelaunayBatch = TGridPointBatch | TGridDataValueBatch; + +export type TIrregularDelaunayWorkerRequest = { + requestId: number; + type: "build"; + latitudes: Float32Array; + longitudes: Float32Array; + latitudeShape: number[]; + longitudeShape: number[]; + data: Float32Array; + batchSize: number; + forceGeometryRebuild: boolean; +}; + +export type TIrregularDelaunayWorkerMetadata = { + totalBatches: number; + rebuildGeometry: boolean; +}; + +export type TIrregularDelaunayWorkerResponse = TGridGeometryWorkerResponse< + TIrregularDelaunayWorkerMetadata, + TIrregularDelaunayBatch +>; diff --git a/src/lib/grids/irregularWorkerClient.ts b/src/lib/grids/irregularWorkerClient.ts new file mode 100644 index 0000000..02794cf --- /dev/null +++ b/src/lib/grids/irregularWorkerClient.ts @@ -0,0 +1,70 @@ +import { + copyGridWorkerArray, + createGridGeometryWorkerClient, +} from "./gridGeometryWorkerClient.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import type { TGridPointBatch } from "./gridWorkerTypes.ts"; +import type { + TIrregularWorkerMetadata, + TIrregularWorkerRequest, +} from "./irregularWorkerProtocol.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +type TIrregularBuildRequest = { + latitudes: Float32Array; + longitudes: Float32Array; + latitudeShape: number[]; + longitudeShape: number[]; + data: Float32Array; + batchSize: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +const client = createGridGeometryWorkerClient< + TIrregularWorkerRequest, + TIrregularWorkerMetadata, + TGridPointBatch +>( + () => + new Worker(new URL("./irregular.worker.ts", import.meta.url), { + type: "module", + }) +); + +export function buildIrregularGrid( + request: TIrregularBuildRequest, + callbacks: { + onMetadata: (metadata: TIrregularWorkerMetadata) => void; + onBatch: (batch: TGridPointBatch) => void; + } +) { + return client.build((requestId) => { + const latitudes = copyGridWorkerArray(request.latitudes); + const longitudes = copyGridWorkerArray(request.longitudes); + const data = copyGridWorkerArray(request.data); + return { + message: { + ...request, + requestId, + type: GridGeometryWorkerMessageType.BUILD, + latitudes, + longitudes, + data, + latitudeShape: [...request.latitudeShape], + longitudeShape: [...request.longitudeShape], + projectionCenter: { + lat: request.projectionCenter.lat, + lon: request.projectionCenter.lon, + }, + }, + transfer: [latitudes.buffer, longitudes.buffer, data.buffer], + }; + }, callbacks); +} + +export const terminateIrregularWorker = client.terminate; diff --git a/src/lib/grids/irregularWorkerProtocol.ts b/src/lib/grids/irregularWorkerProtocol.ts new file mode 100644 index 0000000..b20d67d --- /dev/null +++ b/src/lib/grids/irregularWorkerProtocol.ts @@ -0,0 +1,30 @@ +import type { TGridGeometryWorkerResponse } from "./gridGeometryWorkerProtocol.ts"; +import type { TGridPointBatch } from "./gridWorkerTypes.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +export type TIrregularWorkerRequest = { + requestId: number; + type: "build"; + latitudes: Float32Array; + longitudes: Float32Array; + latitudeShape: number[]; + longitudeShape: number[]; + data: Float32Array; + batchSize: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +export type TIrregularWorkerMetadata = { + totalBatches: number; + estimatedSpacing: number; +}; + +export type TIrregularWorkerResponse = TGridGeometryWorkerResponse< + TIrregularWorkerMetadata, + TGridPointBatch +>; diff --git a/src/lib/grids/serializedGeoSampleIndex.ts b/src/lib/grids/serializedGeoSampleIndex.ts new file mode 100644 index 0000000..e927a54 --- /dev/null +++ b/src/lib/grids/serializedGeoSampleIndex.ts @@ -0,0 +1,37 @@ +import { around } from "geokdbush"; +import KDBush from "kdbush"; + +import type { TSerializedGeoSampleIndexData } from "./gridWorkerTypes.ts"; + +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +type TSerializedGeoSampleIndex = { + findNearest( + lat: number, + lon: number + ): { lat: number; lon: number; value: number } | null; +}; + +export function createSerializedGeoSampleIndex( + data: TSerializedGeoSampleIndexData +): TSerializedGeoSampleIndex { + const index = KDBush.from(data.indexData); + return { + findNearest(lat, lon) { + const nearest = around( + index, + ProjectionHelper.normalizeLongitude(lon), + Math.max(-90, Math.min(90, lat)), + 1 + )[0]; + if (nearest === undefined) { + return null; + } + return { + lat: data.latitudes[nearest], + lon: data.longitudes[nearest], + value: data.values[nearest], + }; + }, + }; +} diff --git a/src/lib/grids/triangular.worker.ts b/src/lib/grids/triangular.worker.ts new file mode 100644 index 0000000..e463044 --- /dev/null +++ b/src/lib/grids/triangular.worker.ts @@ -0,0 +1,133 @@ +/// + +import { + buildTriangularDataBatch, + buildTriangularGeometryBatch, + buildTriangularGrid, + buildTriangularHoverIndexData, + getTriangularBatchCount, + type TTriangularGrid, +} from "./triangularCalculations.ts"; +import { + TriangularWorkerOperation, + type TTriangularDataWorkerRequest, + type TTriangularGeometryWorkerRequest, + type TTriangularWorkerRequest, + type TTriangularWorkerResponse, +} from "./triangularWorkerProtocol.ts"; + +import { GridGeometryWorkerMessageType } from "@/lib/grids/gridGeometryWorkerProtocol.ts"; +import { + postGridDataValueBatch, + postGridGeometryHoverIndex, + postGridGeometryResponse, + postGridPositionBatch, +} from "@/lib/grids/gridGeometryWorkerUtils.ts"; +import { buildSerializedGeoSampleIndexData } from "@/lib/grids/gridWorkerCalculations.ts"; +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +const workerScope = self as unknown as DedicatedWorkerGlobalScope; +let cachedGrid: TTriangularGrid | null = null; + +function postResponse( + response: TTriangularWorkerResponse, + transfer: Transferable[] = [] +) { + postGridGeometryResponse(workerScope, response, transfer); +} + +function postDone(requestId: number) { + postResponse({ + requestId, + type: GridGeometryWorkerMessageType.DONE, + }); +} + +function buildGeometry(request: TTriangularGeometryWorkerRequest) { + cachedGrid = buildTriangularGrid( + request.vertexOfCell, + request.vertexX, + request.vertexY, + request.vertexZ + ); + const triangleCount = cachedGrid.vertices.length / 9; + const totalBatches = getTriangularBatchCount( + triangleCount, + request.batchSize + ); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { operation: request.operation, totalBatches }, + }); + postGridGeometryHoverIndex( + workerScope, + request.requestId, + buildSerializedGeoSampleIndexData( + new Float64Array(0), + new Float64Array(0), + new Float32Array(0) + ) + ); + const projection = new ProjectionHelper( + request.projectionType, + request.projectionCenter + ); + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + postGridPositionBatch( + workerScope, + request.requestId, + buildTriangularGeometryBatch( + cachedGrid, + batchIndex, + request.batchSize, + projection + ) + ); + } + postDone(request.requestId); +} + +function buildData(request: TTriangularDataWorkerRequest) { + if (!cachedGrid) { + throw new Error("Triangular grid geometry has not been built."); + } + const hoverIndexData = buildTriangularHoverIndexData( + cachedGrid, + request.data + ); + const totalBatches = getTriangularBatchCount( + request.data.length, + request.batchSize + ); + postResponse({ + requestId: request.requestId, + type: GridGeometryWorkerMessageType.METADATA, + metadata: { operation: request.operation, totalBatches }, + }); + postGridGeometryHoverIndex(workerScope, request.requestId, hoverIndexData); + for (let batchIndex = 0; batchIndex < totalBatches; batchIndex++) { + postGridDataValueBatch( + workerScope, + request.requestId, + buildTriangularDataBatch(request.data, batchIndex, request.batchSize) + ); + } + postDone(request.requestId); +} + +workerScope.onmessage = (event: MessageEvent) => { + try { + if (event.data.operation === TriangularWorkerOperation.GEOMETRY) { + buildGeometry(event.data); + } else { + buildData(event.data); + } + } catch (error) { + postResponse({ + requestId: event.data.requestId, + type: GridGeometryWorkerMessageType.ERROR, + message: error instanceof Error ? error.message : String(error), + }); + } +}; diff --git a/src/lib/grids/triangularCalculations.ts b/src/lib/grids/triangularCalculations.ts new file mode 100644 index 0000000..6c0e2a7 --- /dev/null +++ b/src/lib/grids/triangularCalculations.ts @@ -0,0 +1,206 @@ +import { + buildSerializedGeoSampleIndexData, + shouldFlipCartesianTriangle, +} from "./gridWorkerCalculations.ts"; +import type { + TGridDataValueBatch, + TGridPositionBatch, + TSerializedGeoSampleIndexData, +} from "./gridWorkerTypes.ts"; + +import { ProjectionHelper } from "@/lib/projection/projectionUtils.ts"; + +type TCoordinateArray = Float32Array | Float64Array; + +export type TTriangularGrid = { + vertices: Float32Array; + centroidLatitudes: Float64Array; + centroidLongitudes: Float64Array; +}; + +function getVertex( + index: number, + vertexX: TCoordinateArray, + vertexY: TCoordinateArray, + vertexZ: TCoordinateArray +): [number, number, number] { + if (index < 0 || index >= vertexX.length) { + throw new Error(`Triangular grid vertex index ${index + 1} is invalid.`); + } + return [vertexX[index], vertexY[index], vertexZ[index]]; +} + +function writeVertex( + vertices: Float32Array, + triangleIndex: number, + vertexIndex: number, + vertex: [number, number, number] +) { + const offset = triangleIndex * 9 + vertexIndex * 3; + vertices[offset] = vertex[0]; + vertices[offset + 1] = vertex[1]; + vertices[offset + 2] = vertex[2]; +} + +function writeTriangle( + vertices: Float32Array, + triangleIndex: number, + vertex0: [number, number, number], + vertex1: [number, number, number], + vertex2: [number, number, number] +) { + if ( + shouldFlipCartesianTriangle( + vertex0[0], + vertex0[1], + vertex0[2], + vertex1[0], + vertex1[1], + vertex1[2], + vertex2[0], + vertex2[1], + vertex2[2] + ) + ) { + [vertex1, vertex2] = [vertex2, vertex1]; + } + writeVertex(vertices, triangleIndex, 0, vertex0); + writeVertex(vertices, triangleIndex, 1, vertex1); + writeVertex(vertices, triangleIndex, 2, vertex2); +} + +function buildTriangleCentroids(vertices: Float32Array) { + const triangleCount = vertices.length / 9; + const centroidLatitudes = new Float64Array(triangleCount); + const centroidLongitudes = new Float64Array(triangleCount); + for (let triangleIndex = 0; triangleIndex < triangleCount; triangleIndex++) { + const offset = triangleIndex * 9; + const centerX = + vertices[offset] + vertices[offset + 3] + vertices[offset + 6]; + const centerY = + vertices[offset + 1] + vertices[offset + 4] + vertices[offset + 7]; + const centerZ = + vertices[offset + 2] + vertices[offset + 5] + vertices[offset + 8]; + const { lat, lon } = ProjectionHelper.cartesianToLatLon( + centerX, + centerY, + centerZ + ); + centroidLatitudes[triangleIndex] = lat; + centroidLongitudes[triangleIndex] = lon; + } + return { centroidLatitudes, centroidLongitudes }; +} + +export function buildTriangularGrid( + vertexOfCell: Int32Array, + vertexX: TCoordinateArray, + vertexY: TCoordinateArray, + vertexZ: TCoordinateArray +): TTriangularGrid { + if (vertexOfCell.length % 3 !== 0) { + throw new Error("Triangular grid connectivity must contain three rows."); + } + if (vertexX.length !== vertexY.length || vertexX.length !== vertexZ.length) { + throw new Error("Triangular grid vertex coordinate lengths do not match."); + } + const triangleCount = vertexOfCell.length / 3; + const vertices = new Float32Array(triangleCount * 9); + for (let triangleIndex = 0; triangleIndex < triangleCount; triangleIndex++) { + const vertex0 = getVertex( + vertexOfCell[triangleIndex] - 1, + vertexX, + vertexY, + vertexZ + ); + const vertex1 = getVertex( + vertexOfCell[triangleCount + triangleIndex] - 1, + vertexX, + vertexY, + vertexZ + ); + const vertex2 = getVertex( + vertexOfCell[triangleCount * 2 + triangleIndex] - 1, + vertexX, + vertexY, + vertexZ + ); + writeTriangle(vertices, triangleIndex, vertex0, vertex1, vertex2); + } + return { vertices, ...buildTriangleCentroids(vertices) }; +} + +export function getTriangularBatchCount( + triangleCount: number, + batchSize: number +) { + return Math.ceil(triangleCount / batchSize); +} + +export function buildTriangularGeometryBatch( + grid: TTriangularGrid, + batchIndex: number, + batchSize: number, + projection: ProjectionHelper +): TGridPositionBatch { + const triangleCount = grid.vertices.length / 9; + const start = batchIndex * batchSize; + const end = Math.min(start + batchSize, triangleCount); + const sourceVertices = grid.vertices.subarray(start * 9, end * 9); + const positionValues = new Float32Array(sourceVertices.length); + const latLonValues = new Float32Array((sourceVertices.length / 3) * 2); + for ( + let vertexIndex = 0; + vertexIndex < sourceVertices.length / 3; + vertexIndex++ + ) { + const sourceOffset = vertexIndex * 3; + const { lat, lon } = ProjectionHelper.cartesianToLatLon( + sourceVertices[sourceOffset], + sourceVertices[sourceOffset + 1], + sourceVertices[sourceOffset + 2] + ); + projection.projectLatLonToArrays( + lat, + lon, + positionValues, + sourceOffset, + latLonValues, + vertexIndex * 2 + ); + } + return { batchIndex, positionValues, latLonValues }; +} + +export function buildTriangularDataBatch( + data: Float32Array, + batchIndex: number, + batchSize: number +): TGridDataValueBatch { + const start = batchIndex * batchSize; + const end = Math.min(start + batchSize, data.length); + const dataValues = new Float32Array((end - start) * 3); + for (let cellIndex = start; cellIndex < end; cellIndex++) { + const offset = (cellIndex - start) * 3; + dataValues[offset] = data[cellIndex]; + dataValues[offset + 1] = data[cellIndex]; + dataValues[offset + 2] = data[cellIndex]; + } + return { batchIndex, dataValues }; +} + +export function buildTriangularHoverIndexData( + grid: TTriangularGrid, + data: Float32Array +): TSerializedGeoSampleIndexData { + if (data.length !== grid.centroidLatitudes.length) { + throw new Error( + `Triangular grid has ${grid.centroidLatitudes.length} cells but data has ${data.length} values.` + ); + } + return buildSerializedGeoSampleIndexData( + grid.centroidLatitudes.slice(), + grid.centroidLongitudes.slice(), + data.slice() + ); +} diff --git a/src/lib/grids/triangularWorkerClient.ts b/src/lib/grids/triangularWorkerClient.ts new file mode 100644 index 0000000..833404d --- /dev/null +++ b/src/lib/grids/triangularWorkerClient.ts @@ -0,0 +1,83 @@ +import { + copyGridWorkerArray, + createGridGeometryWorkerClient, +} from "./gridGeometryWorkerClient.ts"; +import { GridGeometryWorkerMessageType } from "./gridGeometryWorkerProtocol.ts"; +import { + TriangularWorkerOperation, + type TTriangularGeometryWorkerRequest, + type TTriangularWorkerBatch, + type TTriangularWorkerMetadata, + type TTriangularWorkerRequest, +} from "./triangularWorkerProtocol.ts"; + +type TTriangularGeometryRequest = Omit< + TTriangularGeometryWorkerRequest, + "requestId" | "type" | "operation" +>; + +const client = createGridGeometryWorkerClient< + TTriangularWorkerRequest, + TTriangularWorkerMetadata, + TTriangularWorkerBatch +>( + () => + new Worker(new URL("./triangular.worker.ts", import.meta.url), { + type: "module", + }), + true +); + +export function buildTriangularGeometry( + request: TTriangularGeometryRequest, + callbacks: { + onMetadata: (metadata: TTriangularWorkerMetadata) => void; + onBatch: (batch: TTriangularWorkerBatch) => void; + } +) { + return client.build( + (requestId) => ({ + message: { + ...request, + requestId, + type: GridGeometryWorkerMessageType.BUILD, + operation: TriangularWorkerOperation.GEOMETRY, + projectionCenter: { + lat: request.projectionCenter.lat, + lon: request.projectionCenter.lon, + }, + }, + transfer: [ + request.vertexOfCell.buffer, + request.vertexX.buffer, + request.vertexY.buffer, + request.vertexZ.buffer, + ], + }), + callbacks + ); +} + +export function buildTriangularData( + request: { data: Float32Array; batchSize: number }, + callbacks: { + onMetadata: (metadata: TTriangularWorkerMetadata) => void; + onBatch: (batch: TTriangularWorkerBatch) => void; + } +) { + return client.build((requestId) => { + const data = copyGridWorkerArray(request.data); + return { + message: { + requestId, + type: GridGeometryWorkerMessageType.BUILD, + operation: TriangularWorkerOperation.DATA, + data, + batchSize: request.batchSize, + }, + transfer: [data.buffer], + }; + }, callbacks); +} + +export const terminateTriangularWorker = client.terminate; diff --git a/src/lib/grids/triangularWorkerProtocol.ts b/src/lib/grids/triangularWorkerProtocol.ts new file mode 100644 index 0000000..6a55071 --- /dev/null +++ b/src/lib/grids/triangularWorkerProtocol.ts @@ -0,0 +1,55 @@ +import type { TGridGeometryWorkerResponse } from "./gridGeometryWorkerProtocol.ts"; +import type { + TGridDataValueBatch, + TGridPositionBatch, +} from "./gridWorkerTypes.ts"; + +import type { + TProjectionCenter, + TProjectionType, +} from "@/lib/projection/projectionUtils.ts"; + +export const TriangularWorkerOperation = { + GEOMETRY: "geometry", + DATA: "data", +} as const; + +export type TTriangularWorkerOperation = + (typeof TriangularWorkerOperation)[keyof typeof TriangularWorkerOperation]; + +export type TTriangularGeometryWorkerRequest = { + requestId: number; + type: "build"; + operation: typeof TriangularWorkerOperation.GEOMETRY; + vertexOfCell: Int32Array; + vertexX: Float32Array | Float64Array; + vertexY: Float32Array | Float64Array; + vertexZ: Float32Array | Float64Array; + batchSize: number; + projectionType: TProjectionType; + projectionCenter: TProjectionCenter; +}; + +export type TTriangularDataWorkerRequest = { + requestId: number; + type: "build"; + operation: typeof TriangularWorkerOperation.DATA; + data: Float32Array; + batchSize: number; +}; + +export type TTriangularWorkerRequest = + | TTriangularGeometryWorkerRequest + | TTriangularDataWorkerRequest; + +export type TTriangularWorkerBatch = TGridPositionBatch | TGridDataValueBatch; + +export type TTriangularWorkerMetadata = { + operation: TTriangularWorkerOperation; + totalBatches: number; +}; + +export type TTriangularWorkerResponse = TGridGeometryWorkerResponse< + TTriangularWorkerMetadata, + TTriangularWorkerBatch +>; diff --git a/src/ui/grids/Curvilinear.vue b/src/ui/grids/Curvilinear.vue old mode 100755 new mode 100644 index 848fbed..e4e44fd --- a/src/ui/grids/Curvilinear.vue +++ b/src/ui/grids/Curvilinear.vue @@ -1,13 +1,10 @@ diff --git a/src/ui/grids/GaussianReduced.vue b/src/ui/grids/GaussianReduced.vue index e41e714..79563d7 100755 --- a/src/ui/grids/GaussianReduced.vue +++ b/src/ui/grids/GaussianReduced.vue @@ -1,13 +1,10 @@