Skip to content

fix(geotiff): read projection origin from all equivalent GeoKeys #671

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
kylebarron merged 1 commit into from
Sep 25, 2026
Merged

fix(geotiff): read projection origin from all equivalent GeoKeys #671

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion fixtures/geotiff-test-data
Open in desktop
Submodule geotiff-test-data updated 28 files
+6 −3 .github/workflows/ci.yml
+2 −2 README.md
+18 −11 generate_fixtures.py
+1,379 −181 pixi.lock
+5 −1 pixi.toml
+30 −0 rasterio_generated/fixtures/uint8_1band_sparse_no_nodata.py
+ − rasterio_generated/fixtures/uint8_1band_sparse_no_nodata.tif
+43 −0 rasterio_generated/fixtures/uint8_1band_sparse_no_nodata_info.md
+35 −0 rasterio_generated/fixtures/uint8_1band_sparse_nodata.py
+ − rasterio_generated/fixtures/uint8_1band_sparse_nodata.tif
+43 −0 rasterio_generated/fixtures/uint8_1band_sparse_nodata_info.md
+ − real_data/source-coop-dataforcanada/O2308000_7586000_cog.tif
+61 −0 real_data/source-coop-dataforcanada/O2308000_7586000_cog_info.md
+54 −0 real_data/source-coop-dataforcanada/README.md
+58 −0 tests/test_real_data.py
+7 −0 tifffile_generated/README.md
+1 −0 tifffile_generated/__init__.py
+1 −0 tifffile_generated/fixtures/__init__.py
+24 −0 tifffile_generated/fixtures/subifd_pyramid_bigtiff.py
+ − tifffile_generated/fixtures/subifd_pyramid_bigtiff.tif
+54 −0 tifffile_generated/fixtures/subifd_pyramid_bigtiff_info.md
+14 −0 tifffile_generated/fixtures/subifd_pyramid_classic.py
+ − tifffile_generated/fixtures/subifd_pyramid_classic.tif
+54 −0 tifffile_generated/fixtures/subifd_pyramid_classic_info.md
+33 −0 tifffile_generated/fixtures/trailing_directory.py
+ − tifffile_generated/fixtures/trailing_directory.tif
+52 −0 tifffile_generated/fixtures/trailing_directory_info.md
+73 −0 tifffile_generated/write_utils.py
163 changes: 90 additions & 73 deletions packages/geotiff/src/crs.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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,
),
],
};
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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,
),
],
};
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand All @@ -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),
],
};
}
Expand Down
54 changes: 54 additions & 0 deletions packages/geotiff/tests/crs.test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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", () => {
Expand Down Expand Up @@ -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,
},
});
});
});
Loading