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:
- Parse COG metadata to construct a 2D TileMatrixSet
- Custom TileLayer that places an arbitrary matrix of tiles
- 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.