From 91cd628dac32fb4a22ef8f22ebf8ce2c789cf684 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Fri, 25 Sep 2026 13:31:01 -0400 Subject: [PATCH] fix(geotiff): read projection origin from all equivalent GeoKeys Each projection method read its origin, scale factor and false easting/northing from a single hand-picked GeoKey. Writers disagree on which of the equivalent keys they use: GDAL writes ProjNatOrigin* for Oblique Stereographic, but we read ProjCenter*, so a user-defined New Brunswick Stereographic COG got an origin of (0, 0) and rendered near Null Island (source-cooperative/cog-viewer#40). Resolve each parameter the way libgeotiff's GTIFFetchProjParms does: from the first of its equivalent keys that is present. This also fixes the Stereographic scale factor, which GDAL writes to ProjScaleAtNatOrigin. Albers and LCC 2SP keep checking their false-origin keys first. Bumps geotiff-test-data for the O2308000_7586000_cog.tif fixture. Co-Authored-By: Claude Opus 5.5 (1M context) --- fixtures/geotiff-test-data | 2 +- packages/geotiff/src/crs.ts | 163 ++++++++++++++++------------- packages/geotiff/tests/crs.test.ts | 54 ++++++++++ 3 files changed, 145 insertions(+), 74 deletions(-) diff --git a/fixtures/geotiff-test-data b/fixtures/geotiff-test-data index 6f403786..f2980131 160000 --- a/fixtures/geotiff-test-data +++ b/fixtures/geotiff-test-data @@ -1 +1 @@ -Subproject commit 6f4037861ed474d39ef81d02207b3f252c897d7d +Subproject commit f2980131906d2f762f1c954d693a56cd805a0ccb diff --git a/packages/geotiff/src/crs.ts b/packages/geotiff/src/crs.ts index 6dce837c..4354a599 100644 --- a/packages/geotiff/src/crs.ts +++ b/packages/geotiff/src/crs.ts @@ -243,6 +243,23 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { throw new Error("User-defined projected CRS requires projMethod"); } + // Writers disagree on which of several equivalent keys a method stores its + // origin in: GDAL writes ProjNatOrigin* for Oblique Stereographic but + // ProjCenter* for Stereographic, for example. Like libgeotiff, take each + // parameter from the first of its equivalent keys that is present. + // https://github.com/OSGeo/libgeotiff/blob/75cfca539667c6483e796b24b78ecef72e311a1a/libgeotiff/geo_normalize.c#L1649 + const originLat = + gkd.projNatOriginLat ?? gkd.projFalseOriginLat ?? gkd.projCenterLat; + const originLong = + gkd.projNatOriginLong ?? gkd.projFalseOriginLong ?? gkd.projCenterLong; + const originScale = gkd.projScaleAtNatOrigin ?? gkd.projScaleAtCenter; + const falseEasting = + gkd.projFalseEasting ?? gkd.projCenterEasting ?? gkd.projFalseOriginEasting; + const falseNorthing = + gkd.projFalseNorthing ?? + gkd.projCenterNorthing ?? + gkd.projFalseOriginNorthing; + const angular = ( name: string, value: number | null, @@ -285,11 +302,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - scale("Scale factor at natural origin", gkd.projScaleAtNatOrigin), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + scale("Scale factor at natural origin", originScale), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -303,13 +320,13 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of projection centre", gkd.projCenterLat), - angular("Longitude of projection centre", gkd.projCenterLong), + angular("Latitude of projection centre", originLat), + angular("Longitude of projection centre", originLong), angular("Azimuth of initial line", gkd.projAzimuthAngle), angular("Angle from Rectified to Skew Grid", gkd.projAzimuthAngle), - scale("Scale factor on initial line", gkd.projScaleAtCenter), - linear("Easting at projection centre", gkd.projCenterEasting), - linear("Northing at projection centre", gkd.projCenterNorthing), + scale("Scale factor on initial line", originScale), + linear("Easting at projection centre", falseEasting), + linear("Northing at projection centre", falseNorthing), ], }; } @@ -320,11 +337,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - scale("Scale factor at natural origin", gkd.projScaleAtNatOrigin), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + scale("Scale factor at natural origin", originScale), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -337,21 +354,21 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { parameters: [ angular( "Latitude of false origin", - gkd.projFalseOriginLat ?? gkd.projNatOriginLat, + gkd.projFalseOriginLat ?? originLat, ), angular( "Longitude of false origin", - gkd.projFalseOriginLong ?? gkd.projNatOriginLong, + gkd.projFalseOriginLong ?? originLong, ), angular("Latitude of 1st standard parallel", gkd.projStdParallel1), angular("Latitude of 2nd standard parallel", gkd.projStdParallel2), linear( "Easting at false origin", - gkd.projFalseOriginEasting ?? gkd.projFalseEasting, + gkd.projFalseOriginEasting ?? falseEasting, ), linear( "Northing at false origin", - gkd.projFalseOriginNorthing ?? gkd.projFalseNorthing, + gkd.projFalseOriginNorthing ?? falseNorthing, ), ], }; @@ -363,11 +380,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - scale("Scale factor at natural origin", gkd.projScaleAtNatOrigin), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + scale("Scale factor at natural origin", originScale), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -378,10 +395,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projCenterLat), - angular("Longitude of natural origin", gkd.projCenterLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -394,21 +411,21 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { parameters: [ angular( "Latitude of false origin", - gkd.projFalseOriginLat ?? gkd.projNatOriginLat, + gkd.projFalseOriginLat ?? originLat, ), angular( "Longitude of false origin", - gkd.projFalseOriginLong ?? gkd.projNatOriginLong, + gkd.projFalseOriginLong ?? originLong, ), angular("Latitude of 1st standard parallel", gkd.projStdParallel1), angular("Latitude of 2nd standard parallel", gkd.projStdParallel2), linear( "Easting at false origin", - gkd.projFalseOriginEasting ?? gkd.projFalseEasting, + gkd.projFalseOriginEasting ?? falseEasting, ), linear( "Northing at false origin", - gkd.projFalseOriginNorthing ?? gkd.projFalseNorthing, + gkd.projFalseOriginNorthing ?? falseNorthing, ), ], }; @@ -420,10 +437,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projCenterLat), - angular("Longitude of natural origin", gkd.projCenterLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -434,11 +451,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projCenterLat), - angular("Longitude of natural origin", gkd.projCenterLong), - scale("Scale factor at natural origin", gkd.projScaleAtCenter), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + scale("Scale factor at natural origin", originScale), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -451,14 +468,14 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { parameters: [ angular( "Latitude of standard parallel", - gkd.projNatOriginLat ?? gkd.projStdParallel1, + originLat ?? gkd.projStdParallel1, ), angular( "Longitude of origin", - gkd.projStraightVertPoleLong ?? gkd.projNatOriginLong, + gkd.projStraightVertPoleLong ?? originLong, ), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -469,11 +486,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projCenterLat), - angular("Longitude of natural origin", gkd.projCenterLong), - scale("Scale factor at natural origin", gkd.projScaleAtCenter), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + scale("Scale factor at natural origin", originScale), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -486,11 +503,11 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { parameters: [ angular( "Latitude of 1st standard parallel", - gkd.projStdParallel1 ?? gkd.projCenterLat, + gkd.projStdParallel1 ?? originLat, ), - angular("Longitude of natural origin", gkd.projCenterLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -501,10 +518,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -515,10 +532,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -529,9 +546,9 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Longitude of natural origin", gkd.projCenterLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -542,10 +559,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projCenterLat), - angular("Longitude of natural origin", gkd.projCenterLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } @@ -556,10 +573,10 @@ function _buildConversion(gkd: GeoKeyDirectory): ProjJsonConversion { name, method: { name }, parameters: [ - angular("Latitude of natural origin", gkd.projNatOriginLat), - angular("Longitude of natural origin", gkd.projNatOriginLong), - linear("False easting", gkd.projFalseEasting), - linear("False northing", gkd.projFalseNorthing), + angular("Latitude of natural origin", originLat), + angular("Longitude of natural origin", originLong), + linear("False easting", falseEasting), + linear("False northing", falseNorthing), ], }; } diff --git a/packages/geotiff/tests/crs.test.ts b/packages/geotiff/tests/crs.test.ts index bc322d21..046f0c79 100644 --- a/packages/geotiff/tests/crs.test.ts +++ b/packages/geotiff/tests/crs.test.ts @@ -1,5 +1,6 @@ import { parseWkt } from "@developmentseed/proj"; import { describe, expect, it } from "vitest"; +import { crsFromGeoKeys } from "../src/crs.js"; import { loadGeoTIFF } from "./helpers.js"; describe("test CRS", () => { @@ -131,3 +132,56 @@ describe("test GeoKey CRS parsing", () => { expect(proj.a ?? proj.datum?.a).toBe(6378137); }); }); + +/** + * Parameters of New Brunswick Stereographic (equivalent to EPSG:2953), in the + * order the Oblique Stereographic and Stereographic conversions emit them. + */ +const NEW_BRUNSWICK_STEREOGRAPHIC_PARAMETERS = [ + { name: "Latitude of natural origin", value: 46.5, unit: "degree" }, + { name: "Longitude of natural origin", value: -66.5, unit: "degree" }, + { name: "Scale factor at natural origin", value: 0.999912, unit: "unity" }, + { name: "False easting", value: 2500000, unit: "metre" }, + { name: "False northing", value: 7500000, unit: "metre" }, +]; + +describe("user-defined projection parameters", () => { + it("reads the Oblique Stereographic origin from ProjNatOrigin* keys", async () => { + // https://github.com/source-cooperative/cog-viewer/issues/40 + const geotiff = await loadGeoTIFF( + "O2308000_7586000_cog", + "source-coop-dataforcanada", + ); + + expect(geotiff.crs).toMatchObject({ + conversion: { + method: { name: "Oblique Stereographic" }, + parameters: NEW_BRUNSWICK_STEREOGRAPHIC_PARAMETERS, + }, + }); + }); + + it("reads the Stereographic scale factor from ProjScaleAtNatOriginGeoKey", async () => { + // GDAL writes CT_Stereographic with the origin in ProjCenter{Lat,Long} but + // the scale factor in ProjScaleAtNatOrigin. + const { gkd } = await loadGeoTIFF( + "O2308000_7586000_cog", + "source-coop-dataforcanada", + ); + const crs = crsFromGeoKeys({ + ...gkd, + projMethod: 14, + projCenterLat: gkd.projNatOriginLat, + projCenterLong: gkd.projNatOriginLong, + projNatOriginLat: null, + projNatOriginLong: null, + }); + + expect(crs).toMatchObject({ + conversion: { + method: { name: "Stereographic" }, + parameters: NEW_BRUNSWICK_STEREOGRAPHIC_PARAMETERS, + }, + }); + }); +});