import type OlMap from 'ol/Map'; import type { Extent } from 'ol/extent'; import { getWidth } from 'ol/extent'; import ImageLayer from 'ol/layer/Image'; import ImageCanvas from 'ol/source/ImageCanvas'; import { transform, transformExtent } from 'ol/proj'; import { obtainViewer2D } from '@/gis/ol/map'; import type { ContourWorkerRequest, ContourWorkerResponse, FlexpartContourLayerOptions, FlexpartContourLayerResult, FlexpartGridFrame, SerializedContour, } from './flexpartContourTypes'; import { PYTHON_CONTOUR_COLORS } from './flexpartContourTypes'; const WEB_MERCATOR_MAX_LATITUDE = 85.0511287798066; interface ProjectedContour { color: string; path: Path2D; } interface PendingRequest { resolve: (value: SerializedContour[]) => void; reject: (reason?: unknown) => void; } interface PreparedContourGrid { longitudes: readonly number[]; values: Float64Array; } export function createLogLevels( colorMin: number, colorMax: number, bandCount = 10, ): number[] { if (!Number.isFinite(colorMin) || !Number.isFinite(colorMax) || colorMin <= 0 || colorMax <= colorMin || bandCount < 1) { throw new Error('Invalid logarithmic color range.'); } const minExponent = Math.log10(colorMin); const maxExponent = Math.log10(colorMax); return Array.from({ length: bandCount + 1 }, (_, index) => ( 10 ** ( minExponent + (maxExponent - minExponent) * index / bandCount ) )); } export class FlexpartContourRenderer { private readonly map: OlMap; private readonly source: ImageCanvas; private readonly layer: ImageLayer; private readonly worker: Worker; private readonly pending = new Map(); private requestId = 0; private frameVersion = 0; private projectedContours: ProjectedContour[] = []; private layerExtent?: Extent; private readonly colors: readonly string[]; private readonly wrapX: boolean; private readonly worldWidth?: number; constructor( map: OlMap, options: FlexpartContourLayerOptions = {}, ) { this.map = map; this.colors = options.colors ?? PYTHON_CONTOUR_COLORS; this.wrapX = options.wrapX ?? true; const projection = map.getView().getProjection(); const projectionCode = projection.getCode(); const projectionExtent = projection.getExtent(); this.worldWidth = this.wrapX && projection.canWrapX() && projectionExtent ? getWidth(projectionExtent) : undefined; this.source = new ImageCanvas({ projection: projectionCode, ratio: 1, canvasFunction: ( extent, _resolution, _pixelRatio, size, ) => this.drawCanvas(extent, size), }); this.layer = new ImageLayer({ source: this.source, opacity: options.opacity ?? 0.78, zIndex: options.zIndex ?? 1000, }); if (options.layerName) { this.layer.set('name', options.layerName); } this.worker = new Worker( new URL('./flexpartContour.worker.ts', import.meta.url), { type: 'module' }, ); this.worker.onmessage = ( event: MessageEvent, ) => { const request = this.pending.get(event.data.id); if (!request) { return; } this.pending.delete(event.data.id); if (event.data.error) { request.reject(new Error(event.data.error)); } else { request.resolve(event.data.contours ?? []); } }; this.worker.onerror = (event) => { const error = new Error( event.message || 'Contour worker failed.', ); this.pending.forEach((request) => request.reject(error)); this.pending.clear(); }; map.addLayer(this.layer); } async setFrame( frame: FlexpartGridFrame, ): Promise { this.validateFrame(frame); const version = ++this.frameVersion; const levels = this.resolveLevels(frame); const normalizedValues = this.prepareValues(frame); const contourGrid = this.prepareCyclicLongitudeGrid( frame.longitudes, frame.latitudes.length, normalizedValues, ); const rawContours = await this.calculateContours( contourGrid.longitudes.length, frame.latitudes.length, contourGrid.values, levels.slice(0, -1), ); if (version !== this.frameVersion) { throw new Error('A newer FLEXPART frame superseded this frame.'); } const projection = this.map.getView().getProjection(); this.projectedContours = rawContours.map((contour, index) => ({ color: this.colors[index], path: this.createProjectedPath( contour, contourGrid.longitudes, frame.latitudes, projection.getCode(), ), })); this.layerExtent = this.createLayerExtent( frame.longitudes, frame.latitudes, projection.getCode(), ); // A wrapped global layer must not be clipped to the canonical world. this.layer.setExtent(this.worldWidth ? undefined : this.layerExtent); this.layer.setVisible(true); this.source.changed(); this.map.render(); return { layerExtent: this.layerExtent, levels, validTime: frame.validTime, }; } getLayer(): ImageLayer { return this.layer; } setVisible(visible: boolean): void { this.layer.setVisible(visible); } setOpacity(opacity: number): void { this.layer.setOpacity(Math.max(0, Math.min(1, opacity))); } clear(): void { this.projectedContours = []; this.source.changed(); } destroy(): void { this.worker.terminate(); this.map.removeLayer(this.layer); } private validateFrame(frame: FlexpartGridFrame): void { const width = frame.longitudes.length; const height = frame.latitudes.length; if (width < 2 || height < 2) { throw new Error('Longitude and latitude need at least two points.'); } if (isNestedValues(frame.values)) { if (frame.values.length !== height || frame.values.some((row) => row.length !== width)) { throw new Error( `values must be a ${height}x${width} two-dimensional array.`, ); } } else if (frame.values.length !== width * height) { throw new Error( `values length must be ${width * height}, got ${frame.values.length}.`, ); } this.validateAxis(frame.longitudes, 'longitudes'); this.validateAxis(frame.latitudes, 'latitudes'); if (!Number.isFinite(frame.colorMin) || !Number.isFinite(frame.colorMax) || frame.colorMin <= 0 || frame.colorMax <= frame.colorMin) { throw new Error('Invalid colorMin/colorMax.'); } } private validateAxis(axis: readonly number[], name: string): void { let direction = 0; for (let index = 0; index < axis.length; index++) { if (!Number.isFinite(axis[index])) { throw new Error(`${name} contains a non-finite coordinate.`); } if (index === 0) { continue; } const difference = axis[index] - axis[index - 1]; const currentDirection = Math.sign(difference); if (currentDirection === 0 || (direction !== 0 && currentDirection !== direction)) { throw new Error(`${name} must be strictly monotonic.`); } direction = currentDirection; } } private resolveLevels(frame: FlexpartGridFrame): number[] { const levels = frame.levels ? Array.from(frame.levels) : createLogLevels(frame.colorMin, frame.colorMax); if (levels.length !== this.colors.length + 1) { throw new Error( `levels must contain ${this.colors.length + 1} values.`, ); } for (let index = 0; index < levels.length; index++) { if (!Number.isFinite(levels[index]) || levels[index] <= 0 || (index > 0 && levels[index] <= levels[index - 1])) { throw new Error('levels must be positive and strictly increasing.'); } } return levels; } private prepareValues(frame: FlexpartGridFrame): Float64Array { const width = frame.longitudes.length; const height = frame.latitudes.length; const values = new Float64Array(width * height); for (let row = 0; row < height; row++) { for (let column = 0; column < width; column++) { const sourceValue = isNestedValues(frame.values) ? frame.values[row][column] : frame.values[row * width + column]; const value = Number(sourceValue); // Positive values below colorMin still influence contour interpolation. // The first threshold keeps those cells transparent without log10 data. values[row * width + column] = Number.isFinite(value) && value > 0 ? value : Number.NaN; } } return values; } private prepareCyclicLongitudeGrid( longitudes: readonly number[], height: number, values: Float64Array, ): PreparedContourGrid { if (!this.wrapX || !isGlobalCyclicAxis(longitudes)) { return { longitudes, values }; } const sourceWidth = longitudes.length; const targetWidth = sourceWidth + 1; const direction = Math.sign( longitudes[sourceWidth - 1] - longitudes[0], ); const cyclicLongitudes = [ ...longitudes, longitudes[0] + direction * 360, ]; const cyclicValues = new Float64Array(targetWidth * height); for (let row = 0; row < height; row++) { const sourceOffset = row * sourceWidth; const targetOffset = row * targetWidth; cyclicValues.set( values.subarray(sourceOffset, sourceOffset + sourceWidth), targetOffset, ); cyclicValues[targetOffset + sourceWidth] = values[sourceOffset]; } return { longitudes: cyclicLongitudes, values: cyclicValues, }; } private calculateContours( width: number, height: number, values: Float64Array, thresholds: number[], ): Promise { const id = ++this.requestId; return new Promise((resolve, reject) => { this.pending.set(id, { resolve, reject }); const request: ContourWorkerRequest = { id, width, height, values, thresholds, }; this.worker.postMessage(request, [ values.buffer, ]); }); } private createProjectedPath( contour: SerializedContour, longitudes: readonly number[], latitudes: readonly number[], projectionCode: string, ): Path2D { const path = new Path2D(); contour.coordinates.forEach((polygon) => { polygon.forEach((ring) => { ring.forEach(([gridX, gridY], pointIndex) => { const lon = interpolateGridCoordinate(longitudes, gridX); const lat = interpolateGridCoordinate(latitudes, gridY); const projected = transform( [lon, clampLatitude(lat)], 'EPSG:4326', projectionCode, ); if (pointIndex === 0) { path.moveTo(projected[0], projected[1]); } else { path.lineTo(projected[0], projected[1]); } }); path.closePath(); }); }); return path; } private createLayerExtent( longitudes: readonly number[], latitudes: readonly number[], projectionCode: string, ): Extent { const lonMin = interpolateGridCoordinate(longitudes, 0); const lonMax = interpolateGridCoordinate(longitudes, longitudes.length); const latMin = interpolateGridCoordinate(latitudes, 0); const latMax = interpolateGridCoordinate(latitudes, latitudes.length); return transformExtent( [ Math.min(lonMin, lonMax), clampLatitude(Math.min(latMin, latMax)), Math.max(lonMin, lonMax), clampLatitude(Math.max(latMin, latMax)), ], 'EPSG:4326', projectionCode, ); } private drawCanvas( extent: Extent, size: readonly number[], ): HTMLCanvasElement { const canvas = document.createElement('canvas'); canvas.width = Math.max(1, Math.round(size[0])); canvas.height = Math.max(1, Math.round(size[1])); const context = canvas.getContext('2d'); if (!context || extent[2] <= extent[0] || extent[3] <= extent[1]) { return canvas; } context.clearRect(0, 0, canvas.width, canvas.height); const scaleX = canvas.width / (extent[2] - extent[0]); const scaleY = canvas.height / (extent[3] - extent[1]); context.setTransform( scaleX, 0, 0, -scaleY, -extent[0] * scaleX, extent[3] * scaleY, ); const worldOffsets = this.worldWidth ? [-this.worldWidth, 0, this.worldWidth] : [0]; // Lower thresholds are drawn first; higher thresholds cover them. this.projectedContours.forEach((contour) => { context.fillStyle = contour.color; worldOffsets.forEach((offset) => { context.save(); context.translate(offset, 0); context.fill(contour.path, 'evenodd'); context.restore(); }); }); context.setTransform(1, 0, 0, 1, 0, 0); return canvas; } } let defaultRenderer: FlexpartContourRenderer | undefined; /** Direct adapter for the existing olContainer map. */ export async function addFlexpartContourLayer( frame: FlexpartGridFrame, options: FlexpartContourLayerOptions = {}, ): Promise { const map = obtainViewer2D('olContainer') as OlMap | undefined; if (!map) { throw new Error('OpenLayers map olContainer is not initialized.'); } if (!defaultRenderer) { defaultRenderer = new FlexpartContourRenderer(map, options); } return defaultRenderer.setFrame(frame); } export function clearFlexpartContourLayer(): void { defaultRenderer?.clear(); } export function destroyFlexpartContourLayer(): void { defaultRenderer?.destroy(); defaultRenderer = undefined; } /** Move the existing olContainer view to a Pacific-centered longitude. */ // export function centerFlexpartMap( // longitude = 180, // latitude = 0, // duration = 0, // ): void { // const map = obtainViewer2D('olContainer') as OlMap | undefined; // if (!map) { // throw new Error('OpenLayers map olContainer is not initialized.'); // } // if (!Number.isFinite(longitude) // || !Number.isFinite(latitude) // || !Number.isFinite(duration) // || duration < 0) { // throw new Error('Invalid map center or animation duration.'); // } // const center = transform( // [longitude, clampLatitude(latitude)], // 'EPSG:4326', // map.getView().getProjection(), // ); // if (duration > 0) { // map.getView().animate({ center, duration }); // } // else { // map.getView().setCenter(center); // } // } function interpolateGridCoordinate( coordinates: readonly number[], gridPosition: number, ): number { const centerPosition = gridPosition - 0.5; const lastIndex = coordinates.length - 1; if (centerPosition <= 0) { return coordinates[0] + centerPosition * (coordinates[1] - coordinates[0]); } if (centerPosition >= lastIndex) { return coordinates[lastIndex] + (centerPosition - lastIndex) * (coordinates[lastIndex] - coordinates[lastIndex - 1]); } const lowerIndex = Math.floor(centerPosition); const ratio = centerPosition - lowerIndex; return coordinates[lowerIndex] + (coordinates[lowerIndex + 1] - coordinates[lowerIndex]) * ratio; } function clampLatitude(latitude: number): number { return Math.max( -WEB_MERCATOR_MAX_LATITUDE, Math.min(WEB_MERCATOR_MAX_LATITUDE, latitude), ); } function isNestedValues( values: FlexpartGridFrame['values'], ): values is readonly (readonly number[])[] { return Array.isArray(values) && values.length > 0 && Array.isArray(values[0]); } function isGlobalCyclicAxis(longitudes: readonly number[]): boolean { if (longitudes.length < 3) { return false; } const span = Math.abs( longitudes[longitudes.length - 1] - longitudes[0], ); const averageStep = span / (longitudes.length - 1); const coverage = span + averageStep; const tolerance = Math.max(1e-6, averageStep * 0.1); return Math.abs(coverage - 360) <= tolerance; }