feat: Extract TileMatrixSet definition from GeoTIFF object

packages/deck.gl-cog/src/cog-tile-matrix-set.ts

@@ -0,0 +1,216 @@
+import type {
+  TileMatrix,
+  TileMatrixSet,
+} from "@developmentseed/deck.gl-raster";
+import type { GeoTIFF, GeoTIFFImage } from "geotiff";
+import {
+  extractGeotransform,
+  getGeoTIFFProjection,
+} from "./geotiff-reprojection";
+import proj4, { ProjectionDefinition } from "proj4";
+import Ellipsoid from "./ellipsoids.js";
+import type { PROJJSONDefinition } from "proj4/dist/lib/core";
+
+// 0.28 mm per pixel
+// https://docs.ogc.org/is/17-083r4/17-083r4.html#toc15
+const SCREEN_PIXEL_SIZE = 0.00028;
+
+/**
+ *
+ * Ported from Vincent's work here:
+ * https://github.com/developmentseed/morecantile/pull/187/changes#diff-402eedddfa30af554d03750c8a82a09962b44b044976c321b774b484b98e8f48R665
+ *
+ * @return  {TileMatrixSet}[return description]
+ */
+export async function parseCOGTileMatrixSet(
+  tiff: GeoTIFF,
+): Promise<TileMatrixSet> {
+  const fullResImage = await tiff.getImage();
+
+  if (!fullResImage.isTiled) {
+    throw new Error("COG TileMatrixSet requires a tiled GeoTIFF");
+  }
+
+  const imageCount = await tiff.getImageCount();
+  const bbox = fullResImage.getBoundingBox();
+  const fullImageWidth = fullResImage.getWidth();
+  const fullImageHeight = fullResImage.getHeight();
+
+  const crs = await getGeoTIFFProjection(fullResImage);
+
+  if (crs === null) {
+    throw new Error(
+      "Could not determine coordinate reference system from GeoTIFF geo keys",
+    );
+  }
+
+  const parsedCrs = parseCrs(crs);
+
+  const projectToWgs84 = proj4(crs, "EPSG:4326").forward;
+  const projectTo3857 = proj4(crs, "EPSG:3857").forward;
+
+  const boundingBox: TileMatrixSet["boundingBox"] = {
+    lowerLeft: [bbox[0]!, bbox[1]!],
+    upperRight: [bbox[2]!, bbox[3]!],
+  };
+
+  const transform = extractGeotransform(fullResImage);
+
+  if (transform[1] !== 0 || transform[3] !== 0) {
+    // TileMatrixSet assumes orthogonal axes
+    throw new Error(
+      "COG TileMatrixSet with rotation/skewed geotransform is not supported",
+    );
+  }
+
+  const cellSize = Math.abs(transform[0]);
+
+  const tileWidth = fullResImage.getTileWidth();
+  const tileHeight = fullResImage.getTileHeight();
+
+  const tileMatrices: TileMatrix[] = [
+    {
+      // Set as highest resolution / finest level
+      id: String(imageCount - 1),
+      scaleDenominator:
+        (cellSize * metersPerUnit(parsedCrs)) / SCREEN_PIXEL_SIZE,
+      cellSize,
+      pointOfOrigin: [transform[2], transform[5]],
+      tileWidth: fullResImage.getTileWidth(),
+      tileHeight: fullResImage.getTileHeight(),
+      matrixWidth: Math.ceil(fullImageWidth / tileWidth),
+      matrixHeight: Math.ceil(fullImageHeight / tileHeight),
+      geotransform: transform,
+    },
+  ];
+
+  // Starting from 1 to skip full res image
+  for (let imageIdx = 1; imageIdx < imageCount; imageIdx++) {
+    const image = await tiff.getImage(imageIdx);
+
+    if (!image.isTiled) {
+      throw new Error("COG TileMatrixSet requires a tiled GeoTIFF");
+    }
+
+    const tileMatrix = createOverviewTileMatrix({
+      id: String(imageCount - 1 - imageIdx),
+      image,
+      fullWidth: fullImageWidth,
+      fullHeight: fullImageHeight,
+      baseTransform: transform,
+      parsedCrs,
+    });
+    tileMatrices.push(tileMatrix);
+  }
+
+  // Reverse to have coarsest level first
+  tileMatrices.reverse();
+
+  return {
+    crs,
+    boundingBox,
+    tileMatrices,
+    projectToWgs84,
+    projectTo3857,
+  };
+}
+
+/**
+ * Coefficient to convert the coordinate reference system (CRS)
+ * units into meters (metersPerUnit).
+ *
+ * From note g in http://docs.opengeospatial.org/is/17-083r2/17-083r2.html#table_2:
+ *
+ * > If the CRS uses meters as units of measure for the horizontal dimensions,
+ * > then metersPerUnit=1; if it has degrees, then metersPerUnit=2pa/360
+ * > (a is the Earth maximum radius of the ellipsoid).
+ */
+// https://github.com/developmentseed/morecantile/blob/7c95a11c491303700d6e33e9c1607f2719584dec/morecantile/utils.py#L67-L90
+function metersPerUnit(parsedCrs: ProjectionDefinition): number {
+  switch (parsedCrs.units) {
+    case "metre":
+    case "meter":
+    case "meters":
+      return 1;
+    case "foot":
+      return 0.3048;
+    case "US survey foot":
+      return 1200 / 3937;
+  }
+
+  if (parsedCrs.units === "degree") {
+    // 2 * π * ellipsoid semi-major-axis / 360
+    const { a } = Ellipsoid[parsedCrs.ellps as keyof typeof Ellipsoid];
+    return (2 * Math.PI * a) / 360;
+  }
+
+  throw new Error(`Unsupported CRS units: ${parsedCrs.units}`);
+}
+
+/**
+ * Create tile matrix for COG overview
+ */
+function createOverviewTileMatrix({
+  id,
+  image,
+  fullWidth,
+  baseTransform,
+  parsedCrs,
+}: {
+  id: string;
+  image: GeoTIFFImage;
+  fullWidth: number;
+  fullHeight: number;
+  baseTransform: [number, number, number, number, number, number];
+  parsedCrs: ProjectionDefinition;
+}): TileMatrix {
+  const width = image.getWidth();
+  const height = image.getHeight();
+
+  // For now, just use scaleX
+  // https://github.com/developmentseed/morecantile/pull/187/changes#r2621314673
+  const scaleX = fullWidth / width;
+
+  // const scaleY = fullHeight / height;
+  // if (Math.abs(scaleX - scaleY) > 1e-3) {
+  //   throw new Error("Non-uniform overview scaling detected (X/Y differ)");
+  // }
+
+  const scale = scaleX;
+
+  const geotransform: [number, number, number, number, number, number] = [
+    baseTransform[0] * scale,
+    baseTransform[1] * scale,
+    baseTransform[2], // x origin stays the same
+    baseTransform[3] * scale,
+    baseTransform[4] * scale,
+    baseTransform[5], // y origin stays the same
+  ];
+  const cellSize = Math.abs(geotransform[0]);
+
+  const tileWidth = image.getTileWidth();
+  const tileHeight = image.getTileHeight();
+
+  return {
+    id,
+    scaleDenominator: (cellSize * metersPerUnit(parsedCrs)) / SCREEN_PIXEL_SIZE,
+    cellSize,
+    pointOfOrigin: [geotransform[2], geotransform[5]],
+    tileWidth,
+    tileHeight,
+    matrixWidth: Math.ceil(width / tileWidth),
+    matrixHeight: Math.ceil(height / tileHeight),
+    geotransform,
+  };
+}
+
+function parseCrs(crs: PROJJSONDefinition): ProjectionDefinition {
+  // If you pass proj4.defs a projjson, it doesn't parse it; it just returns the
+  // input.
+  //
+  // Instead, you need to assign it to an alias and then retrieve it.
+
+  const key = "__deck.gl-cog-internal__";
+  proj4.defs(key, crs);
+  return proj4.defs(key);
+}

packages/deck.gl-cog/src/ellipsoids.ts

@@ -0,0 +1,154 @@
+/**
+ * Vendored and edited from proj4.js.
+ *
+ * In the implementation of metersPerUnit while generated a COG TileMatrixSet,
+ * we need to know the size of the semi major axis in the case the CRS is in
+ * degrees.
+ *
+ * https://github.com/proj4js/proj4js/blob/e90e5fa6872a1ffc40edb161cbeb4bd5e3bd9db5/lib/constants/Ellipsoid.js
+ */
+
+const ellipsoids = {
+  MERIT: {
+    a: 6378137,
+  },
+  SGS85: {
+    a: 6378136,
+  },
+  GRS80: {
+    a: 6378137,
+  },
+  IAU76: {
+    a: 6378140,
+  },
+  airy: {
+    a: 6377563.396,
+    b: 6356256.91,
+  },
+  APL4: {
+    a: 6378137,
+  },
+  NWL9D: {
+    a: 6378145,
+  },
+  mod_airy: {
+    a: 6377340.189,
+    b: 6356034.446,
+  },
+  andrae: {
+    a: 6377104.43,
+  },
+  aust_SA: {
+    a: 6378160,
+  },
+  GRS67: {
+    a: 6378160,
+  },
+  bessel: {
+    a: 6377397.155,
+  },
+  bess_nam: {
+    a: 6377483.865,
+  },
+  clrk66: {
+    a: 6378206.4,
+    b: 6356583.8,
+  },
+  clrk80: {
+    a: 6378249.145,
+  },
+  clrk80ign: {
+    a: 6378249.2,
+    b: 6356515,
+  },
+  clrk58: {
+    a: 6378293.645208759,
+  },
+  CPM: {
+    a: 6375738.7,
+  },
+  delmbr: {
+    a: 6376428,
+  },
+  engelis: {
+    a: 6378136.05,
+  },
+  evrst30: {
+    a: 6377276.345,
+  },
+  evrst48: {
+    a: 6377304.063,
+  },
+  evrst56: {
+    a: 6377301.243,
+  },
+  evrst69: {
+    a: 6377295.664,
+  },
+  evrstSS: {
+    a: 6377298.556,
+  },
+  fschr60: {
+    a: 6378166,
+  },
+  fschr60m: {
+    a: 6378155,
+  },
+  fschr68: {
+    a: 6378150,
+  },
+  helmert: {
+    a: 6378200,
+  },
+  hough: {
+    a: 6378270,
+  },
+  intl: {
+    a: 6378388,
+  },
+  kaula: {
+    a: 6378163,
+  },
+  lerch: {
+    a: 6378139,
+  },
+  mprts: {
+    a: 6397300,
+  },
+  new_intl: {
+    a: 6378157.5,
+    b: 6356772.2,
+  },
+  plessis: {
+    a: 6376523,
+  },
+  krass: {
+    a: 6378245,
+  },
+  SEasia: {
+    a: 6378155,
+    b: 6356773.3205,
+  },
+  walbeck: {
+    a: 6376896,
+    b: 6355834.8467,
+  },
+  WGS60: {
+    a: 6378165,
+  },
+  WGS66: {
+    a: 6378145,
+  },
+  WGS7: {
+    a: 6378135,
+  },
+  WGS84: {
+    a: 6378137,
+  },
+  sphere: {
+    a: 6370997,
+    b: 6370997,
+  },
+};
+
+export default ellipsoids;