diff --git a/src/main/scala/grib/GribParser.scala b/src/main/scala/grib/GribParser.scala index 158c9aa..da1f7fe 100644 --- a/src/main/scala/grib/GribParser.scala +++ b/src/main/scala/grib/GribParser.scala @@ -45,6 +45,8 @@ object GribParser { (conversion, bitsPerDataPoint, len5) = res5 ptr6 = ptr5 + len5 len6 <- parse6(path, ptr6) + ptr7 = ptr6 + len6 + len7 <- parse7(path, ptr7) sections = List( GribSection(0, ptr, len0), GribSection(1, ptr1, len1), @@ -52,6 +54,7 @@ object GribParser { GribSection(4, ptr4, len4), GribSection(5, ptr5, len5), GribSection(6, ptr6, len6), + GribSection(7, ptr7, len7), ) title = Codes.codesToString(meteo.discipline, meteo.category, meteo.product) grib = Grib(version, gribLength, title, grid, meteo, time, conversion, bitsPerDataPoint, sections) @@ -147,6 +150,13 @@ object GribParser { } yield length } + private def parse7(path: Path, ptr: Long): IO[Int] = { + for { + bytes <- readBytes(path, ptr, 64) + length = ByteBuffer.wrap(bytes.slice(0, 4)).getInt + } yield length + } + private def readBytes(path: Path, ptr: Long, len: Int): IO[Array[Byte]] = { Files[IO].readRange(path, 1024, ptr, ptr + len) .compile diff --git a/src/main/scala/server/Server.scala b/src/main/scala/server/Server.scala index f1bb7c2..1dbef52 100644 --- a/src/main/scala/server/Server.scala +++ b/src/main/scala/server/Server.scala @@ -7,7 +7,7 @@ import data.DataService import db.PostgresService import fetch.FetchService import parse.Aggregate -import server.ValidateRoutes.{AggFieldList, AggKey, CityList, DateTimeRange, Granularity, ValidateDate, ValidateDateTime, ValidateMonths, ValidateZonedDateTime} +import server.ValidateRoutes.{AggFieldList, AggKey, CityList, DateTimeRange, Granularity, ValidateDate, ValidateDateTime, ValidateMonths, ValidateInt, ValidateZonedDateTime} import io.circe.{Json, Printer} import org.http4s._ import org.http4s.dsl.io._ @@ -62,6 +62,9 @@ class Server(postgresService: PostgresService, fetch: FetchService, log: Logger[ case GET -> Root / "show" / "grib-list" => DataService.getFileList().flatMap(fileList => Ok(fileList.asJson)) + // TODO implement binary-chunk/ get request + case GET -> Root / "grib" / "binary-chunk" / ValidateInt(binaryOffset) / ValidateInt(binaryLength) / fileName => + DataService.getBinaryChunk(binaryOffset, binaryLength, fileName).flatMap(buffer => Ok(buffer)) diff --git a/src/main/scala/server/ValidateRoutes.scala b/src/main/scala/server/ValidateRoutes.scala index e700703..98c128e 100644 --- a/src/main/scala/server/ValidateRoutes.scala +++ b/src/main/scala/server/ValidateRoutes.scala @@ -89,4 +89,10 @@ object ValidateRoutes { } } } + + object ValidateInt { + def unapply(str: String): Option[Int] = { + Option(str.toInt) + } + } } diff --git a/web/src/harmonie/GribFile.tsx b/web/src/harmonie/GribFile.tsx index 2ae134d..5598fa3 100644 --- a/web/src/harmonie/GribFile.tsx +++ b/web/src/harmonie/GribFile.tsx @@ -2,13 +2,20 @@ import { Accessor, Component, createSignal } from 'solid-js' import styles from './harmonie.module.css' import { apiHost } from '../consts' +import { GribMessage } from './interfaces' +import { fetchBuffer } from '../helpers/fetch' +import { drawGrib } from './draw/drawGrib' + +const CROP_BOUNDS = { x: 1906-1-400, y: 950, width: 400, height: 300 } export const GribFile: Component<{ name: string, isActive: Accessor, + getCanvas: Accessor, onClick: (name: string) => void, -}> = ({ name, isActive, onClick }) => { - const [getStructure, setStructure] = createSignal([]) +}> = ({ name, isActive, getCanvas, onClick }) => { + const [getStructure, setStructure] = createSignal([]) + const isCrop = createSignal(false) function onFileClick() { onClick(name) @@ -20,14 +27,55 @@ export const GribFile: Component<{ } } + let cachedMessages: GribMessage[] = [] + let cachedBuffers: Uint8Array[] = [] + let cachedBitmasks: Uint8Array[] = [] + + function onParamClick(paramId: number) { + const grib = getStructure()[paramId] + console.log("onParamClick", grib) + const bitmaskSection = grib.sections.find(section => section.id === 6) + const binarySection = grib.sections.find(section => section.id === 7) + if (!bitmaskSection || !binarySection) return; + console.log("sections 6,7") + + const bitmaskOffset = bitmaskSection.offset + 6 + const bitmaskLength = bitmaskSection.size - 6 + const bitmaskPromise = bitmaskSection.size > 6 + ? fetchBuffer(`${apiHost}/api/grib/binary-chunk/${bitmaskOffset}/${bitmaskLength}/${name}`).then(b=>[b]) + : Promise.resolve([]) + + const binaryOffset = binarySection.offset + 5 + const binaryLength = binarySection.size - 5 + const fetchPromise: Promise<[GribMessage[], ArrayBuffer[], ArrayBuffer[]]> = + Promise.all([ + [grib], + fetchBuffer(`${apiHost}/api/grib/binary-chunk/${binaryOffset}/${binaryLength}/${name}`).then(b=>[b]), + bitmaskPromise, + ]) + + fetchPromise.then(([messages, binaryBuffers, bitmasks]) => { + cachedMessages = messages + cachedBuffers = binaryBuffers.map(b => new Uint8Array(b)) + cachedBitmasks = bitmasks.map(b => new Uint8Array(b)) + const colors: [string, string] = ['#0000ff', '#ffff00'] + const cropBounds = isCrop[0]() ? CROP_BOUNDS : undefined + drawGrib(getCanvas()!, cachedMessages, cachedBuffers, cachedBitmasks, colors, cropBounds) + }) + .catch(err => console.warn(err.message)) + // .finally(() => setIsLoading(false)) + } + return
  • { trimName(name) }
      - { getStructure().map(grib => -
    • { grib.title }
      • + { getStructure() + .sort((a, b) => a.title > b.title ? 1 : -1) + .map((grib, i) => +
    • onParamClick(i)}>{ grib.title }
      • )}
    • diff --git a/web/src/harmonie/Harmonie.tsx b/web/src/harmonie/Harmonie.tsx index 5500fe0..b0a528a 100644 --- a/web/src/harmonie/Harmonie.tsx +++ b/web/src/harmonie/Harmonie.tsx @@ -7,6 +7,7 @@ import { GribFile } from './GribFile' export const Harmonie: Component<{}> = () => { const [getFileList, setFileList] = createSignal([]) const [getActiveGrib, setActiveGrib] = createSignal('') + const [getCanvas, setCanvas] = createSignal() fetch(`${apiHost}/api/show/grib-list`) .then(re => re.json()) @@ -22,11 +23,14 @@ export const Harmonie: Component<{}> = () => { getActiveGrib() === fileName} + getCanvas={getCanvas} onClick={() => setActiveGrib(fileName)} /> )} -
    here goes pictures
    +
    + +
    } diff --git a/web/src/harmonie/draw/bitmask.ts b/web/src/harmonie/draw/bitmask.ts new file mode 100644 index 0000000..72629ce --- /dev/null +++ b/web/src/harmonie/draw/bitmask.ts @@ -0,0 +1,34 @@ +import { u8ToBits } from '../../helpers/u8ToBits.js' +import { MeteoGrid } from '../interfaces.js' + +export function applyBitmask( + grid: MeteoGrid, + buffer: Uint8Array, + bitmask: Uint8Array, + bytesPerPoint: number, +): Uint8Array { + const newBuffer = new Uint8Array(grid.rows * grid.cols * bytesPerPoint) + + let i=0, bufferI=0 + for (; i= newBuffer.length) { + break; + } + if (bits[bitI]) { + newBuffer[newI] = buffer[bufferI] + newBuffer[newI+1] = buffer[bufferI+1] + newBuffer[newI+2] = buffer[bufferI+2] + bufferI += bytesPerPoint + } else { + newBuffer[newI] = 255 + newBuffer[newI+1] = 255 + newBuffer[newI+2] = 255 + } + } + } + + return newBuffer +} diff --git a/web/src/harmonie/draw/bounds.ts b/web/src/harmonie/draw/bounds.ts new file mode 100644 index 0000000..f8ce4aa --- /dev/null +++ b/web/src/harmonie/draw/bounds.ts @@ -0,0 +1,32 @@ +import { GribMessage } from '../../interfaces/interfaces' +import { CropBounds } from './drawGrib' + +export function extractFromBounds( + grib: GribMessage, + source: Uint8Array, + cropBounds: CropBounds, +): Uint8Array { + const { grid, bitsPerDataPoint } = grib + const bytesPerPoint = bitsPerDataPoint / 8 + const { x, y, width, height } = cropBounds + + if ( + x < 0 + || y < 0 + || x + width > grid.cols-1 + || y + height > grid.rows-1 + ) { + throw new Error('Extract bbox out of grid bounds') + } + + const output = new Uint8Array(width*height*bytesPerPoint) + for (let row=y, i=0; row < y+height; row++) { + const inputOffset = (row*grid.cols + x)*bytesPerPoint + const readBytes = width*bytesPerPoint + const inputBuffer = source.slice(inputOffset, inputOffset+readBytes) + output.set(inputBuffer, i) + i += readBytes + } + + return output +} diff --git a/web/src/harmonie/draw/categoricalRain.ts b/web/src/harmonie/draw/categoricalRain.ts new file mode 100644 index 0000000..3770ab5 --- /dev/null +++ b/web/src/harmonie/draw/categoricalRain.ts @@ -0,0 +1,43 @@ +type RGBAu8 = [number, number, number, number] + +const FOG: RGBAu8 = [254, 248, 0, 255] // 'fef800' +const FOG_FREEZING: RGBAu8 = [194, 188, 0, 255] // 'c2bc00' +const RAIN_LIGHT: RGBAu8 = [0, 255, 0, 255] // '00ff00' +const RAIN_MODERATE: RGBAu8 = [0, 197, 0, 255] // '00c500' +const RAIN_HEAVY: RGBAu8 = [0, 132, 0, 255] // '008400' +const ICE_LIGHT: RGBAu8 = [255, 32, 55, 255] // 'ff2037' +const ICE_HEAVY: RGBAu8 = [220, 0, 0, 255] // 'dc0000' +export function categoricalRainColors(value: number): RGBAu8 { + switch (value) { + case 0: + return FOG + case 1*32: + return FOG_FREEZING + case 2*32: + return RAIN_LIGHT + case 3*32: + return RAIN_MODERATE + case 4*32: + return RAIN_HEAVY + case 5*32: + return ICE_LIGHT + case 6*32: + return ICE_HEAVY + default: + return [255, 255, 255, 0] + } +} + +/* FROM CHATGPT +0 = No rain +1 * 32 = 32 (Drizzle) +2 * 32 = 64 (Light rain) +3 * 32 = 96 (Moderate rain) +4 * 32 = 128 (Heavy rain) +5 * 32 = 160 (Very heavy rain, if applicable) +6 * 32 = 192 (Extreme rain, if applicable) +*/ + +export function hexToU8(hex: string): [number, number, number] { + return [parseInt('0x'+hex.slice(0, 2)), parseInt('0x'+hex.slice(2, 4)), parseInt('0x'+hex.slice(4, 6))] +} diff --git a/web/src/harmonie/draw/constants.ts b/web/src/harmonie/draw/constants.ts new file mode 100644 index 0000000..f196b85 --- /dev/null +++ b/web/src/harmonie/draw/constants.ts @@ -0,0 +1,53 @@ +export type ColorEntry = { value: number, color: [number, number, number, number]} + +export const WIND_SPEED: ColorEntry[] = [ + { value: 40, color: [170, 0, 190, 255] }, + { value: 35, color: [170, 0, 190, 255] }, + { value: 30, color: [225, 20, 0, 255] }, + { value: 25, color: [225, 20, 0, 255] }, + { value: 20, color: [255, 160, 0, 255] }, + { value: 15, color: [255, 250, 170, 255] }, + { value: 10, color: [80, 240, 80, 255] }, + { value: 5, color: [45, 155, 150, 255] }, + { value: 0, color: [180, 240, 250, 255] }, +] + +export const PRECIPITATION: ColorEntry[] = [ + { value: 30, color: [126, 26, 99, 255] }, + { value: 25, color: [120, 27, 131, 255] }, + { value: 20, color: [84, 20, 130, 255] }, + { value: 15, color: [49, 16, 129, 255] }, + { value: 10, color: [9, 15, 129, 255] }, + { value: 6, color: [0, 24, 150, 255] }, + { value: 4, color: [0, 43, 186, 255] }, + { value: 2, color: [10, 70, 220, 255] }, + { value: 1, color: [40, 109, 246, 255] }, + { value: 0.5, color: [80, 162, 248, 255] }, + { value: 0.2, color: [118, 202, 249, 255] }, + { value: 0.1, color: [152, 233, 252, 255] }, + { value: 0.05, color: [255, 255, 255, 255] }, +] + +// from ltv +export const TEMPERATURES: ColorEntry[] = [ + { value: 35, color: [155, 30, 30, 255] }, + { value: 10, color: [250, 225, 5, 255] }, + { value: 0, color: [80, 190, 240, 255] }, + { value: -15, color: [30, 70, 155, 255] }, + { value: -30, color: [140, 30, 190, 255] }, +] + +// from yr.no +// export const TEMPERATURES: ColorEntry[] = [ +// { value: 50, color: [133, 0, 62, 255] }, +// { value: 40, color: [195, 0, 0, 255] }, +// { value: 30, color: [255, 76, 56, 255] }, +// { value: 20, color: [255, 175, 111, 255] }, +// { value: 10, color: [255, 243, 81, 255] }, +// { value: 0, color: [195, 246, 215, 255] }, +// { value: -10, color: [94, 231, 240, 255] }, +// { value: -20, color: [63, 201, 243, 255] }, +// { value: -30, color: [79, 157, 232, 255] }, +// { value: -40, color: [0, 81, 163, 255] }, +// { value: -50, color: [79, 15, 134, 255] }, +// ] diff --git a/web/src/harmonie/draw/drawGrib.ts b/web/src/harmonie/draw/drawGrib.ts new file mode 100644 index 0000000..31410a6 --- /dev/null +++ b/web/src/harmonie/draw/drawGrib.ts @@ -0,0 +1,163 @@ +import { interpolateColors } from '../../helpers/interpolateColors.ts' +import { GribMessage, MeteoParam } from '../../interfaces/interfaces.ts' +import { applyBitmask } from './bitmask.ts' +import { extractFromBounds } from './bounds.ts' +import { categoricalRainColors } from './categoricalRain.ts' +import { precipitationColors } from './precipitation.ts' +import { temperatureColors } from './temperature.ts' +import { windDirectionArrows, windDirectionColors, windSpeedColors } from './windDirection.ts' + +export type CropBounds = { x: number, y: number, width: number, height: number } + +export function drawGrib( + canvas: HTMLCanvasElement, + messages: GribMessage[], + buffers: Uint8Array[], + bitmasks: Uint8Array[], + colors: [string, string], + cropBounds: CropBounds | undefined, +): void { + // normally we have one message/buffer/bitmask?. special cases have multiple like wind direction + const [grib] = messages + + let { grid } = grib + let { cols, rows } = grid + + let modifiedBuffers = buffers.map((buffer, i) => { + const bytesPerPoint = messages[i].bitsPerDataPoint / 8 + return bitmasks[i] ? applyBitmask(grid, buffer, bitmasks[i], bytesPerPoint) : buffer + }) + + if (cropBounds) { + modifiedBuffers = modifiedBuffers.map(buffer => extractFromBounds(grib, buffer, cropBounds)) + cols = cropBounds.width + rows = cropBounds.height + } + + canvas.width = cols + canvas.height = rows + canvas.style.width = '100%' + canvas.style.minWidth = '1280px' + const ctx = canvas.getContext('2d')! + let imgData = ctx.createImageData(cols, rows) + + fillImageData(imgData, messages, modifiedBuffers, colors) + if (grib.meteo.discipline === 0 && grib.meteo.category === 2 && grib.meteo.product === 192) { + imgData = windDirectionArrows(imgData, messages, modifiedBuffers) + } + + const tempCanvas = document.createElement('canvas') + const tempCtx = tempCanvas.getContext('2d')! + tempCanvas.width = imgData.width + tempCanvas.height = imgData.height + tempCtx.putImageData(imgData, 0, 0) + + ctx.save() + ctx.scale(1, -1) + ctx.drawImage(tempCanvas, 0, -canvas.height) + // ctx.drawImage(tempCanvas, 0, 0) + ctx.restore() +} + +const CATEGORICAL_RAIN = [0, 1, 192] +const TOTAL_PRECIPITATION = [0, 1, 52] +const RAIN_PRECIPITATION = [0, 1, 65] +const TEMPERATURE = [0, 0, 0] +const WIND_DIRECTION = [0, 2, 192] +const WIND_SPEED = [0, 2, 1] +const WIND_SPEED_GUST = [0, 2, 22] + +function fillImageData( + imgData: ImageData, + messages: GribMessage[], + buffers: Uint8Array[], + colors: [string, string], +) { + const [grib] = messages + const [buffer] = buffers + + const { meteo, conversion, bitsPerDataPoint } = grib + const bytesPerPoint = grib.bitsPerDataPoint / 8 + const fromColor = rgbHexToU8(colors[0]) + const toColor = rgbHexToU8(colors[1]) + + const cols = imgData.width + const rows = imgData.height + + for (let row = 0; row < rows; row++) { + for (let col = 0; col < cols; col++) { + + const bufferI = (row * cols + col) * bytesPerPoint + const index = (row * cols + col) * 4 + const firstByte = buffer[bufferI] + const encodedValue = toInt(buffer.slice(bufferI, bufferI+bitsPerDataPoint/8)) + + let color = [255, 255, 255, 255] + if (isMeteoEqual(meteo, CATEGORICAL_RAIN)) { + color = categoricalRainColors(firstByte) + } + else if (isMeteoEqual(meteo, TOTAL_PRECIPITATION)) { + color = precipitationColors(encodedValue, conversion) + } + else if (isMeteoEqual(meteo, RAIN_PRECIPITATION)) { + color = precipitationColors(encodedValue, conversion) + } + else if (isMeteoEqual(meteo, TEMPERATURE)) { + color = temperatureColors(encodedValue, conversion) + } + else if (isMeteoEqual(meteo, WIND_DIRECTION)) { + const [, bufferU, bufferV] = buffers + const encodedValU = toInt(bufferU.slice(bufferI, bufferI+bitsPerDataPoint/8)) + const encodedValV = toInt(bufferV.slice(bufferI, bufferI+bitsPerDataPoint/8)) + const [, metaU, metaV] = messages // first message fake one 0-2-192 + color = windDirectionColors(encodedValU, encodedValV, metaU!.conversion, metaV!.conversion) + } + else if (isMeteoEqual(meteo, WIND_SPEED)) { + color = windSpeedColors(encodedValue, conversion) + } + else if (isMeteoEqual(meteo, WIND_SPEED_GUST)) { + color = windSpeedColors(encodedValue, conversion) + } + else { + color = interpolateColors(firstByte, fromColor, toColor) + } + + imgData.data[index] = color[0] + imgData.data[index + 1] = color[1] + imgData.data[index + 2] = color[2] + imgData.data[index + 3] = color[3] + } + } +} + + +function rgbHexToU8(hex: string): RGBu8 { + return [ + parseInt(`0x${hex.slice(1, 3)}`), + parseInt(`0x${hex.slice(3, 5)}`), + parseInt(`0x${hex.slice(5, 7)}`), + ] +} + +type RGBu8 = [number, number, number] + +export function isMeteoEqual(meteo: MeteoParam, arr: number[]): boolean { + const arr2 = [meteo.discipline, meteo.category, meteo.product] + return arr.length === arr2.length && arr.every((value, index) => value === arr2[index]) +} + +function toInt(bytes: Uint8Array): number { + return bytes.reduce((acc, curr) => acc * 256 + curr) +} + +export function toSignedInt(bytes: Uint8Array): number { + const unsigned = toInt(bytes) + + const signBit = 1 << (bytes.length * 8 - 1) // Example: 16-bit -> 0x8000 + if (unsigned & signBit) { + // If the sign bit is set, compute the two's complement + return unsigned - (1 << (bytes.length * 8)) + } + + return unsigned // If the sign bit is not set, return as is +} diff --git a/web/src/harmonie/draw/precipitation.ts b/web/src/harmonie/draw/precipitation.ts new file mode 100644 index 0000000..093e4d5 --- /dev/null +++ b/web/src/harmonie/draw/precipitation.ts @@ -0,0 +1,17 @@ +import { valueToColorInterpolated, valueToColorThreshold } from '../../helpers/interpolateColors' +import { MeteoConversion } from '../../interfaces/interfaces' +import { PRECIPITATION } from './constants' + + +export function precipitationColors( + encodedValue: number, + { reference, binaryScale, decimalScale}: MeteoConversion, + isInterpolated = true, +): [number, number, number, number] { + // const rainMM = (reference + encodedValue * Math.pow(2, binaryScale)) * Math.pow(10, -decimalScale) / 1000 + const rainMM = (reference + encodedValue * Math.pow(2, binaryScale)) * Math.pow(10, -decimalScale) + + return isInterpolated + ? valueToColorInterpolated(rainMM, PRECIPITATION) + : valueToColorThreshold(rainMM, PRECIPITATION) +} diff --git a/web/src/harmonie/draw/temperature.ts b/web/src/harmonie/draw/temperature.ts new file mode 100644 index 0000000..839f730 --- /dev/null +++ b/web/src/harmonie/draw/temperature.ts @@ -0,0 +1,16 @@ +import { valueToColorInterpolated, valueToColorThreshold } from '../../helpers/interpolateColors' +import { MeteoConversion } from '../../interfaces/interfaces' +import { TEMPERATURES } from './constants' + + +export function temperatureColors( + encodedValue: number, + { reference, binaryScale, decimalScale}: MeteoConversion, + isInterpolated = true, +): [number, number, number, number] { + const temperatureC = (reference + encodedValue * Math.pow(2, binaryScale)) * Math.pow(10, -decimalScale) - 273.15 + + return isInterpolated + ? valueToColorInterpolated(temperatureC, TEMPERATURES) + : valueToColorThreshold(temperatureC, TEMPERATURES) +} diff --git a/web/src/harmonie/draw/windDirection.ts b/web/src/harmonie/draw/windDirection.ts new file mode 100644 index 0000000..3c06ff3 --- /dev/null +++ b/web/src/harmonie/draw/windDirection.ts @@ -0,0 +1,174 @@ +import { apiHost } from '../../consts' +import { fetchBuffer } from '../../helpers/fetch' +import { valueToColorInterpolated } from '../../helpers/interpolateColors' +import { GribMessage, MeteoConversion } from '../interfaces' +import { WIND_SPEED } from './constants' +import { rotateWind } from './windRotate' + +const CELL_SIZE = 12 + +export function windDirectionArrows( + imgData: ImageData, + messages: GribMessage[], + buffers: Uint8Array[], +) { + const [, metaU, metaV] = messages + const { conversion: convU } = metaU + const { conversion: convV } = metaV + const [, bufferU, bufferV] = buffers + const cols = imgData.width + const rows = imgData.height + const canvas = document.createElement('canvas') + canvas.width = imgData.width + canvas.height = imgData.height + const ctx = canvas.getContext('2d')! + ctx.putImageData(imgData, 0, 0) // draw wind speed color below direction + const directions: number[][] = [] + const bytesPerPoint = metaU.bitsPerDataPoint/8 + + const lambert = (messages[0].grid as any).lambert + let rot_lat = lambert[0] / 1_000_000 + let rot_lon = lambert[1] / 1_000_000 + // rot_lat = -9999.0 + // rot_lon = -9999.0 + let reg_lat = 56.530592 + let reg_lon = -2.918742 + + for (let row = 0; row < rows; row++) { + directions[row] = [] + for (let col = 0; col < cols; col++) { + const bufferI = (row * cols + col) * bytesPerPoint + + const encodedU = toInt(bufferU.slice(bufferI, bufferI+bytesPerPoint)) + const encodedV = toInt(bufferV.slice(bufferI, bufferI+bytesPerPoint)) + const windSpeedU = (convU.reference + encodedU * Math.pow(2, convU.binaryScale)) * Math.pow(10, -convU.decimalScale) + const windSpeedV = (convV.reference + encodedV * Math.pow(2, convV.binaryScale)) * Math.pow(10, -convV.decimalScale) + + // const directionDeg = (Math.atan2(windSpeedU, windSpeedV)*180/Math.PI + 360 +45) % 360 + const directionDeg = rotateWind(rot_lat, rot_lon, reg_lat, reg_lon, windSpeedU, windSpeedV)[0] + const directionRad = (Math.PI / 180) * directionDeg + + directions[row][col] = directionRad + } + } + + const gridH = Math.floor(directions.length/CELL_SIZE) + const gridW = Math.floor(directions[0].length/CELL_SIZE) + for (let row = 0; row < gridH; row++) { + for (let col = 0; col < gridW; col++) { + const directionAvg = true + ? getAvgDirection(directions, row, col) + : getDirection(directions, row, col) + + const centerX = col * CELL_SIZE + CELL_SIZE/2 + const centerY = row * CELL_SIZE + CELL_SIZE/2 + const arrowLength = CELL_SIZE*0.9 + + ctx.save() + ctx.translate(centerX, centerY) + ctx.rotate(directionAvg) + + ctx.beginPath(); + ctx.moveTo(-arrowLength / 2, 0) + ctx.lineTo(arrowLength / 2, 0) + ctx.stroke() + + const arrowheadSize = CELL_SIZE/4 + ctx.beginPath() + ctx.moveTo(arrowLength / 2, 0) + ctx.lineTo(arrowLength / 2 - arrowheadSize, -arrowheadSize / 2) + ctx.lineTo(arrowLength / 2 - arrowheadSize, arrowheadSize / 2) + ctx.closePath() + ctx.fill() + + ctx.restore() + } + } + + return ctx.getImageData(0, 0, imgData.width, imgData.height) +} + +function getDirection(directions: number[][], gridRow: number, gridCol: number) { + const directionRow = gridRow * CELL_SIZE + Math.round(CELL_SIZE/2) + const directionCol = gridCol * CELL_SIZE + Math.round(CELL_SIZE/2) + return directions[directionRow][directionCol] +} + +// in radians +function getAvgDirection(directions: number[][], gridRow: number, gridCol: number) { + let sumSin = 0; // Sum of sine components + let sumCos = 0; // Sum of cosine components + + for (let row = 0; row < CELL_SIZE; row++) { + for (let col = 0; col < CELL_SIZE; col++) { + const directionRow = gridRow * CELL_SIZE + row; + const directionCol = gridCol * CELL_SIZE + col; + + const directionRad = directions[directionRow][directionCol] + + sumSin += Math.sin(directionRad); + sumCos += Math.cos(directionRad); + } + } + + // Compute the average direction + const avgDirection = Math.atan2(sumSin, sumCos); // Result is in radians + return (avgDirection + 2 * Math.PI) % (2 * Math.PI); // Ensure the result is in [0, 2π) +} + +// actually calculates and draws wind speed +export function windDirectionColors( + encodedU: number, + encodedV: number, + convU: MeteoConversion, + convV: MeteoConversion, +) { + const windSpeedU = (convU.reference + encodedU * Math.pow(2, convU.binaryScale)) * Math.pow(10, -convU.decimalScale) + const windSpeedV = (convV.reference + encodedV * Math.pow(2, convV.binaryScale)) * Math.pow(10, -convV.decimalScale) + const windSpeed = Math.sqrt(Math.pow(windSpeedU, 2) + Math.pow(windSpeedV, 2)) + + return valueToColorInterpolated(windSpeed, WIND_SPEED) +} + +export function windSpeedColors( + encodedValue: number, + { reference, binaryScale, decimalScale}: MeteoConversion, +) { + const windSpeed = (reference + encodedValue * Math.pow(2, binaryScale)) * Math.pow(10, -decimalScale) + return valueToColorInterpolated(windSpeed, WIND_SPEED) +} + +export function fetchWindData( + customMessage: GribMessage, + gribArr: GribMessage[], +): Promise<[GribMessage[], ArrayBuffer[], ArrayBuffer[]]> { + const windU = gribArr.find(m => m.meteo.discipline===0 && m.meteo.category===2 && m.meteo.product===2 && m.meteo.levelType===103 && m.meteo.levelValue===10) + const windV = gribArr.find(m => m.meteo.discipline===0 && m.meteo.category===2 && m.meteo.product===3 && m.meteo.levelType===103 && m.meteo.levelValue===10) + if (!windU || !windV) throw new Error('Didnt found u/v components of wind') + + const section7u = windU.sections.find(section => section.id === 7) + const section7v = windV.sections.find(section => section.id === 7) + if (!section7u || !section7v) throw new Error('Didnt found binaru section for wind u/v') + + const uBinaryOffset = section7u.offset + 5 + const uBinaryLength = section7u.size - 5 + + const vBinaryOffset = section7v.offset + 5 + const vBinaryLength = section7v.size - 5 + + return Promise.all([ + fetchBuffer(`${apiHost}/api/binary-chunk/${uBinaryOffset}/${uBinaryLength}`), + fetchBuffer(`${apiHost}/api/binary-chunk/${vBinaryOffset}/${vBinaryLength}`), + ]).then(([bufferU, bufferV]) => { + const buffer = new Uint8Array(bufferU.byteLength + bufferV.byteLength) + buffer.set(new Uint8Array(bufferU)) + buffer.set(new Uint8Array(bufferV), bufferU.byteLength) + const messages = [customMessage, windU, windV] + const buffers = [bufferU, bufferU, bufferV] + return [messages, buffers, []] + }) +} + +function toInt(bytes: Uint8Array): number { + return bytes.reduce((acc, curr) => acc * 256 + curr) +} diff --git a/web/src/harmonie/draw/windRotate.ts b/web/src/harmonie/draw/windRotate.ts new file mode 100644 index 0000000..69704a0 --- /dev/null +++ b/web/src/harmonie/draw/windRotate.ts @@ -0,0 +1,135 @@ +// converted from C code https://opendatadocs.dmi.govcloud.dk/Data/Forecast_Data_Weather_Model_HARMONIE +export function rotateWind( + rot_lat: number, + rot_lon: number, + reg_lat: number, + reg_lon: number, + u_in: number, + v_in: number, + southpole_lat = 26.5, + southpole_lon = -40, +): [direction: number, strength: number, u_out: number, v_out: number] +/* Given either a point in the regular grid (set `*rot_lat' <= -999.0) +* or a point in the rotated grid, calculate the corresponding point +* in the opposite grid, change the (u, v)-vector from rotated to +* regular grid and calculate the wind force (`*strength') and the +* wind direction in the regular grid. `southpole_lat' and +* `southpole_lon' defines the coordinate of the southpole in +* the roated grid */ +{ + /* Find the missing point, whether is is the rotated or the regular */ + if (rot_lat <= -999.0) [rot_lat, rot_lon] = reg2rot(reg_lat, reg_lon, southpole_lat, southpole_lon) + else [reg_lat, reg_lon] = rot2reg(rot_lat, rot_lon, southpole_lat, southpole_lon); + + /* Calculate the wind strength */ + const strength = Math.sqrt(u_in*u_in + v_in*v_in); + + /* Add a small distance in the direction of the wind to the rotated + * grid point, changing the distance into degrees */ + const rot_lat2 = rot_lat + 0.1*v_in/(strength); + let clat = Math.cos(rot_lat*Math.PI/180.0) + if (0.0001 > clat && clat > -0.0001) { + throw new Error("Internal error: Too close to pole to calculate rotated wind") + } + const rot_lon2 = rot_lon + 0.1*u_in/(strength * clat); + + /* Translate new rotated grid point to regular grid */ + const [reg_lat2, reg_lon2] = rot2reg(rot_lat2, rot_lon2, southpole_lat, southpole_lon) + + /* Transform offset in lat-lon to offset in x-y */ + clat = Math.cos(reg_lat*Math.PI/180.0) + const dx = clat*(reg_lon2 - reg_lon) + + /* Calculate the direction of the wind vector in the regular grid */ + const direc = Math.atan2(reg_lat2 - reg_lat, dx); + + /* Regular direction in degrees */ + let direction = 630.0 - direc*180.0 / Math.PI + while (direction > 360.0) direction -= 360.0 + + const u_out = Math.cos(direc) * strength + const v_out = Math.sin(direc) * strength + + return [direction, strength, u_out, v_out] +} + + +function rot2reg( + rot_lat: number, + rot_lon: number, + southpole_lat: number, + southpole_lon: number, +): [reg_lat: number, reg_lon: number] +/* Convert from rotated latitude-longitude to regular latitude-longitude +with the transformation defined by the southpole coordinates. +Coordinates are given in degrees N (negative for S) and degrees E +(negative for W). */ +{ + const to_rad = Math.PI/180.0 + const to_deg = 1.0/to_rad + + const sin_y_cen = Math.sin(to_rad*(southpole_lat + 90.0)) + const cos_y_cen = Math.cos(to_rad*(southpole_lat + 90.0)) + + const sin_x_rot = Math.sin(to_rad*rot_lon) + const cos_x_rot = Math.cos(to_rad*rot_lon) + const sin_y_rot = Math.sin(to_rad*rot_lat) + const cos_y_rot = Math.cos(to_rad*rot_lat) + let sin_y_reg = cos_y_cen*sin_y_rot + sin_y_cen*cos_y_rot*cos_x_rot + if (sin_y_reg < -1.0) sin_y_reg = -1.0 + if (sin_y_reg > 1.0) sin_y_reg = 1.0 + + const reg_lat = to_deg*Math.asin(sin_y_reg) + + const cos_y_reg = Math.cos(reg_lat*to_rad); + let cos_lon_rad = (cos_y_cen*cos_y_rot*cos_x_rot - sin_y_cen*sin_y_rot)/cos_y_reg; + if (cos_lon_rad < -1.0) cos_lon_rad = -1.0; + if (cos_lon_rad > 1.0) cos_lon_rad = 1.0; + const sin_lon_rad = cos_y_rot*sin_x_rot/cos_y_reg; + let lon_rad = Math.acos(cos_lon_rad); + if (sin_lon_rad < 0.0) lon_rad = -lon_rad; + + const reg_lon = to_deg*lon_rad + southpole_lon; + + return [reg_lat, reg_lon] +} + + +function reg2rot( + reg_lat: number, + reg_lon: number, + southpole_lat: number, + southpole_lon: number, +): [rot_lat: number, rot_lon: number] +/* Convert from regular latitude-longitude to rotated latitude-longitude +with the transformation defined by the southpole coordinates. +Coordinates are given in degrees N (negative for S) and degrees E +(negative for W). */ +{ + const to_rad = Math.PI/180.0 + const to_deg = 1.0/to_rad + const sin_y_cen = Math.sin(to_rad*(southpole_lat + 90.0)); + const cos_y_cen = Math.cos(to_rad*(southpole_lat + 90.0)); + + const lon_rad = to_rad*(reg_lon - southpole_lon) + const sin_lon_rad = Math.sin(lon_rad) + const cos_lon_rad = Math.cos(lon_rad) + const sin_y_reg = Math.sin(to_rad*reg_lat) + const cos_y_reg = Math.cos(to_rad*reg_lat) + let sin_y_rot = cos_y_cen*sin_y_reg - sin_y_cen*cos_y_reg*cos_lon_rad + if (sin_y_rot < -1.0) sin_y_rot = -1.0 + if (sin_y_rot > 1.0) sin_y_rot = 1.0 + + const rot_lat = Math.asin(sin_y_rot)*to_deg + + const cos_y_rot = Math.cos(rot_lat*to_rad) + let cos_x_rot = (cos_y_cen*cos_y_reg*cos_lon_rad + sin_y_cen*sin_y_reg)/cos_y_rot + if (cos_x_rot < -1.0) cos_x_rot = -1.0 + if (cos_x_rot > 1.0) cos_x_rot = 1.0 + const sin_x_rot = cos_y_reg*sin_lon_rad/cos_y_rot + + let rot_lon = Math.acos(cos_x_rot)*to_deg + if (sin_x_rot < 0.0) rot_lon = -rot_lon + + return [rot_lat, rot_lon] +} diff --git a/web/src/harmonie/interfaces.ts b/web/src/harmonie/interfaces.ts new file mode 100644 index 0000000..382dd6f --- /dev/null +++ b/web/src/harmonie/interfaces.ts @@ -0,0 +1,39 @@ +export type MeteoParam = { + discipline: number, // 0=meteo, 1=hydro, 2=land surface, 3=space products + category: number, // 0=temperature, 1=moisture, 6=cloud, 19=atmospheric + product: number, + subType: string, // now or over time avg/sum + levelType: number, // 102 - entire atmosphere, 103 - above ground, above sea level + levelValue: number, // meters above ground/sea level +} + +export type MeteoConversion = { + reference: number, // float + binaryScale: number, // int + decimalScale: number, // int +} + +export type MeteoGrid = { + cols: number, + rows: number, + template: number, // 0 - regular lat/lon +} + +export type GribMessage = { + offset: number, + size: number, + version: number, + title: string, + meteo: MeteoParam, + grid: MeteoGrid, + bitsPerDataPoint: number, + subType: string, + conversion: MeteoConversion, + sections: GribSection[], +} + +export type GribSection = { + offset: number, + size: number, + id: number, +} \ No newline at end of file diff --git a/web/src/helpers/fetch.ts b/web/src/helpers/fetch.ts new file mode 100644 index 0000000..7cd0501 --- /dev/null +++ b/web/src/helpers/fetch.ts @@ -0,0 +1,49 @@ +export const fetchJson = (url: string, options = {}) => { + return fetch(url, options).then(async response => { + let responseCode = '' + if (!response.ok) { + responseCode = `error ${response.status}` + } + + let data + let responseMessage = '' + try { + data = await response.json() + if (data.error) responseMessage = data.error + } catch (err) { + responseMessage = 'Invalid JSON response' + } + + if (responseCode) { + if (responseMessage) throw new Error(`${responseCode}: ${responseMessage}`) + else throw new Error(responseCode) + } else if (responseMessage) { + throw new Error(responseMessage) + } + + return data + }) +} + + +export const fetchBuffer = (url: string, options = {}) => { + return fetch(url, options).then(async response => { + let responseCode = '' + if (!response.ok) { + responseCode = `error ${response.status}` + } + + let data + try { + data = await response.arrayBuffer() + } catch (err) { + throw new Error('Invalid ArrayBuffer response') + } + + if (responseCode) { + throw new Error(responseCode) + } + + return data + }) +} diff --git a/web/src/helpers/interpolateColors.ts b/web/src/helpers/interpolateColors.ts new file mode 100644 index 0000000..af2523e --- /dev/null +++ b/web/src/helpers/interpolateColors.ts @@ -0,0 +1,47 @@ +import { ColorEntry } from '../grib/draw/constants' + +type RGBu8 = [number, number, number] +type RGBAu8 = [number, number, number, number] + +export function interpolateColors(value: number, a: RGBu8, b: RGBu8): RGBAu8 { + const color = a.slice(0).map((from, i) => { + const to = b[i] + const delta = (to - from) * (value/255) + return from + delta + }) + + return [...color, 255] as RGBAu8 +} + +export function valueToColorThreshold(value: number, colorArr: ColorEntry[]): [number, number, number, number] { + const entries = colorArr.filter(t => t.value <= value) + + return entries.length === 0 ? colorArr[colorArr.length-1].color : entries[0].color +} + +export function valueToColorInterpolated(value: number, colorArr: ColorEntry[]): [number, number, number, number] { + const minIdx = colorArr.length-1 + if (value >= colorArr[0].value) return colorArr[0].color + if (value <= colorArr[minIdx].value) return colorArr[minIdx].color + + let closeMax: ColorEntry = colorArr[0] + let closeMin: ColorEntry = colorArr[minIdx] + for (let i=0; i= value && currentDeg < closeMax.value) { + closeMax = colorArr[i] + } + + if (currentDeg <= value && currentDeg > closeMin.value) { + closeMin = colorArr[i] + } + } + + const maxDeg = closeMax.color.slice(0, 3) as [number, number, number] + const minDeg = closeMin.color.slice(0, 3) as [number, number, number] + const scale = Math.abs(closeMax.value - closeMin.value) + const delta = 255*(value - closeMin.value)/scale + const rgba = interpolateColors(delta, minDeg, maxDeg) + + return rgba +} diff --git a/web/src/helpers/u8ToBits.ts b/web/src/helpers/u8ToBits.ts new file mode 100644 index 0000000..87ea190 --- /dev/null +++ b/web/src/helpers/u8ToBits.ts @@ -0,0 +1,19 @@ +export type U8bits = [boolean, boolean, boolean, boolean, boolean, boolean, boolean, boolean] + +export function u8ToBits(dec: number): U8bits { + const str = (dec >>> 0).toString(2) + const str8 = leftPad(str) + const arr: U8bits = [false, false, false, false, false, false, false, false] + for (let i=0; i<8; i++) { + arr[i] = str8[i] === '1' ? true : false + } + return arr +} + +function leftPad(str: string): string { + const padCount = 8 - str.length + for(let i=0; i