Handle rotated, skewed or flipped GeoTIFF tile grids

examples/cog-modeltransformation.html

@@ -0,0 +1,11 @@
+---
+layout: example.html
+title: COG with ModelTransformation
+shortdesc: COG with ModelTransformation.
+docs: >
+  Example of using a COG with ModelTransformation.
+  A ModelTransformation may rotate or skew the tile grid relative to the reported projection.
+  This can be reprojected to that or any other supported projection.
+tags: "cog, projection, ModelTransformation"
+---
+<div id="map" class="map"></div>

examples/cog-modeltransformation.js

@@ -0,0 +1,34 @@
+import GeoTIFF from '../src/ol/source/GeoTIFF.js';
+import Map from '../src/ol/Map.js';
+import OSM from '../src/ol/source/OSM.js';
+import TileLayer from '../src/ol/layer/WebGLTile.js';
+import proj4 from 'proj4';
+import {fromEPSGCode, register} from '../src/ol/proj/proj4.js';
+
+register(proj4);
+
+const cogSource = new GeoTIFF({
+  sources: [
+    {
+      url: 'https://umbra-open-data-catalog.s3.amazonaws.com/sar-data/tasks/Tanna%20Island,%20Vanuatu/9c76a918-9247-42bf-b9f6-3b4f672bc148/2023-02-12-21-33-56_UMBRA-04/2023-02-12-21-33-56_UMBRA-04_GEC.tif',
+    },
+  ],
+});
+
+const map = new Map({
+  target: 'map',
+  layers: [
+    new TileLayer({
+      source: new OSM(),
+    }),
+    new TileLayer({
+      source: cogSource,
+      opacity: 0.8,
+    }),
+  ],
+  view: cogSource
+    .getView()
+    .then((viewConfig) =>
+      fromEPSGCode(viewConfig.projection.getCode()).then(() => viewConfig),
+    ),
+});

package.json

@@ -48,7 +48,7 @@
     "color-rgba": "^3.0.0",
     "color-space": "^2.0.1",
     "earcut": "^2.2.3",
-    "geotiff": "^2.0.7",
+    "geotiff": "^2.1.3",
     "pbf": "3.2.1",
     "rbush": "^3.0.1"
   },

src/ol/proj.js

@@ -71,9 +71,11 @@ import {
   clear as clearTransformFuncs,
   get as getTransformFunc,
 } from './proj/transforms.js';
+import {apply as applyMatrix, invert as invertMatrix} from './transform.js';
 import {applyTransform, getWidth} from './extent.js';
 import {clamp, modulo} from './math.js';
 import {equals, getWorldsAway} from './coordinate.js';
+import {equals as equalsMatrix} from './array.js';
 import {getDistance} from './sphere.js';
 import {warn} from './console.js';
 
@@ -440,6 +442,14 @@ export function equivalent(projection1, projection2) {
   if (projection1 === projection2) {
     return true;
   }
+  const matrix1 = projection1.getMatrix();
+  const matrix2 = projection2.getMatrix();
+  if (
+    (matrix1 || matrix2) &&
+    (!matrix1 || !matrix2 || !equalsMatrix(matrix1, matrix2))
+  ) {
+    return false;
+  }
   const equalUnits = projection1.getUnits() === projection2.getUnits();
   if (projection1.getCode() === projection2.getCode()) {
     return equalUnits;
@@ -463,11 +473,55 @@ export function getTransformFromProjections(
 ) {
   const sourceCode = sourceProjection.getCode();
   const destinationCode = destinationProjection.getCode();
-  let transformFunc = getTransformFunc(sourceCode, destinationCode);
-  if (!transformFunc) {
-    transformFunc = identityTransform;
+  const transformFunc = getTransformFunc(sourceCode, destinationCode);
+  const sourceMatrix = sourceProjection.getMatrix();
+  const destinationMatrix = destinationProjection.getMatrix();
+  if (!sourceMatrix && !destinationMatrix) {
+    return transformFunc || identityTransform;
+  }
+
+  let sourceTransform, destinationTransform;
+  if (sourceMatrix) {
+    const matrix = invertMatrix(sourceMatrix.slice());
+    sourceTransform = createTransformFromCoordinateTransform(
+      function (coordinate) {
+        return applyMatrix(matrix, coordinate);
+      },
+    );
   }
-  return transformFunc;
+  if (destinationMatrix) {
+    const matrix = destinationMatrix;
+    destinationTransform = createTransformFromCoordinateTransform(
+      function (coordinate) {
+        return applyMatrix(matrix, coordinate);
+      },
+    );
+  }
+  return (
+    /**
+     * @param {Array<number>} input Input.
+     * @param {Array<number>} [output] Output.
+     * @param {number} [dimension] Dimension.
+     * @return {Array<number>} Output.
+     */
+    function (input, output, dimension) {
+      const length = input.length;
+      output = output !== undefined ? output : new Array(length);
+      for (let i = 0; i < length; ++i) {
+        output[i] = input[i];
+      }
+      if (sourceTransform) {
+        sourceTransform(output, output, dimension);
+      }
+      if (transformFunc) {
+        transformFunc(output, output, dimension);
+      }
+      if (destinationTransform) {
+        destinationTransform(output, output, dimension);
+      }
+      return output;
+    }
+  );
 }
 
 /**

src/ol/proj/Projection.js

@@ -124,6 +124,13 @@ class Projection {
      * @type {number|undefined}
      */
     this.metersPerUnit_ = options.metersPerUnit;
+
+    /**
+     * Transform matrix to rotate or skew the projection.
+     * @private
+     * @type {import("../transform.js").Transform|undefined}
+     */
+    this.matrix_ = undefined;
   }
 
   /**
@@ -266,6 +273,20 @@ class Projection {
   getPointResolutionFunc() {
     return this.getPointResolutionFunc_;
   }
+
+  /**
+   * @return {import("../transform.js").Transform|undefined} Transform matrix.
+   */
+  getMatrix() {
+    return this.matrix_;
+  }
+
+  /**
+   * @param {import("../transform.js").Transform|undefined} matrix Transform matrix.
+   */
+  setMatrix(matrix) {
+    this.matrix_ = matrix;
+  }
 }
 
 export default Projection;

src/ol/source/GeoTIFF.js

@@ -15,10 +15,12 @@ import {
   get as getCachedProjection,
   toUserCoordinate,
   toUserExtent,
+  transformExtent,
 } from '../proj.js';
 import {clamp} from '../math.js';
 import {getCenter, getIntersection} from '../extent.js';
 import {error as logError} from '../console.js';
+import {multiply as multiplyTransform} from '../transform.js';
 import {fromCode as unitsFromCode} from '../proj/Units.js';
 
 /**
@@ -131,7 +133,7 @@ function getWorkerPool() {
  */
 function getBoundingBox(image) {
   try {
-    return image.getBoundingBox();
+    return image.getBoundingBox(true);
   } catch (_) {
     return [0, 0, image.getWidth(), image.getHeight()];
   }
@@ -171,43 +173,69 @@ function getResolutions(image, referenceImage) {
 
 /**
  * @param {GeoTIFFImage} image A GeoTIFF.
+ * @param {import("../proj/Projection.js").default} [defaultProjection] Default projection.
  * @return {import("../proj/Projection.js").default} The image projection.
  */
-function getProjection(image) {
-  const geoKeys = image.geoKeys;
-  if (!geoKeys) {
-    return null;
-  }
+function getProjection(image, defaultProjection) {
+  let projection = defaultProjection;
+
+  if (!projection) {
+    const geoKeys = image.geoKeys;
+    if (!geoKeys) {
+      return null;
+    }
 
-  if (
-    geoKeys.ProjectedCSTypeGeoKey &&
-    geoKeys.ProjectedCSTypeGeoKey !== 32767
-  ) {
-    const code = 'EPSG:' + geoKeys.ProjectedCSTypeGeoKey;
-    let projection = getCachedProjection(code);
-    if (!projection) {
-      const units = unitsFromCode(geoKeys.ProjLinearUnitsGeoKey);
-      if (units) {
+    let code, units;
+    if (
+      geoKeys.ProjectedCSTypeGeoKey &&
+      geoKeys.ProjectedCSTypeGeoKey !== 32767
+    ) {
+      code = 'EPSG:' + geoKeys.ProjectedCSTypeGeoKey;
+      units = unitsFromCode(geoKeys.ProjLinearUnitsGeoKey);
+    } else if (
+      geoKeys.GeographicTypeGeoKey &&
+      geoKeys.GeographicTypeGeoKey !== 32767
+    ) {
+      code = 'EPSG:' + geoKeys.GeographicTypeGeoKey;
+      units = unitsFromCode(geoKeys.GeogAngularUnitsGeoKey);
+    }
+
+    if (code) {
+      projection = getCachedProjection(code);
+      if (units && !projection) {
         projection = new Projection({
           code: code,
           units: units,
         });
       }
     }
-    return projection;
   }
 
-  if (geoKeys.GeographicTypeGeoKey && geoKeys.GeographicTypeGeoKey !== 32767) {
-    const code = 'EPSG:' + geoKeys.GeographicTypeGeoKey;
-    let projection = getCachedProjection(code);
-    if (!projection) {
-      const units = unitsFromCode(geoKeys.GeogAngularUnitsGeoKey);
-      if (units) {
-        projection = new Projection({
-          code: code,
-          units: units,
-        });
-      }
+  if (projection) {
+    const modelTransformation = image.fileDirectory.ModelTransformation;
+    if (modelTransformation) {
+      // eslint-disable-next-line no-unused-vars
+      const [a, b, c, d, e, f, g, h] = modelTransformation;
+      const matrix = multiplyTransform(
+        multiplyTransform(
+          [
+            1 / Math.sqrt(a * a + e * e),
+            0,
+            0,
+            -1 / Math.sqrt(b * b + f * f),
+            d,
+            h,
+          ],
+          [a, e, b, f, 0, 0],
+        ),
+        [1, 0, 0, 1, -d, -h],
+      );
+
+      projection = new Projection({
+        code: projection.getCode(),
+        units: projection.getUnits(),
+      });
+      projection.setMatrix(matrix);
     }
     return projection;
   }
@@ -514,9 +542,21 @@ class GeoTIFFSource extends DataTile {
     const firstSource = sources[0];
     for (let i = firstSource.length - 1; i >= 0; --i) {
       const image = firstSource[i];
-      const projection = getProjection(image);
+      const projection = getProjection(image, this.projection);
       if (projection) {
-        this.projection = projection;
+        const matrix = projection.getMatrix();
+        if (!this.projection) {
+          this.projection = projection;
+        } else if (matrix) {
+          this.projection = new Projection({
+            code: this.projection.getCode(),
+            units: this.projection.getUnits(),
+          });
+          this.projection.setMatrix(matrix);
+        }
+        if (matrix) {
+          this.addAlpha_ = true;
+        }
         break;
       }
     }
@@ -683,9 +723,7 @@ class GeoTIFFSource extends DataTile {
       }
     }
 
-    if (!this.getProjection()) {
-      this.determineProjection(sources);
-    }
+    this.determineProjection(sources);
 
     this.samplesPerPixel_ = samplesPerPixel;
     this.nodataValues_ = nodataValues;
@@ -753,12 +791,25 @@ class GeoTIFFSource extends DataTile {
       resolutions = [resolutions[0] * 2, resolutions[0], resolutions[0] / 2];
     }
 
+    let viewProjection = this.projection;
+    let viewExtent = extent;
+    if (viewProjection && viewProjection.getMatrix()) {
+      viewProjection = getCachedProjection(this.projection.getCode());
+      if (!viewProjection) {
+        viewProjection = new Projection({
+          code: this.projection.getCode(),
+          units: this.projection.getUnits(),
+        });
+      }
+      viewExtent = transformExtent(viewExtent, this.projection, viewProjection);
+    }
+
     this.viewResolver({
       showFullExtent: true,
-      projection: this.projection,
+      projection: viewProjection,
       resolutions: resolutions,
-      center: toUserCoordinate(getCenter(extent), this.projection),
-      extent: toUserExtent(extent, this.projection),
+      center: toUserCoordinate(getCenter(viewExtent), viewProjection),
+      extent: toUserExtent(viewExtent, viewProjection),
       zoom: zoom,
     });
   }