fix: scale geostationary scan angles (rad) to PROJ metres before inverse projection

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"
     }
   ]
 }

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.
   }

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<string> {
   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,211 @@ 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: 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(discriminant)) / quadraticCoeff;
+  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),
+  ]);
+
+  // 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 satelliteOrbitalRadius = h + a to convert.
+  const satelliteOrbitalRadius = 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 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, yRaw, ...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 };
+}