Openlayers raster reprojection 柵格重投影原理分析

背景

OpenLayers 3 新版本可以以不一樣於服務器提供的座標系顯示來自WMS，WMTS，靜態圖像和許多其餘來源的柵格數據。圖像的地圖投影的轉換直接在Web瀏覽器中進行。任何Proj4js支持的座標參考系統中的視圖都是可能的，而且之前不兼容的圖層如今能夠組合和覆蓋。

Raster Reproject意義

能夠在瀏覽器端實現不一樣資源在不一樣投影下的轉換，再也不依賴於服務端處理，好比在Geoserver中指定投影類型，仍是至關給力的

效果圖

還能夠看如下效果對比圖

3857使用web Mercator 中國區域是方形，Reproject到4326時，就有些扁了，在相同的中心點，縮放級別下，地圖的視察範圍也不太一致，4326通過壓扁後，能顯示更多的範圍。

再來一張放大後的，變形明顯。

分析過程

以reprojection-image爲例說明，使用的是三角仿射變換triangle affine transformation。

三角形動態計算，以正好鋪滿當前視口

不一樣的縮放級別1

不一樣的縮放級別2

不一樣的縮放級別3

不一樣的縮放級別4

不一樣的比例尺下會更新三角格網的大小，以剛好平分地圖視口，至少兩個三角

每一個三角形單獨進行Reproject，不一樣三角之間互不影響

對三角格網渲染進行Degbu過濾後示意

渲染其中一部分三角

這個圖能方便明白Reproject原理

• 首先劃分三角格網
• 根據位置，計算每一個三角網對應的原始raster位置和區域，內部使用數學方法，求出轉換參數，對每一個三角網進行轉換，這裏其實並無對每一個像素進行Reproject，而是以三角爲最小單位，進行Reproject, 理論上只要三角足夠小，經過積分是能夠達到一樣的效果的，相比像素角度的處理效率也會高不少
• 將多個三角進行無緣拼接

核心實現

三角Reproject

源碼ol/reproj.js, render funcction

/** * Renders the source data into new canvas based on the triangulation. * * @param {number} width Width of the canvas. * @param {number} height Height of the canvas. * @param {number} pixelRatio Pixel ratio. * @param {number} sourceResolution Source resolution. * @param {import("./extent.js").Extent} sourceExtent Extent of the data source. * @param {number} targetResolution Target resolution. * @param {import("./extent.js").Extent} targetExtent Target extent. * @param {import("./reproj/Triangulation.js").default} triangulation * Calculated triangulation. * @param {Array<{extent: import("./extent.js").Extent, * image: (HTMLCanvasElement|HTMLImageElement|HTMLVideoElement)}>} sources * Array of sources. * @param {number} gutter Gutter of the sources. * @param {boolean=} opt_renderEdges Render reprojection edges. * @return {HTMLCanvasElement} Canvas with reprojected data. */
export function render(width, height, pixelRatio, sourceResolution, sourceExtent, targetResolution, targetExtent, triangulation, sources, gutter, opt_renderEdges) {

  var context = createCanvasContext2D(Math.round(pixelRatio * width),
    Math.round(pixelRatio * height));

  if (sources.length === 0) {
    return context.canvas;
  }

  context.scale(pixelRatio, pixelRatio);

  var sourceDataExtent = createEmpty();
  sources.forEach(function(src, i, arr) {
    extend(sourceDataExtent, src.extent);
  });

  var canvasWidthInUnits = getWidth(sourceDataExtent);
  var canvasHeightInUnits = getHeight(sourceDataExtent);
  
  // 加載image的內部canvas
  var stitchContext = createCanvasContext2D(
    Math.round(pixelRatio * canvasWidthInUnits / sourceResolution),
    Math.round(pixelRatio * canvasHeightInUnits / sourceResolution));

  var stitchScale = pixelRatio / sourceResolution;

  sources.forEach(function(src, i, arr) {
    var xPos = src.extent[0] - sourceDataExtent[0];
    var yPos = -(src.extent[3] - sourceDataExtent[3]);
    var srcWidth = getWidth(src.extent);
    var srcHeight = getHeight(src.extent);

    stitchContext.drawImage(
      src.image,
      gutter, gutter,
      src.image.width - 2 * gutter, src.image.height - 2 * gutter,
      xPos * stitchScale, yPos * stitchScale,
      srcWidth * stitchScale, srcHeight * stitchScale);
  });

  var targetTopLeft = getTopLeft(targetExtent);

  // 對三角進行循環處理
  triangulation.getTriangles().forEach(function(triangle, i, arr) {

    /* Calculate affine transform (src -> dst) * Resulting matrix can be used to transform coordinate * from `sourceProjection` to destination pixels. * * To optimize number of context calls and increase numerical stability, * we also do the following operations: * trans(-topLeftExtentCorner), scale(1 / targetResolution), scale(1, -1) * here before solving the linear system so [ui, vi] are pixel coordinates. * * Src points: xi, yi * Dst points: ui, vi * Affine coefficients: aij * * | x0 y0 1 0 0 0 | |a00| |u0| * | x1 y1 1 0 0 0 | |a01| |u1| * | x2 y2 1 0 0 0 | x |a02| = |u2| * | 0 0 0 x0 y0 1 | |a10| |v0| * | 0 0 0 x1 y1 1 | |a11| |v1| * | 0 0 0 x2 y2 1 | |a12| |v2| */
    var source = triangle.source;
    var target = triangle.target;
    var x0 = source[0][0], y0 = source[0][1];
    var x1 = source[1][0], y1 = source[1][1];
    var x2 = source[2][0], y2 = source[2][1];
    var u0 = (target[0][0] - targetTopLeft[0]) / targetResolution;
    var v0 = -(target[0][1] - targetTopLeft[1]) / targetResolution;
    var u1 = (target[1][0] - targetTopLeft[0]) / targetResolution;
    var v1 = -(target[1][1] - targetTopLeft[1]) / targetResolution;
    var u2 = (target[2][0] - targetTopLeft[0]) / targetResolution;
    var v2 = -(target[2][1] - targetTopLeft[1]) / targetResolution;

    // Shift all the source points to improve numerical stability
    // of all the subsequent calculations. The [x0, y0] is used here.
    // This is also used to simplify the linear system.
    var sourceNumericalShiftX = x0;
    var sourceNumericalShiftY = y0;
    x0 = 0;
    y0 = 0;
    x1 -= sourceNumericalShiftX;
    y1 -= sourceNumericalShiftY;
    x2 -= sourceNumericalShiftX;
    y2 -= sourceNumericalShiftY;

    var augmentedMatrix = [
      [x1, y1, 0, 0, u1 - u0],
      [x2, y2, 0, 0, u2 - u0],
      [0, 0, x1, y1, v1 - v0],
      [0, 0, x2, y2, v2 - v0]
    ];
    var affineCoefs = solveLinearSystem(augmentedMatrix);
    if (!affineCoefs) {
      return;
    }
    context.save();
    context.beginPath();
    var centroidX = (u0 + u1 + u2) / 3;
    var centroidY = (v0 + v1 + v2) / 3;
    var p0 = enlargeClipPoint(centroidX, centroidY, u0, v0);
    var p1 = enlargeClipPoint(centroidX, centroidY, u1, v1);
    var p2 = enlargeClipPoint(centroidX, centroidY, u2, v2);

    // 設置三角的切割範圍，只顯示這個範圍內的圖片，多個三角拼接起來就是總體圖像
    context.moveTo(p1[0], p1[1]);
    context.lineTo(p0[0], p0[1]);
    context.lineTo(p2[0], p2[1]);
    context.clip();

    // 最關鍵的三個方法，經過一系列轉換，保證在相應的比例尺下在當前範圍內
    // 顯示正確的Reproject圖像
    // 直接設置下面3個參數，很差容易能達到想要的效果
    // 能夠經過將相同的參數輸出在本地，進行調試，觀察三角內的圖像
    context.transform(
      affineCoefs[0], affineCoefs[2], affineCoefs[1], affineCoefs[3], u0, v0);

    context.translate(sourceDataExtent[0] - sourceNumericalShiftX,
      sourceDataExtent[3] - sourceNumericalShiftY);

    context.scale(sourceResolution / pixelRatio,
      -sourceResolution / pixelRatio);

    context.drawImage(stitchContext.canvas, 0, 0);
    context.restore();
  });

  // 調試使用，是否顯示三角
  if (opt_renderEdges) {
    context.save();

    context.strokeStyle = 'black';
    context.lineWidth = 1;

    triangulation.getTriangles().forEach(function(triangle, i, arr) {
      var target = triangle.target;
      var u0 = (target[0][0] - targetTopLeft[0]) / targetResolution;
      var v0 = -(target[0][1] - targetTopLeft[1]) / targetResolution;
      var u1 = (target[1][0] - targetTopLeft[0]) / targetResolution;
      var v1 = -(target[1][1] - targetTopLeft[1]) / targetResolution;
      var u2 = (target[2][0] - targetTopLeft[0]) / targetResolution;
      var v2 = -(target[2][1] - targetTopLeft[1]) / targetResolution;

      context.beginPath();
      context.moveTo(u1, v1);
      context.lineTo(u0, v0);
      context.lineTo(u2, v2);
      context.closePath();
      context.stroke();
    });

    context.restore();
  }
  return context.canvas;
}
複製代碼

我以其中一個三角的計算參數爲示例，能夠看到參數仍是比較複雜的

context.moveTo(1136, 199);
context.lineTo(850, 200);
context.lineTo(1136, 401);
context.clip();

context.transform(
  0.011129369548469296, 0.00006403305449769778, 0.00009456173415458175, -0.0111338055978892, 851.2500000000001, 200.00000000000063);

context.translate(-413988.7802610396,
  835088.5497620346);

context.scale(351.73160173160176,
  -351.73160173160176);

context.drawImage(stitchContext.canvas, 0, 0);
context.restore();
複製代碼

其中，進行仿射變換的轉換過程及數學原理，還須要進一步研究。

參考

• tutorials-raster-reprojection
• reprojection-image