Handle rotated, skewed or flipped GeoTIFF tile grids

[mike-000]

Fixes #15294
Fixes #15370
Fixes #15428

Uses rotated/skewed/flipped variants of projections to handle GeoTIFF tile grids with rotation, skew or flip specified by ModelTransformation.

See https://deploy-preview-15402--ol-site.netlify.app/en/latest/examples/cog-modeltransformation.html

Uses a local alternative to geotiffjs/geotiff.js#403 which could be removed if the replacement geotiffjs/geotiff.js#413 (with corrected Y axis scale) is merged. The code could be further simplified if geotiffjs/geotiff.js#414 is merged.

[github-actions]

📦 Preview the website for this branch here: https://deploy-preview-15402--ol-site.netlify.app/.

[mike-000]

Ready for review as I do not have any files small enough to use for a rotated or skewed rendering test, although the test file from #15428 would be suitable. The Projection setMatrix method should work equally well for a rotated ImageStatic, etc. but I have not documented that as part of the API.

[github-actions]

📦 Preview the website for this branch here: https://deploy-preview-15402--ol-site.netlify.app/.

[loganwilliams]

@mike-000 Thanks for this! I'm happy to provide more test files, though the smallest I have is ~80MB.

[mike-000]

@loganwilliams There is a large selection of open data available via https://radiantearth.github.io/stac-browser/#/external/s3.us-west-2.amazonaws.com/umbra-open-data-catalog/stac/catalog.json typically around that size suitable for use in examples but too big to use as local test data. I have now found a reduced version of the file used in the example at https://github.com/GeoTIFF/test-data which is under 4MB so might be suitable for another rendering test.

[loganwilliams]

This is really great, thanks for your work on it @mike-000. Hope it gets merged soon.

[ahocevar]

The approach developed in this pull request is good. My suggestions are only about where parts of the implementation belong to, so more users can benefit from it.

[ahocevar]

Wouldn't it make more sense to contribute this to geotiff.js?

[mike-000]

I agree, if geotiffjs/geotiff.js#413 and geotiffjs/geotiff.js#414 are accepted and included in a release then these copies could be replaced by the real thing.

[ahocevar]

Both are now merged and released.

[ahocevar]

I wouldn't bake matrix transforms into projections. It does not feel right. Instead, a matrix transform should be defined with custom transform functions using addCoordinateTransforms(), like e.g. shown in https://openlayers.org/en/latest/examples/wms-custom-proj.html.

[mike-000]

The problem with that approach is transforms would need to be maintained between every matrix variation and every projection added (past and future) via proj4. Given that this is a relatively unusual situation (although it could also be applied to KML ground overlays) I think setting an additional property on the projection and calculating the transforms dynamically as required would be easier.

[ahocevar]

When you create a transform function, you can mix in the matrix transform from a separate function or variable, so you can still handle matrix variations outside the once created transform function. Basically the same thing you're currently doing in src/ol/proj.js.

[mike-000]

The existing transform functions are indexed by projection codes - the matrixes are assigned to projection objects with those codes. Doing it this way avoids adding complexity to the indexing of the primary transforms.

[ahocevar]

Can you point me to the code section where that happens?

[mike-000]

The new code calls projection.setMatrix(matrix); or this.projection.setMatrix(matrix); when determining the projection https://github.com/mike-000/openlayers/blob/7cc890b34747feb99173b974d2e6e143c88e0504/src/ol/source/GeoTIFF.js#L297
and if a matrix is detected when an override projection has been specified
https://github.com/mike-000/openlayers/blob/7cc890b34747feb99173b974d2e6e143c88e0504/src/ol/source/GeoTIFF.js#615

It is handled in existing reprojection code via getTransform and transform calls which need no changes as they are passed the projection objects (not codes) and therefore pass them on unchanged to getTransformFromProjections.

[mike-000]

@ahocevar This is ready for further review. Workarounds have been removed in favour of the latest geotiff.js methods. The handling of matrixes in projections now uses entire transform functions (in series if necessary) instead of very inefficiently extracting a single coordinate transform from a projection transform function to use in a new transform function. This has made the rendering noticeably faster. In many cases (as in the two examples and current tests) the matrix transform is independent of projection transforms and could potentially be cached and/or used directly, e.g.

  if ((!sourceTransform || !destinationTransform) && !transformFunc) {
    return sourceTransform || destinationTransform;
  }

although I suspect there would be little further performance gain from this.

Matrix and projection transforms used in series will be needed for generic viewers such as attempted in https://radiantearth.github.io/stac-browser/#/external/s3.us-west-2.amazonaws.com/umbra-open-data-catalog/stac/2023/2023-02/2023-02-12/6d584f33-0489-47dd-9412-14d2c83532fc/6d584f33-0489-47dd-9412-14d2c83532fc.json where the rotated UTM is further transformed to 3857 as in https://jsfiddle.net/w4qxs50k/ - a test is still needed for that (potentially using an ImageStatic as in #7255 (comment)).

[mike-000]

A single transform could be used directly here, but any performance gain from the extra code would be minimal

  if ((!sourceTransform || !destinationTransform) && !transformFunc) {
    return sourceTransform || destinationTransform;
  }

[ahocevar]

Thanks for your continued effort here, @mike-000. The one thing that still does not look right is that you're using Projection as container for the transform matrix, and the projection transform function to apply a transform matrix. Like I said already, that does not feel like the right place.

Tthinking about it more, my previous suggestion to use transform functions directly does not feel right.

The GeoTIFF source looks like the right place to maintain a mapping image -> matrix, and to make use of it in an overridden getTile() method.

[mike-000]

Separating matrix and projection transforms has the potential for double reprojection. I would regard the image as being in a rotated version of a standard projection - in the example case a UTM projection - defined by a standard projection code and a rotation matrix, so to render it in for example EPSG:3857 only one reprojection should be needed (no reprojection of reproj tiles) the transform function being a combination of the UTM to 3857 transform and the matrix transform which can be easily done by getTransformFromProjections using only its current two parameters (and if source and destination codes are the same only the matrix is used). Having the overall projection defined by a code in the projection and a matrix in the source would be difficult for the current renderer to manage.

Similarly the technique could be used for KML ground overlays, by regarding them as being in rotated/scaled version of EPSG:4326, they could be easily rendered in any view projection which has a transform to EPSG:4326.

But I do think setting a matrix on a projection object should not be exposed in the API, it should be done only inside the GeoTIFF source code (and potentially inside a ground overlay parser) to avoid any future breaking changes.

[ahocevar]

I agree that a single transform of raster data is preferable. But still, this can be done in ol/source/DataTile#getReprojTile_(), if the ReporjDataTile is configured with a transform and a projection.

[mike-000]

A solution which is not specific to DataTile would be preferable, so it could be easily reused elsewhere, e.g. as part of a fix for #4949/#7255 where rotation transforms can be combined projection transforms without the need to maintain large numbers of individual transforms for multiple projection code/rotation permutations.

Also reprojection inevitably reduces the output quality, sometimes as in #15579 that can be unacceptable, so it would nice to allow the view to be in the tiles native projection and reproject any less important background or vector data. Currently getView() returns a standard projection familiar to the user, and the fact that the source may have a non-standard projection with the same code is not documented, but for best quality that can be used as the view projection as in https://jsfiddle.net/jo9wav6y/ something I think should be retained and documented.

[ahocevar]

@mike-000 To make it easier for me to review, please answer one question: in Projection, is it now either matrix or transform function, or both? If both, then how does that work exactly?

[mike-000]

@mike-000 To make it easier for me to review, please answer one question: in Projection, is it now either matrix or transform function, or both? If both, then how does that work exactly?

The projection code together with an optional transform matrix defines a projection. getTransformFromProjections considers the proj4 transform between the codes together with each of the transform matrixes (if present) to return a complete transform between any pair of projections.

[ahocevar]

@mike-000 Thanks for the explanation. I wish it was the other way, i.e. either numeric transform or matrix transform. Now that it's not like that, what exactly does

getTransformFromProjections considers the proj4 transform between the codes together with each of the transform matrixes (if present) to return a complete transform between any pair of projections.

mean?

I stell feel that in this case, registering transforms with transform functions (without ol/proj/proj4) that do whatever you do now to achieve the rotation, skew or flip or whatever.

[mike-000]

what exactly does

getTransformFromProjections considers the proj4 transform between the codes together with each of the transform matrixes (if present) to return a complete transform between any pair of projections.

mean?

The process necessary to display a rotated or flipped UTM projection tile grid as in the example or in #15923 in a generic viewer such as http:https://app.geotiff.io/ which uses a global projection (or in some cases the opposite as in https://jsfiddle.net/jo9wav6y/ to reproject a background such as OSM to the camera projection). When the source is loaded we know its projection code and ModelTransformation, but we do not know the view projection and may need look up the source projection proj4 definition to obtain that part of the transform. The adjustments for ModelTransformation do not need a lookup and can easily be added on demand.

[ahocevar]

Well, like you're saying - the whole transform needed here is a "process", not a matter of projection. If you manage to change your pull request to do the transform with a chain of transform functions (i.e. a process) instead of trying to squeeze this into Projection, I'll review again.

[mike-000]

So how do you propose to do that? Since all methods which take a source projection and destination projection eventually call getTransformFromProjections that seems the obvious place to do it. If the matrix is not included in the projection then each method which takes source and destination projection would also need source and destination matrix options.

[mike-000]

Would you prefer to create unique projections with codes such as

EPSG:32759?[-0.24064389665799932,0.9706134735317022,-0.9706134735316878,-0.24064389665796512,8034048.497358791,9399383.005560089]

(where we would still need to wait until transforms for EPSG:32759 were added before creating transforms for the above)?

[ahocevar]

I can answer that once I fully understand what you're aiming for, which is not the case. So one more question: why are you saying the transform matrix is applied after reprojecting? I would think that if a GeoTIFF is rotated, skewed or flipped, that applying a transform matrix would be the first thing to do. Then we're done, unless the result is still in a different projection than the View. If so, a source -> view projection transform will be applied.

[mike-000]

Yes, that is the problem situation. Consider an OpenLayers version of https://app.geotiff.io/ where the code might include

register(proj4);
const map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new OSM(),
    }),
  ],
  view: new View({
    center: [0, 0],
    zoom: 0,
  }),
});

// call this when URL is entered in form

function addGeoTIFF(url) {
  const cogSource = new GeoTIFF({
    sources: [{url: url}],
  });

  cogSource
    .getView()
    .then((viewConfig) =>
      fromEPSGCode(viewConfig.projection.getCode()).then(() => {
        map.getView().fit(
          transformExtent(
            viewConfig.extent,
            viewConfig.projection,
            map.getView().getProjection(),
          ),
        );
        map.addlayer(
          new TileLayer({
            source: cogSource,
          }),
        );
      })
    );
}

To avoid reprojecting a reprojection, which would affect output quality, the tiles should be reprojected directly from the rotated UTM projection to the view projection - the overall transform for doing that is obtained by combining the matrix and proj4 transforms (only the proj4 transform is cached, transform functions for a matrix are easy to create as required).

That is more complex than in the example where only the matrix transform is used to reproject the tiles (to an unrotated reprojection tile grid) and the proj4 transform is only needed to reproject the background OSM to UTM.

[ahocevar]

Ok, so the matrix transform comes first, then the projection transform. So once you have a transform function that combines both, you don't have to create a transformed canvas from a transformed canvas. Try to think less in what we currently have (i.e. Projection), think more in what would make most sense.

For example, a renderer could get data from a source, obtain a transform matrix from the source, get the projection from the view, and then do what we currently call "reprojection" using a function that combines the matrix transform and the reprojection transform functions.

[mike-000]

I would consider the projection with matrix to be a projection in its own right as a tile grid must be aligned with the source projection axes, so it does make sense to me. And as noted above

If the matrix is not included in the projection then each method which takes source and destination projection would also need source and destination matrix options.

so any other method is likely to be more complex to implement (especially for something which will be relatively uncommon) and if done at renderer level would not be reusable for other layer types (e.g. a KML GroundOverlay could be easily implemented as an ImageStatic in a scaled and rotated EPSG:4326 based projection and reprojected into any view projection with minimal additional effect).

Also it would be desirable to be able use the rotated source projection directly at times as in https://jsfiddle.net/jo9wav6y/ as any reprojection will affect output quality and in a rotated version of a COG such as in #15579 changing the pixel alignment would make the output even more meaningless than the unwanted resizing of pixels.

[ahocevar]

Again, I'm not talking about multiple transforms on a canvas. Just chaining of transform functions. And my suggestion is to separate affine transforms (source matrix transform) from projections (source to view projection). An API for the user could look like this:

new Source({
  projection: sourceProjection, 
  transform: transformMatrix, // transform data into source projection
});

For GeoTIFF, both can of course be determined from the GeoTIFF metadata.

[mike-000]

Why make it unnecessarily complicated for the rendererers? getTransformFromProjections is already chaining transforms, and it does that using the two existing parameters which can now contain both the code for the projection and the transform matrix.

[ahocevar]

We should aim for making it easy for the users. I don't see where this would make it more complicated in the renderer code.

The problem with your approach is that projections work based on a projection code lookup. And with an arbitrary matrix in addition to a projection that can uniquely be identified by a code, we get into a situation where we need to set up a code (or a projection object) for each transform+projection combination.

[mike-000]

Yes, we need new projection objects for custom projections - before each call to setMatrix() the projection on which is based needs to be cloned (ol/source/GeoTIFF does that). I do not see any problem with that as long as setMatrix() is not exposed in the API, but maybe a non-API createProjectionWithMatrix method would be safer?

[ahocevar]

I think I have made valid suggestions and also shown a way how these could be implemented. Please let's bury your selected approach of trying to squeeze all that into Projection objects.

[mike-000]

Your proposal works well when the source is reprojected to another projection, e.g. UTM source with a rotation matrix reprojected to 3857. The matrix needs to be passed to ReprojDataTile and on to Triangulation to be included in the transformInv. Ideally it should also be passed to calculateSourceExtentResolution and getPointResolution in case the matrix has significant scaling. But when attempting to use reprojection with the same projection (source with matrix and view without matrix) it starts getting messy. A reprojection must be forced in getTile() if the source has a matrix, but if the same getTile() is called from ReprojDataTile it must return a source tile and not recursively create reproj tiles. I managed that but still cannot get any output. I suspect the need for separate caches, tile grid and gutters for source and reproj tiles, indexed by projection object uid, may be the problem.

[ahocevar]

Maybe some things can be simplified along the way:

• Reprojected tiles in view projection can use a standard tile grid
• Reprojected tiles in view projection do not need gutters
• Tiles in view projection are cached by the renderer since New image tile source #15905, so no need to keep a lookup of caches by projection in the source any more.

It sounds like that would resolve most or your issues.

[mike-000]

Yes, I have it mostly working now, by taking the source transform matrix into account in getTileGridForProjection, getTileCacheForProjection and getGutterForProjection. One of the test samples is not rendering the full extent, which I don't understand as the overall transform should be the same regardless of where the components are obtained from:

[ahocevar]

Yes, I have it mostly working now, by taking the source transform matrix into account in getTileGridForProjection, getTileCacheForProjection and getGutterForProjection.

Interesting. Would you mind sharing your code? I'm unable to come up with a reason why the transform would matter for the view projection.