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..79d643c 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,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 }; +} diff --git a/src/ui/grids/Curvilinear.vue b/src/ui/grids/Curvilinear.vue index 3624b5a..c1b22eb 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 = 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(