From cda1d25d90991b8e2df2ee245f72facb7aac59d9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 28 May 2026 03:51:53 +0000 Subject: [PATCH 1/5] Initial plan From fe953674216f8df9be7ca2842c513ef01786de33 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 28 May 2026 04:01:45 +0000 Subject: [PATCH 2/5] feat: add geostationary projection support (Himawari-8) --- public/static/catalog.json | 10 ++ src/lib/data/gridTypeDetector.ts | 9 ++ src/lib/data/zarrUtils.ts | 194 ++++++++++++++++++++++++++++++- src/ui/grids/Curvilinear.vue | 83 +++++++++---- 4 files changed, 272 insertions(+), 24 deletions(-) diff --git a/public/static/catalog.json b/public/static/catalog.json index b7f7b71..109f824 100644 --- a/public/static/catalog.json +++ b/public/static/catalog.json @@ -441,6 +441,16 @@ "layout": "grouped", "grid": "utm", "crs": "EPSG:32637" + }, + { + "title": "Himawari-8 AHI Full Disk Reflectance", + "url": "https://data.source.coop/bkr/geo/himawari_1000m.icechunk", + "format": "Zarr v3", + "access": "icechunk", + "convention": null, + "layout": "simple", + "grid": "geostationary", + "crs": "geos" } ] } diff --git a/src/lib/data/gridTypeDetector.ts b/src/lib/data/gridTypeDetector.ts index 6c76346..9cbf7a1 100644 --- a/src/lib/data/gridTypeDetector.ts +++ b/src/lib/data/gridTypeDetector.ts @@ -4,6 +4,7 @@ import { ZarrDataManager } from "./ZarrDataManager.ts"; import { getCRSStringForXYVariable, getLatLonVariableInfo, + isGeostationaryCRS, isLatitudeName, isLongitudeName, isPolarStereographicCRS, @@ -125,6 +126,11 @@ async function determineGridTypeFromCRS( if (checkRegularRotatedGrid(crs)) { return GRID_TYPES.REGULAR_ROTATED; } + // Geostationary datasets are routed to CURVILINEAR so that + // computeGeostatLatLon2D can produce proper 2-D lat/lon arrays. + if (crs.attrs?.grid_mapping_name === "geostationary") { + return GRID_TYPES.CURVILINEAR; + } // Polar stereographic datasets are routed to CURVILINEAR so that // computePolarStereoLatLon2D can produce proper 2-D lat/lon arrays. if (crs.attrs?.grid_mapping_name === "polar_stereographic") { @@ -144,6 +150,9 @@ async function determineGridTypeFromCRS( if (isPolarStereographicCRS(crsStr)) { return GRID_TYPES.CURVILINEAR; } + if (isGeostationaryCRS(crsStr)) { + return GRID_TYPES.CURVILINEAR; + } } catch { // No CRS info available. } diff --git a/src/lib/data/zarrUtils.ts b/src/lib/data/zarrUtils.ts index 9aa0597..3764fe4 100644 --- a/src/lib/data/zarrUtils.ts +++ b/src/lib/data/zarrUtils.ts @@ -105,11 +105,11 @@ export function isLongitudeName(name: string) { } export function isXName(name: string) { - return name === "x"; + return name === "x" || name === "x_geostationary"; } export function isYName(name: string) { - return name === "y"; + return name === "y" || name === "y_geostationary"; } export function isLatitudeVariable(name: string, attrs: unknown) { @@ -619,6 +619,9 @@ export async function getCRSStringForXYVariable( ): Promise { try { const crs = await ZarrDataManager.getCRSInfo(datasources, currentVarname); + if (crs.attrs?.grid_mapping_name === "geostationary") { + return "geostationary"; + } if (crs.attrs?.grid_mapping_name === "polar_stereographic") { return "polar_stereographic"; } @@ -816,3 +819,190 @@ export async function computePolarStereoLatLon2D( return { latitudes2D, longitudes2D, ny, nx, isNorthPole }; } + +// --------------------------------------------------------------------------- +// Geostationary projection support (PROJ `geos`) +// --------------------------------------------------------------------------- + +/** + * Returns true when the string describes a geostationary projection. + */ +export function isGeostationaryCRS(str: string): boolean { + const lower = str.toLowerCase(); + return lower.includes("geostationary") || lower.includes("+proj=geos"); +} + +/** + * Extract the projection parameters for a geostationary dataset from the + * CRS variable attached to `currentVarname`. + * + * @returns `{ lon0, h, a }` where `lon0` is the sub-satellite longitude (°), + * `h` is the perspective point height above the surface (m), and + * `a` is the semi-major axis (m). + */ +export async function getGeostatCRSParams( + datasources: TSources, + currentVarname: string +): Promise<{ lon0: number; h: number; a: number }> { + const crs = await ZarrDataManager.getCRSInfo(datasources, currentVarname); + const lon0 = Number( + crs.attrs?.longitude_of_projection_origin ?? crs.attrs?.lon_0 ?? 0 + ); + const h = Number( + crs.attrs?.perspective_point_height ?? crs.attrs?.h ?? 35785831 // nominal geostationary orbit height + ); + const a = Number( + crs.attrs?.semi_major_axis ?? crs.attrs?.a ?? 6378137 // WGS-84 semi-major + ); + return { lon0, h, a }; +} + +/** + * Inverse geostationary (PROJ `geos`) projection — spherical Earth. + * + * Converts a single (x, y) point in the projected coordinate system + * (units: metres) to geographic (lat, lon) in degrees. + * + * Derivation uses PROJ4 `geos` forward formula: + * Let rg = H/a, rg1 = h/a = rg - 1 + * x/a = rg * cos(φ) * sin(λ) / D, y/a = sin(φ) / D + * where D = rg1 + cos(φ)*cos(λ) + * + * Substituting C = cos(φ)*cos(λ), S = sin(φ), T = cos(φ)*sin(λ), + * with C² + S² + T² = 1 gives a quadratic in D, solved for the + * root that maps (0,0) to the sub-satellite point. + * + * @returns `{ lat, lon }` in degrees, or `null` when the point falls + * outside the visible Earth disk. + */ +function invGeostatPoint( + x: number, + y: number, + h: number, + a: number, + lon0: number +): { lat: number; lon: number } | null { + const rg = (h + a) / a; // satellite orbital radius / Earth radius + const rg1 = h / a; // = rg - 1 + + const xn = x / a; + const yn = y / a; + + // Quadratic coefficients: a_q * D² - 2*rg1 * D + (rg1² - 1) = 0 + const aq = 1 + yn * yn + (xn * xn) / (rg * rg); + const disc = rg1 * rg1 - aq * (rg1 * rg1 - 1); + if (disc < 0) { + return null; // outside Earth disk + } + + const D = (rg1 + Math.sqrt(disc)) / aq; + const C = D - rg1; // = cos(φ)*cos(λ) + const S = yn * D; // = sin(φ) + const T = (xn * D) / rg; // = cos(φ)*sin(λ) + + const lat = (Math.atan2(S, Math.sqrt(C * C + T * T)) * 180) / Math.PI; + let lon = (Math.atan2(T, C) * 180) / Math.PI + lon0; + if (lon > 180) { + lon -= 360; + } + if (lon <= -180) { + lon += 360; + } + + return { lat, lon }; +} + +/** + * Load the 1-D x and y coordinate arrays for a geostationary grid along with + * the CRS projection parameters. + */ +async function loadGeostatCoords( + datasources: TSources, + currentVarname: string +): Promise<{ + xRaw: Float32Array; + yRaw: Float32Array; + lon0: number; + h: number; + a: number; +}> { + const dimensions = await ZarrDataManager.getDimensionNames( + datasources, + currentVarname + ); + const xDimName = dimensions.find(isXName) ?? "x_geostationary"; + const yDimName = dimensions.find(isYName) ?? "y_geostationary"; + + const xRef = ZarrDataManager.resolveVariableReference( + datasources, + currentVarname, + xDimName + ); + const yRef = ZarrDataManager.resolveVariableReference( + datasources, + currentVarname, + yDimName + ); + + const [xArray, yArray] = await Promise.all([ + ZarrDataManager.getVariableInfo(xRef.datasource, xRef.variable), + ZarrDataManager.getVariableInfo(yRef.datasource, yRef.variable), + ]); + + const [xData, yData, params] = await Promise.all([ + ZarrDataManager.getVariableDataFromArray(xArray), + ZarrDataManager.getVariableDataFromArray(yArray), + getGeostatCRSParams(datasources, currentVarname), + ]); + + return { + xRaw: castDataVarToFloat32(xData.data), + yRaw: castDataVarToFloat32(yData.data), + ...params, + }; +} + +/** + * Compute geographic (lat/lon) 2-D coordinate arrays for a geostationary + * x/y grid by applying the inverse PROJ `geos` projection to every + * (xᵢ, yⱼ) grid point. + * + * The returned flat arrays have length `ny × nx` and are laid out in + * row-major order (j-major, i-minor), matching the data array layout + * expected by the Curvilinear grid renderer. Points outside the visible + * Earth disk are stored as NaN. + */ +export async function computeGeostatLatLon2D( + datasources: TSources, + currentVarname: string +): Promise<{ + latitudes2D: Float64Array; + longitudes2D: Float64Array; + ny: number; + nx: number; +}> { + const { xRaw, yRaw, lon0, h, a } = await loadGeostatCoords( + datasources, + currentVarname + ); + const nx = xRaw.length; + const ny = yRaw.length; + + const latitudes2D = new Float64Array(ny * nx); + const longitudes2D = new Float64Array(ny * nx); + + for (let j = 0; j < ny; j++) { + for (let i = 0; i < nx; i++) { + const result = invGeostatPoint(xRaw[i], yRaw[j], h, a, lon0); + if (result !== null) { + latitudes2D[j * nx + i] = result.lat; + longitudes2D[j * nx + i] = result.lon; + } else { + latitudes2D[j * nx + i] = NaN; + longitudes2D[j * nx + i] = NaN; + } + } + } + + return { latitudes2D, longitudes2D, ny, nx }; +} diff --git a/src/ui/grids/Curvilinear.vue b/src/ui/grids/Curvilinear.vue index 3624b5a..b7e49d7 100755 --- a/src/ui/grids/Curvilinear.vue +++ b/src/ui/grids/Curvilinear.vue @@ -15,13 +15,16 @@ import { ZarrDataManager } from "@/lib/data/ZarrDataManager.ts"; import { applyDisplayTransformToData, castDataVarToFloat32, + computeGeostatLatLon2D, computePolarStereoLatLon2D, createMissingOrFillPredicate, getDataBounds, getCRSStringForXYVariable, + getGeostatCRSParams, getLatLonData, getMissingAndFillValues, getPolarStereoCRSParams, + isGeostationaryCRS, isPolarStereographicCRS, mapMissingAndFillToNaN, } from "@/lib/data/zarrUtils.ts"; @@ -70,6 +73,8 @@ const updatingData = ref(false); /** True when the loaded dataset uses a polar stereographic CRS. */ const isPolarStereoData = ref(false); +/** True when the loaded dataset uses a geostationary CRS. */ +const isGeostatData = ref(false); /** Aspect ratio (width / height = nx / ny) for the polar stereo grid canvas. */ const polarAspectRatio = ref(1); /** Inline style for the canvas box — sets aspect-ratio when polar data is loaded. */ @@ -193,8 +198,10 @@ const colormapMaterial = computed(() => { async function datasourceUpdate() { clearHoverLookup(); isPolarStereoData.value = false; + isGeostatData.value = false; if (props.datasources !== undefined) { - // Detect polar stereographic CRS and auto-configure projection/mask. + // Detect polar stereographic or geostationary CRS and auto-configure + // projection and mask. try { const crsStr = await getCRSStringForXYVariable( props.datasources, @@ -212,10 +219,22 @@ async function datasourceUpdate() { varnameSelector.value ); store.projectionCenter = { lat: isNorthPole ? 90 : -90, lon: 0 }; + } else if (isGeostationaryCRS(crsStr)) { + isGeostatData.value = true; + // Use nearside-perspective projection centred on the sub-satellite + // point so the full disk is immediately visible. + store.projectionMode = PROJECTION_TYPES.NEARSIDE_PERSPECTIVE; + // The global land/sea mask is not meaningful for a geostationary domain. + store.landSeaMaskChoice = LAND_SEA_MASK_MODES.OFF; + const { lon0 } = await getGeostatCRSParams( + props.datasources, + varnameSelector.value + ); + store.projectionCenter = { lat: 0, lon: lon0 }; } } catch { // CRS lookup may fail for datasets without a CRS variable or group-level - // projection attributes. Treat as a non-polar curvilinear grid. + // projection attributes. Treat as a standard curvilinear grid. } await Promise.all([getData()]); updateLandSeaMask(); @@ -225,45 +244,65 @@ async function datasourceUpdate() { const BATCH_SIZE = 30; +type TLatLon2DResult = { + latitudesData: Float64Array; + longitudesData: Float64Array; + nj: number; + ni: number; +}; + +/** Reshape a `computePolarStereoLatLon2D` / `computeGeostatLatLon2D` result. */ +function reshapeProjectedLatLon(result: { + latitudes2D: Float64Array; + longitudes2D: Float64Array; + ny: number; + nx: number; +}): TLatLon2DResult { + polarAspectRatio.value = result.nx / result.ny; + return { + latitudesData: result.latitudes2D, + longitudesData: result.longitudes2D, + nj: result.ny, + ni: result.nx, + }; +} + /** * Resolve the 2-D lat/lon coordinate arrays for a grid variable. * Returns proper geographic coordinates whether the variable uses named - * lat/lon arrays (standard curvilinear) or x/y arrays with a polar - * stereographic CRS (computed via inverse projection). + * lat/lon arrays (standard curvilinear), x/y arrays with a polar + * stereographic CRS, or x/y arrays with a geostationary CRS. */ async function resolveLatLon2D( datavar: zarr.Array -): Promise<{ - latitudesData: Float64Array; - longitudesData: Float64Array; - nj: number; - ni: number; -}> { - // Detect polar stereographic CRS. +): Promise { + // Detect projected CRS type. let isPolarStereo = false; + let isGeostat = false; try { const crsStr = await getCRSStringForXYVariable( props.datasources!, varnameSelector.value ); isPolarStereo = isPolarStereographicCRS(crsStr); + isGeostat = !isPolarStereo && isGeostationaryCRS(crsStr); } catch { // No CRS info available — treat as regular curvilinear. } if (isPolarStereo) { - const result = await computePolarStereoLatLon2D( - props.datasources!, - varnameSelector.value + return reshapeProjectedLatLon( + await computePolarStereoLatLon2D( + props.datasources!, + varnameSelector.value + ) + ); + } + + if (isGeostat) { + return reshapeProjectedLatLon( + await computeGeostatLatLon2D(props.datasources!, varnameSelector.value) ); - // Update aspect ratio to match actual grid dimensions (nx / ny). - polarAspectRatio.value = result.nx / result.ny; - return { - latitudesData: result.latitudes2D, - longitudesData: result.longitudes2D, - nj: result.ny, - ni: result.nx, - }; } const { latitudes, longitudes } = await getLatLonData( From bf37f553117c58f29f286d80428777c776e6c451 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 28 May 2026 04:03:34 +0000 Subject: [PATCH 3/5] fix: rename aq/disc to quadraticCoeff/discriminant; drop redundant CRS check --- src/lib/data/zarrUtils.ts | 10 +++++----- src/ui/grids/Curvilinear.vue | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/lib/data/zarrUtils.ts b/src/lib/data/zarrUtils.ts index 3764fe4..c4647ae 100644 --- a/src/lib/data/zarrUtils.ts +++ b/src/lib/data/zarrUtils.ts @@ -888,14 +888,14 @@ function invGeostatPoint( const xn = x / a; const yn = y / a; - // Quadratic coefficients: a_q * D² - 2*rg1 * D + (rg1² - 1) = 0 - const aq = 1 + yn * yn + (xn * xn) / (rg * rg); - const disc = rg1 * rg1 - aq * (rg1 * rg1 - 1); - if (disc < 0) { + // Quadratic coefficients: quadraticCoeff * D² - 2*rg1 * D + (rg1² - 1) = 0 + const quadraticCoeff = 1 + yn * yn + (xn * xn) / (rg * rg); + const discriminant = rg1 * rg1 - quadraticCoeff * (rg1 * rg1 - 1); + if (discriminant < 0) { return null; // outside Earth disk } - const D = (rg1 + Math.sqrt(disc)) / aq; + const D = (rg1 + Math.sqrt(discriminant)) / quadraticCoeff; const C = D - rg1; // = cos(φ)*cos(λ) const S = yn * D; // = sin(φ) const T = (xn * D) / rg; // = cos(φ)*sin(λ) diff --git a/src/ui/grids/Curvilinear.vue b/src/ui/grids/Curvilinear.vue index b7e49d7..c1b22eb 100755 --- a/src/ui/grids/Curvilinear.vue +++ b/src/ui/grids/Curvilinear.vue @@ -285,7 +285,7 @@ async function resolveLatLon2D( varnameSelector.value ); isPolarStereo = isPolarStereographicCRS(crsStr); - isGeostat = !isPolarStereo && isGeostationaryCRS(crsStr); + isGeostat = isGeostationaryCRS(crsStr); } catch { // No CRS info available — treat as regular curvilinear. } From 77224692eba339f60d8a57f202831d4f18ec7bf2 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 28 May 2026 19:43:21 +0000 Subject: [PATCH 4/5] fix: scale geostationary scan angles (rad) to PROJ metres before inverse projection --- src/lib/data/zarrUtils.ts | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/lib/data/zarrUtils.ts b/src/lib/data/zarrUtils.ts index c4647ae..06c490a 100644 --- a/src/lib/data/zarrUtils.ts +++ b/src/lib/data/zarrUtils.ts @@ -955,9 +955,24 @@ async function loadGeostatCoords( getGeostatCRSParams(datasources, currentVarname), ]); + // CF-convention geostationary datasets (e.g. Himawari-8 AHI, GOES-R ABI) + // store x_geostationary / y_geostationary as scan angles in radians. + // The PROJ `geos` inverse formula in invGeostatPoint expects coordinates + // in metres. Multiply by H = h + a (satellite orbital radius) to convert. + const H = params.h + params.a; + const xUnitsRad = + hasUnits(xArray.attrs) && + xArray.attrs.units.toLowerCase().startsWith("rad"); + const yUnitsRad = + hasUnits(yArray.attrs) && + yArray.attrs.units.toLowerCase().startsWith("rad"); + + const xRawBase = castDataVarToFloat32(xData.data); + const yRawBase = castDataVarToFloat32(yData.data); + return { - xRaw: castDataVarToFloat32(xData.data), - yRaw: castDataVarToFloat32(yData.data), + xRaw: xUnitsRad ? xRawBase.map((v) => v * H) : xRawBase, + yRaw: yUnitsRad ? yRawBase.map((v) => v * H) : yRawBase, ...params, }; } From bebd7da6ada36c2cb2630c09b40ca2bbda66fdde Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 28 May 2026 19:45:12 +0000 Subject: [PATCH 5/5] fix: rename H to satelliteOrbitalRadius; scale coordinate arrays in-place --- src/lib/data/zarrUtils.ts | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/lib/data/zarrUtils.ts b/src/lib/data/zarrUtils.ts index 06c490a..79d643c 100644 --- a/src/lib/data/zarrUtils.ts +++ b/src/lib/data/zarrUtils.ts @@ -958,8 +958,8 @@ async function loadGeostatCoords( // CF-convention geostationary datasets (e.g. Himawari-8 AHI, GOES-R ABI) // store x_geostationary / y_geostationary as scan angles in radians. // The PROJ `geos` inverse formula in invGeostatPoint expects coordinates - // in metres. Multiply by H = h + a (satellite orbital radius) to convert. - const H = params.h + params.a; + // in metres. Multiply by satelliteOrbitalRadius = h + a to convert. + const satelliteOrbitalRadius = params.h + params.a; const xUnitsRad = hasUnits(xArray.attrs) && xArray.attrs.units.toLowerCase().startsWith("rad"); @@ -967,14 +967,20 @@ async function loadGeostatCoords( hasUnits(yArray.attrs) && yArray.attrs.units.toLowerCase().startsWith("rad"); - const xRawBase = castDataVarToFloat32(xData.data); - const yRawBase = castDataVarToFloat32(yData.data); + const xRaw = castDataVarToFloat32(xData.data); + const yRaw = castDataVarToFloat32(yData.data); + if (xUnitsRad) { + for (let i = 0; i < xRaw.length; i++) { + xRaw[i] *= satelliteOrbitalRadius; + } + } + if (yUnitsRad) { + for (let i = 0; i < yRaw.length; i++) { + yRaw[i] *= satelliteOrbitalRadius; + } + } - return { - xRaw: xUnitsRad ? xRawBase.map((v) => v * H) : xRawBase, - yRaw: yUnitsRad ? yRawBase.map((v) => v * H) : yRawBase, - ...params, - }; + return { xRaw, yRaw, ...params }; } /**