feat(geotiff): generate TileMatrixSet from GeoTIFF instance

packages/affine/src/affine.ts

@@ -83,3 +83,33 @@ export function invert([sa, sb, sc, sd, se, sf]: Affine): Affine {
 
   return [ra, rb, -sc * ra - sf * rb, rd, re, -sc * rd - sf * re];
 }
+
+/** Get the 'a' component of an Affine transform. */
+export function a(affine: Affine): number {
+  return affine[0];
+}
+
+/** Get the 'b' component of an Affine transform. */
+export function b(affine: Affine): number {
+  return affine[1];
+}
+
+/** Get the 'c' component of an Affine transform. */
+export function c(affine: Affine): number {
+  return affine[2];
+}
+
+/** Get the 'd' component of an Affine transform. */
+export function d(affine: Affine): number {
+  return affine[3];
+}
+
+/** Get the 'e' component of an Affine transform. */
+export function e(affine: Affine): number {
+  return affine[4];
+}
+
+/** Get the 'f' component of an Affine transform. */
+export function f(affine: Affine): number {
+  return affine[5];
+}

packages/geotiff/package.json

@@ -47,6 +47,8 @@
   },
   "dependencies": {
     "@cogeotiff/core": "^9.1.2",
-    "@developmentseed/affine": "workspace:^"
+    "@developmentseed/affine": "workspace:^",
+    "@developmentseed/morecantile": "workspace:^",
+    "uuid": "^13.0.0"
   }
 }

packages/geotiff/src/geotiff.ts

@@ -147,16 +147,6 @@ export class GeoTIFF {
     return this._crs;
   }
 
-  /** EPSG code from GeoTIFF tags, or null if not set.
-   *
-   * See also {@link GeoTIFF.crs} for the full PROJJSON definition, which should
-   * always be available, even when an EPSG code is not explicitly stored in the
-   * tags.
-   */
-  get epsg(): number | null {
-    return this.image.epsg;
-  }
-
   /** Image width in pixels. */
   get width(): number {
     return this.image.size.width;

packages/geotiff/src/index.ts

@@ -5,3 +5,4 @@ export { decode, registry } from "./decode/api.js";
 export { GeoTIFF } from "./geotiff.js";
 export { Overview } from "./overview.js";
 export type { Tile } from "./tile.js";
+export { generateTileMatrixSet } from "./tile-matrix-set.js";

packages/geotiff/src/tile-matrix-set.ts

@@ -0,0 +1,166 @@
+import type { Affine } from "@developmentseed/affine";
+import * as affine from "@developmentseed/affine";
+import type {
+  BoundingBox,
+  CRS,
+  TileMatrix,
+  TileMatrixSet,
+} from "@developmentseed/morecantile";
+import { metersPerUnit } from "@developmentseed/morecantile";
+import { v4 as uuidv4 } from "uuid";
+import type { ProjJson } from "./crs.js";
+import type { GeoTIFF } from "./geotiff.js";
+
+/**
+ * A minimal projection definition compatible with what wkt-parser returns.
+ *
+ * This type extracts only the partial properties we need from the full
+ * wkt-parser output.
+ */
+interface ProjectionDefinition {
+  datum?: {
+    /** Semi-major axis of the ellipsoid. */
+    a: number;
+  };
+  a?: number;
+  to_meter?: number;
+  units?: string;
+}
+
+const SCREEN_PIXEL_SIZE = 0.28e-3;
+
+function buildCrs(crs: number | ProjJson): CRS {
+  if (typeof crs === "number") {
+    return {
+      uri: `http://www.opengis.net/def/crs/EPSG/0/${crs}`,
+    };
+  }
+
+  // @ts-expect-error - typing issues between different projjson definitions.
+  return {
+    wkt: crs,
+  };
+}
+
+/**
+ * Build a TileMatrix entry for a single resolution level.
+ */
+function buildTileMatrix(
+  id: string,
+  transform: Affine,
+  mpu: number,
+  cornerOfOrigin: "bottomLeft" | "topLeft",
+  tileWidth: number,
+  tileHeight: number,
+  width: number,
+  height: number,
+): TileMatrix {
+  return {
+    id,
+    scaleDenominator: (affine.a(transform) * mpu) / SCREEN_PIXEL_SIZE,
+    cellSize: affine.a(transform),
+    cornerOfOrigin,
+    pointOfOrigin: [affine.c(transform), affine.f(transform)],
+    tileWidth,
+    tileHeight,
+    matrixWidth: Math.ceil(width / tileWidth),
+    matrixHeight: Math.ceil(height / tileHeight),
+  };
+}
+
+/**
+ * Generate a Tile Matrix Set from a GeoTIFF file.
+ *
+ * Produces one TileMatrix per overview (coarsest first) plus a final entry
+ * for the full-resolution level. The GeoTIFF must be tiled.
+ *
+ * This requires a crs definition that includes a `units` property, so that we
+ * can convert pixel sizes to physical screen units. Use [`wkt-parser`] to parse
+ * a WKT string or PROJJSON object, then pass the result as the `crs` argument.
+ *
+ * [`wkt-parser`]: https://github.com/proj4js/wkt-parser
+ *
+ * @see https://docs.ogc.org/is/17-083r4/17-083r4.html
+ */
+export function generateTileMatrixSet(
+  geotiff: GeoTIFF,
+  crs: ProjectionDefinition,
+  { id = uuidv4() }: { id?: string } = {},
+): TileMatrixSet {
+  const bbox = geotiff.bbox;
+  const tr = geotiff.transform;
+
+  // Perhaps we should allow metersPerUnit to take any string
+  const crsUnit = crs.units as
+    | "m"
+    | "metre"
+    | "meter"
+    | "meters"
+    | "foot"
+    | "us survey foot"
+    | "degree"
+    | undefined;
+
+  if (!crsUnit) {
+    throw new Error(`CRS definition must include "units" property`);
+  }
+
+  const semiMajorAxis = crs.a || crs.datum?.a;
+  const mpu = metersPerUnit(crsUnit, { semiMajorAxis });
+  const cornerOfOrigin: "bottomLeft" | "topLeft" =
+    affine.e(tr) > 0 ? "bottomLeft" : "topLeft";
+
+  const tileMatrices: TileMatrix[] = [];
+
+  // Overviews are sorted finest-to-coarsest; reverse to emit coarsest first.
+  const overviewsCoarseFirst = [...geotiff.overviews].reverse();
+
+  for (let idx = 0; idx < overviewsCoarseFirst.length; idx++) {
+    const overview = overviewsCoarseFirst[idx]!;
+    tileMatrices.push(
+      buildTileMatrix(
+        String(idx),
+        overview.transform,
+        mpu,
+        cornerOfOrigin,
+        overview.tileWidth,
+        overview.tileHeight,
+        overview.width,
+        overview.height,
+      ),
+    );
+  }
+
+  // Full-resolution level is appended last.
+  if (!geotiff.isTiled) {
+    throw new Error("GeoTIFF must be tiled to generate a TMS.");
+  }
+
+  tileMatrices.push(
+    buildTileMatrix(
+      String(geotiff.overviews.length),
+      tr,
+      mpu,
+      cornerOfOrigin,
+      geotiff.tileWidth,
+      geotiff.tileHeight,
+      geotiff.width,
+      geotiff.height,
+    ),
+  );
+
+  const tmsCrs = buildCrs(geotiff.crs);
+  const boundingBox: BoundingBox = {
+    lowerLeft: [bbox[0], bbox[1]],
+    upperRight: [bbox[2], bbox[3]],
+    crs: tmsCrs,
+  };
+
+  return {
+    title: "Generated TMS",
+    id,
+    crs: tmsCrs,
+    boundingBox,
+    tileMatrices,
+  };
+}

packages/geotiff/tests/crs.test.ts

@@ -3,16 +3,13 @@ import wktParser from "wkt-parser";
 import { loadGeoTIFF } from "./helpers.js";
 
 describe("test CRS", () => {
-  it("can fetch EPSG CRS from epsg.io", async () => {
+  it("returns EPSG code", async () => {
     const geotiff = await loadGeoTIFF(
       "uint8_rgb_deflate_block64_cog",
       "rasterio",
     );
     const crs = geotiff.crs;
     expect(crs).toEqual(4326);
-
-    const epsg = geotiff.epsg;
-    expect(epsg).toBe(4326);
   });
 });
 

packages/geotiff/tests/geotiff.test.ts

@@ -178,19 +178,6 @@ describe("GeoTIFF", () => {
     expect(geo.nodata).toBe(-9999);
   });
 
-  it("exposes epsg", async () => {
-    const primary = mockImage({
-      width: 100,
-      height: 100,
-      origin: [0, 0, 0],
-      resolution: [1, -1, 0],
-      epsg: 4326,
-    });
-    const tiff = mockTiff([primary]);
-    const geo = await GeoTIFF.fromTiff(tiff);
-    expect(geo.epsg).toBe(4326);
-  });
-
   it("exposes bbox", async () => {
     const primary = mockImage({
       width: 100,

packages/geotiff/tests/tile-matrix-set.test.ts

@@ -0,0 +1,128 @@
+import { describe, expect, it } from "vitest";
+import wktParser from "wkt-parser";
+import { generateTileMatrixSet } from "../src/tile-matrix-set.js";
+import { loadGeoTIFF } from "./helpers.js";
+
+const EPSG_4326 = {
+  $schema: "https://proj.org/schemas/v0.7/projjson.schema.json",
+  type: "GeographicCRS",
+  name: "WGS 84",
+  datum_ensemble: {
+    name: "World Geodetic System 1984 ensemble",
+    members: [
+      {
+        name: "World Geodetic System 1984 (Transit)",
+        id: { authority: "EPSG", code: 1166 },
+      },
+      {
+        name: "World Geodetic System 1984 (G730)",
+        id: { authority: "EPSG", code: 1152 },
+      },
+      {
+        name: "World Geodetic System 1984 (G873)",
+        id: { authority: "EPSG", code: 1153 },
+      },
+      {
+        name: "World Geodetic System 1984 (G1150)",
+        id: { authority: "EPSG", code: 1154 },
+      },
+      {
+        name: "World Geodetic System 1984 (G1674)",
+        id: { authority: "EPSG", code: 1155 },
+      },
+      {
+        name: "World Geodetic System 1984 (G1762)",
+        id: { authority: "EPSG", code: 1156 },
+      },
+      {
+        name: "World Geodetic System 1984 (G2139)",
+        id: { authority: "EPSG", code: 1309 },
+      },
+      {
+        name: "World Geodetic System 1984 (G2296)",
+        id: { authority: "EPSG", code: 1383 },
+      },
+    ],
+    ellipsoid: {
+      name: "WGS 84",
+      semi_major_axis: 6378137,
+      inverse_flattening: 298.257223563,
+    },
+    accuracy: "2.0",
+    id: { authority: "EPSG", code: 6326 },
+  },
+  coordinate_system: {
+    subtype: "ellipsoidal",
+    axis: [
+      {
+        name: "Geodetic latitude",
+        abbreviation: "Lat",
+        direction: "north",
+        unit: "degree",
+      },
+      {
+        name: "Geodetic longitude",
+        abbreviation: "Lon",
+        direction: "east",
+        unit: "degree",
+      },
+    ],
+  },
+  scope: "Horizontal component of 3D system.",
+  area: "World.",
+  bbox: {
+    south_latitude: -90,
+    west_longitude: -180,
+    north_latitude: 90,
+    east_longitude: 180,
+  },
+  id: { authority: "EPSG", code: 4326 },
+};
+
+describe("test TMS", () => {
+  it("can generate TMS from EPSG CRS", async () => {
+    const geotiff = await loadGeoTIFF(
+      "uint8_rgb_deflate_block64_cog",
+      "rasterio",
+    );
+    const crs = geotiff.crs;
+    expect(crs).toEqual(4326);
+
+    const parsedCrs = wktParser(EPSG_4326);
+
+    const tms = generateTileMatrixSet(geotiff, parsedCrs, { id: "test-tms" });
+
+    expect(tms.crs).toEqual({
+      uri: "http://www.opengis.net/def/crs/EPSG/0/4326",
+    });
+    expect(tms.boundingBox).toEqual({
+      lowerLeft: [0.0, -1.28],
+      upperRight: [1.28, 0.0],
+      crs: { uri: "http://www.opengis.net/def/crs/EPSG/0/4326" },
+    });
+    expect(tms.tileMatrices).toEqual([
+      {
+        id: "0",
+        scaleDenominator: 7951392.199519542,
+        cellSize: 0.02,
+        cornerOfOrigin: "topLeft",
+        pointOfOrigin: [0.0, 0.0],
+        tileWidth: 64,
+        tileHeight: 64,
+        matrixWidth: 1,
+        matrixHeight: 1,
+      },
+      {
+        id: "1",
+        scaleDenominator: 3975696.099759771,
+        cellSize: 0.01,
+        cornerOfOrigin: "topLeft",
+        pointOfOrigin: [0.0, 0.0],
+        tileWidth: 64,
+        tileHeight: 64,
+        matrixWidth: 2,
+        matrixHeight: 2,
+      },
+    ]);
+  });
+});