GeoTiff reprojection from WGS84 to Mercator

[maxschmi]

Describe the bug

Openlayers reprojection of a GeoTiff TileLayer from WGS84 (EPS:4326) to Mercator (EPSG:3857) is not working as expected.
The resulting grid has different cell widths. They somehow togle from wide to thin. So one column is wide and the next one is thin.

See the screenshot:

To Reproduce

Steps to reproduce the behavior:

Go to codesandbox of minimal example or checkout the minimal example repository

Expected behavior

The grid cells should have the same width as seen in this screenshot from QGIS of the same layer and the same map projection:

[mike-000]

You need to prevent excessive overzooming of the original tiles during reprojection by limiting the maxZoom of the reprojected tiles, for example

tif_layer.getSource().setTileGridForProjection('EPSG:3857', createXYZ({maxZoom: 4}))

works in your CodeSandbox.

You could attempt to calculate an appropriate setting based on tif_view.resolutions.

[maxschmi]

Thanks for the quick response. At first the sugested solution looked good, but having a closer look I realized it interpolates between cells to get the pixel value.
This should absolutly not be the case, as I really want to show the grid cells value in my application.

I forked the minimal example and added @mike-000's sugested solution.
Furthermore I added a hover to show the pixels value on mouse hover. The values in the test GeoTiffs are all integers, but as openalyer interpolates the values, they become float numbers.
See the codeSandbox for the updated version.

Unfortunatly I don't really get how setTileGridForProjection and createXYZ work, so I couldn't fix this problem. I could only play around and change the values, but it never happened to work out as excpected.

You could attempt to calculate an appropriate setting based on tif_view.resolutions.

How could I calculate this?

[mike-000]

The unwanted interpolation is only seen in Firefox, see #14099. It can be fixed with a polyfill such as

/*
 * Polyfill to make Gecko browser `CanvasRenderingContext2D.drawImage()` behavior when 
 * downscaling with `imageSmoothingEnabled` set to `false` consistent with other browsers
 * (no smoothing or interpolation of pixel values) and not as updated in
 * https://github.com/mdn/content/blob/main/files/en-us/mozilla/firefox/releases/56/index.md#canvas-and-webgl
 */
(function () {

    const canvas1 = document.createElement('canvas');
    canvas1.width = 2;
    canvas1.height = 2;
    const ctx1 = canvas1.getContext('2d');
    const imageData = ctx1.createImageData(2, 2);
    imageData.data[0] = 255;
    imageData.data[3] = 255;
    imageData.data[7] = 255;
    imageData.data[11] = 255;
    imageData.data[12] = 255;
    imageData.data[15] = 255;
    ctx1.putImageData(imageData, 0, 0);
    const canvas2 = document.createElement('canvas');
    canvas2.width = 1;
    canvas2.height = 1;
    const ctx2 = canvas2.getContext('2d');
    ctx2.imageSmoothingEnabled = false;
    ctx2.drawImage(canvas1, 0, 0, 2, 2, 0, 0, 1, 1);
    const data = ctx2.getImageData(0, 0, 1, 1).data;

    if (data[0] !== 0 && data[0] !== 255) {
      // browser interpolates, polyfill is needed
      const defaultDrawImage = CanvasRenderingContext2D.prototype.drawImage;

      CanvasRenderingContext2D.prototype.drawImage = function (
        image,
        ...params
      ) {
        // Calculate scales in case workaround is needed
        let scaleX = 1;
        let scaleY = 1;

        let sx = 0;
        let sy = 0;
        let sWidth = image.width;
        let sHeight = image.height;
        let dx;
        let dy;
        let dWidth = sWidth;
        let dHeight = sHeight;

        if (!this.imageSmoothingEnabled) {
          const transform = this.getTransform();
          scaleX = Math.sqrt(
            transform.a * transform.a + transform.b * transform.b
          );
          scaleY = Math.sqrt(
            transform.c * transform.c + transform.d * transform.d
          );

          if (params.length === 2) {
            [dx, dy] = params;
          } else if (params.length === 4) {
            [dx, dy, dWidth, dHeight] = params;
          } else if (params.length === 8) {
            [sx, sy, sWidth, sHeight, dx, dy, dWidth, dHeight] = params;
          }

          scaleX *= Math.abs(dWidth / sWidth);
          scaleY *= Math.abs(dHeight / sHeight);
        }

        if (dx !== undefined && dy !== undefined && scaleX < 1 && scaleY < 1) {
          // Image reduction workaround to avoid downscaling by
          // Gecko browsers when interpolation is not required
          let toKeepX = Math.abs(sWidth);
          let toKeepY = Math.abs(sHeight);
          let sWidthFinal = toKeepX;
          let sHeightFinal = toKeepY;
          // Reduce whichever dimension needs less reduction
          // keeping at least one row or column
          if (scaleX > scaleY) {
            toKeepX = Math.max(Math.floor(scaleX * toKeepX), 1);
            sWidthFinal = Math.min(toKeepX, scaleX * sWidthFinal);
          } else {
            toKeepY = Math.max(Math.floor(scaleY * toKeepY), 1);
            sHeightFinal = Math.min(toKeepY, scaleY * sHeightFinal);
          }

          const drawImage = document.createElement('canvas');
          drawImage.width = toKeepX;
          drawImage.height = toKeepY;
          const drawContext = drawImage.getContext('2d');
          drawContext.imageSmoothingEnabled = false;
          defaultDrawImage.call(
            drawContext,
            image,
            sx,
            sy,
            sWidth,
            sHeight,
            0,
            0,
            toKeepX,
            toKeepY
          );
          defaultDrawImage.call(
            this,
            drawImage,
            0,
            0,
            sWidthFinal,
            sHeightFinal,
            dx,
            dy,
            dWidth,
            dHeight
          );

          drawImage.width = 1;
          drawImage.height = 1;
        } else {
          defaultDrawImage.call(this, image, ...params);
        }
      };
    }

})();

[maxschmi]

Thanks a lot this did the deal.

[mike-000]

You could attempt to calculate an appropriate setting based on tif_view.resolutions.

How could I calculate this?

The aim would be to avoid having target (reprojected) tile grid resolutions which would result in source resolutions calculated by https://github.com/openlayers/openlayers/blob/main/src/ol/reproj.js#L104 being smaller than those available in the source.

Maybe this could/should be attempted for all reprojections? However as in https://openlayers.org/en/latest/examples/geotiff-reprojection.html what might work well at the equator may not be suitable for polar regions or vice versa.

[maxschmi]

Hmm even after hours of playing around with createXYZ I couldn't get the right parameters to prevent having columns where 2 grid cells get merged. So I allways end up having some merged columns or overzooming efects as in the beginning.

Here is my actual setup: codesandbox

I tried setting maxResolution, tileSize, minZoom, maxZoom, extent parameters of the createXYZ function, but it feels more like guessing around and maybe being lucky one day. Isn't there a more systematic approach on how to calculate this? Especially as the horizontal difference in the cell width between Mercator and WGS84 in my region is not very large I think that this should somehow be possible.

[mike-000]

In your case your source has non-square-pixels i.e. the rendered aspect ratio of tiles in the tile grid does not match the 256 : 256 ratio of their data. Reprojection only produces tiles with square pixels and scales up one dimension to achieve that.

This code will maintain alignment of columns as in the unreprojected source, but the row alignment is not precise

const source = tif_layer.getSource();
const tileGrid = source.getTileGrid();
const resolutions = tileGrid.getResolutions();
const maxZoom = resolutions.length - 1;
const resolution =
  (resolutions[maxZoom] * getWidth(map_proj.getExtent())) /
  getWidth(tif_view.projection.getExtent());

const reprojTilePixelRatio = Math.max.apply(
  null,
  tileGrid.getResolutions().map((r, z) => {
    const tileSize = toSize(tileGrid.getTileSize(z));
    const textureSize = source.getTileSize(z);
    return Math.max(textureSize[0] / tileSize[0], textureSize[1] / tileSize[1]);
  })
);

source.setTileGridForProjection(
  map_proj,
  createXYZ({
    extent: transformExtent(tif_view.extent, tif_view.projection, map_proj),
    maxResolution: resolution * reprojTilePixelRatio,
    maxZoom: maxZoom,
  })
);

[maxschmi]

Thanks a lot @mike-000
I think now I get the problem, that the tiles have to be somehow squares. I tested your code but the alignment offset is too big for my project. So with my tiff-layer the tiles can't be splitted into squares, unless they would be very small -> resulting in a huge TileGrid.

Thats a bit disapointing, but then maybe openlayer was not the tool to go. :-(

[maxschmi]

Just to show the resulting problem, I added some Polygons of the Grid shapes, where there border should be to the CodeSandbox.
As you can see the Borders ar not well aligned:

[mike-000]

Yes, that is to be expected as your non-square pixels are being drawn as square pixels in a canvas stitch context optimised for one dimension only before reprojection. Ideally your GeoTIFF should have better resolution, e.g. instead of one tile with 256 x 256 pixels have 2560 x 2560 pixels in 100 tiles.

While OpenLayers could potentially take an option to scale up before reprojecting (e.g. drawing 256 x 256 pixels as 2560 x 2560 pixels in a larger stitch context) that would use a lot of system resources compared with tiling for which GeoTIFFs can be optimised.

[maxschmi]

Oh wow. I'm not completly sure, why this works now, but this was the missing tip. Thanks @mike-000
After changing the Tiffs tiling and pixel resolution I found one setup that works perfectly: CodeSandbox

Interestingly it did only work with this specific pixel resolution (used now 4 times as many pixels as before). If I went up It sstoped working again.