COGLayer

[kylebarron]

COGLayer: A deck.gl layer to render a COG directly

Overview

A COG is a single GeoTIFF file that holds multiple resolutions of raster data. See https://cogeo.org for more info.

COGs are internally tiled. Thus they're a great fit for the TileLayer, which renders a tiled dataset. However the TileLayer can't be used directly, because it assumes a web-mercator tiling system. Even COGs that are in the web mercator projection can't be used directly, because the internal indexing doesn't match the global OSM indexing.

Hence there are three main steps I need to do to get this working:

1. Parse COG metadata to construct a 2D TileMatrixSet
2. Custom TileLayer that places an arbitrary matrix of tiles
3. Reproject data in the BitmapLayer

Parse COG metadata to construct TileMatrixSet

Basically working code:

Code

var fs = require('fs')
var geotiff = require('./node_modules/geotiff/dist/geotiff.bundle')
var buf = fs.readFileSync('m_3411861_ne_11_060_20180723_20190208.tif')

var tif;
geotiff.fromArrayBuffer(buf.buffer).then(x => {
  tif = x
})
var image;
tif.getImage().then(x => {
  image = x
})


/**
 * Get the geotransform of the image at a specific overview level
 *
 * @param  {Object} tif GeoTiff object
 * @param  {Number} ovrLevel overview level
 * @return {Object}          [description]
 */
function getGeotransform(tif, ovrLevel = 0) {
  var ifds = tif.fileDirectories;
  var firstIfd = ifds[0][0];

  // Calculate overview for source image
  var gt = [
    firstIfd.ModelPixelScale[0],
    0.0,
    firstIfd.ModelTiepoint[3],
    0.0,
    -firstIfd.ModelPixelScale[1],
    firstIfd.ModelTiepoint[4],
  ]

  // Decimate the geotransform if an overview is requested
  if (ovrLevel > 0) {
    var bounds = getBounds(gt, firstIfd)
    var ifd = ifds[ovrLevel][0];
    // TODO implement the translation and scale
    // Maybe use math.gl?
    gt = affine.Affine.translation(bounds[0], bounds[3]) * affine.Affine.scale(
        (bounds[2] - bounds[0]) / ifd.ImageWidth.value,
        (bounds[1] - bounds[3]) / ifd.ImageHeight.value,
    )
  }
  return gt;
}

function getBounds(gt, firstIfd) {
  tlx = gt[2]
  tly = gt[5]
  brx = tlx + (gt[0] * firstIfd.ImageWidth)
  bry = tly + (gt[4] * firstIfd.ImageLength)
  return [tlx, bry, brx, tly];
}


function create_tile_matrix_set(tif) {
  var matrices = [];
  var ifds = tif.fileDirectories;
  var firstIfd = ifds[0]

  for (var i = 0; i < ifds.length; i++) {
    var ifd = ifds[i][0];
    var gt = getGeotransform(tif, i);

    var matrix = {
      identifier: ifds.length - i - 1,
      topLeftCorner: [gt[2], gt[5]],
      tileWidth: ifd.TileWidth,
      tileHeight: ifd.TileLength,
      matrixWidth: Math.ceil(ifd.ImageWidth / ifd.TileWidth),
      matrixHeight: Math.ceil(ifd.ImageLength / ifd.TileLength),
      scaleDenominator: gt[0] / 0.28e-3,
    }
    matrices.push(matrix);
  }

  // Reverse order of array so that identifier 0 (lowest resolution) is first
  matrices.reverse();

  // Is this stable?
  var geoKeys = firstIfd[1];
  var epsg = geoKeys.ProjectedCSTypeGeoKey || geoKeys.GeographicTypeGeoKey

  var tms = {
    title: "Tile matrix for GeoTIFF",
    identifier: 'identifier or str(uuid.uuid4())',
    supportedCRS: "http://www.opengis.net" + `/def/crs/EPSG/0/${epsg}`,
    tileMatrix: matrices
  }

  return tms;
}

Note: you should be able to use math.gl where affine is used. If you look into the affine code, the source for the transformation matrices is pretty simple.

Custom TileLayer that supports TMS

While this could conceivably be included in deck.gl upstream, for now it would be better to work on it locally. The TMS provides a mapping from x, y, z values to locations in space. Note these are arbitrary locations, and should not be confused with the x, y, z indexing in the OSM system.

There are a few parts to this:

• Determining reprojection function
• Finding visibility

Determining reprojection function

First you need to find the reprojection string. So if the EPSG code in the geotiff is EPSG:26911, send a request to

https://epsg.io/26911.proj4

that returns

+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Then pass that to proj4.

Finding visibility

Deck.gl uses frustum culling to determine the visibility of tiles. Hence I need to modify the tile-2d-traversal.js script to work with an arbitrary TMS.

I'll need to overwrite getBoundingVolume to find the tile's bounding box in the local CRS, then reproject it. Note that getBoundingVolume deals in deck.gl's common space, this is a web mercator projection of world size 512, with the origin in the top left.

Note that one caveat to TMS is that child tiles don't have strict nesting. Hence I believe it’s not always true that for any tile t, it has only one parent. This makes visibility computations harder, as finding a tile's children is less straightforward. With OSM indexing, given any tile t, there are always exactly four children tiles. With non-strict children, you essentially need to make a pass over all tiles at that level, checking if the tile's bounding box intersects with the parent.

Therefore you'll also have to overwrite the children method to determine the tile's children.

Since Xiaoji is such a good engineer, you might have to overwrite only those two methods!

Reproject data in renderSubLayers

This is explored a bit in #37. Basically you'll need to reproject the vertices within createMesh from the source CRS to one deck can understand.

Note that the mesh already is created at a certain resolution that's derived from the zoom level.

Note that this will only be possible with the RasterLayer, not the RasterMeshLayer, as the mesh layer needs to combine elevation data, and that would only work if you could fetch elevation data in the same local coordinate system.

[kylebarron]

Note that one caveat to TMS is that child tiles don't have strict nesting

Note however that you have access to the coordinates of each tile in the local coordinate system. So if tile t has (local) coordinates [2, 2, 4, 4], then to find child tiles, you want to find tiles at level ("identifier" in TMS speak) n + 1 that intersect that box. Each level of the TMS has data:

{
    "type": "TileMatrixType",
    "identifier": "4",
    "scaleDenominator": 34884088.0626143,
    "topLeftCorner": [
        -9501965.72931276,
        20003931.4586255
    ],
    "tileWidth": 256,
    "tileHeight": 256,
    "matrixWidth": 8,
    "matrixHeight": 16
}

Since at each level the tile width and height must be regular, it should be possible to do a line or two of math to find the grid of tiles desired, without doing a loop. You just always need to know the coordinates in the local coordinate system of the current tile.

[kylebarron]

Note you'll probably also need masking to discard areas of the tiff that are nodata.

Note you'll need to pass tile width/height to the Texture constructor, and make mipmaps: false for non-power of two textures