Edit

Interpolating Contours from a DEM

image-interpolation1 contours1 raster6 mouse-position2 mapzen2

Example of WebGL data interpolation

Example of data resampling when using raster DEM (digital elevation model) data with WebGL. Mapzen Terrarium tiles with gutter are used which allow interpolation across tile edges. One buffered tile is produced by loading nine adjacent tiles. The data is interpolated by calculating the values represented by the RGB bands of the pixels (Float32 data tiles are used to improve precision on mobile devices) and the interpolated values are styled back as RGB so they can be accessed by a ol/source/Raster. One raster source is used to plot contours where adjacent interpolated pixels cross the contour interval. A second unseen raster source allows a custom mouse position control to display interpolated elevations obtained from the rendered pixels.

The DEM source has a maximum zoom level of 15 while the view and background layer can zoom in much further. The data values of the original source are shown for comparison. Notice how, particularly on steep ground, the interpolated values and contours change smoothly while the original data values are stepped.

main.js
import Map from 'ol/Map.js';
import View from 'ol/View.js';
import MousePosition from 'ol/control/MousePosition.js';
import {toStringHDMS} from 'ol/coordinate.js';
import {Image as ImageLayer, WebGLTile as TileLayer} from 'ol/layer.js';
import {fromLonLat, toLonLat} from 'ol/proj.js';
import {DataTile, OSM, Raster as RasterSource} from 'ol/source.js';

const attribution =
  '<a href="https://github.com/tilezen/joerd/blob/master/docs/attribution.md" target="_blank">Data sources and attribution</a>';

const calculateElevation = function (pixel) {
  if (pixel[3]) {
    return -32768 + (pixel[0] * 256 + pixel[1] + pixel[2] / 256);
  }
};

// Use Float32 interpolation where supported for best results on mobile devices.
const calculateElevationFromData = function (pixel) {
  return pixel[0];
};
const elevation = ['band', 1];
const tileSize = 256;
const gutter = 1;
const canvas = document.createElement('canvas');
canvas.width = tileSize * 3;
canvas.height = tileSize * 3;
const context = canvas.getContext('2d', {willReadFrequently: true});

const source = new DataTile({
  attributions: attribution,
  tileSize: tileSize,
  gutter: gutter,
  maxZoom: 15,
  interpolate: true,
  wrapX: true,
  loader: (z, x, y) => {
    const promises = [];
    for (let i = 0; i < 3; ++i) {
      for (let j = 0; j < 3; ++j) {
        promises.push(
          new Promise((resolve, reject) => {
            const maxY = 2 ** z;
            const yy = y + j - 1;
            if (yy < 0 || yy >= maxY) {
              return resolve();
            }
            const maxX = 2 ** z;
            const xx = (((x + i - 1) % maxX) + maxX) % maxX;
            const image = new Image();
            image.crossOrigin = '';
            image.addEventListener('error', () => reject());
            image.addEventListener('load', () => resolve(image));
            image.src = `https://s3.amazonaws.com/elevation-tiles-prod/terrarium/${z}/${xx}/${yy}.png`;
          }),
        );
      }
    }
    return Promise.all(promises).then((images) => {
      context.clearRect(0, 0, canvas.width, canvas.height);
      for (let i = 0; i < 3; ++i) {
        for (let j = 0; j < 3; ++j) {
          const image = images.shift();
          if (image) {
            context.drawImage(image, i * tileSize, j * tileSize);
          }
        }
      }
      const data = context.getImageData(
        tileSize - gutter,
        tileSize - gutter,
        tileSize + 2 * gutter,
        tileSize + 2 * gutter,
      ).data;
      const pixels = data.length / 4;
      const floatData = new Float32Array(data.buffer);
      for (let i = 0, j = 0; i < pixels; ) {
        floatData[i++] = calculateElevation(data.slice(j, (j += 4)));
      }
      return floatData;
    });
  },
});
const pixelValue = ['*', ['+', elevation, 32768], 256];

const style = {
  color: [
    'array',
    ['/', ['floor', ['/', pixelValue, 256 * 256]], 255],
    ['/', ['floor', ['/', ['%', pixelValue, 256 * 256], 256]], 255],
    ['/', ['%', pixelValue, 256], 255],
    1,
  ],
};

/// duplicate layers as one layer shared by two raster sources causes rendering issues
const elevation1 = new TileLayer({
  source: source,
  style: style,
});
const elevation2 = new TileLayer({
  source: source,
  style: style,
});

const contours = function (inputs, data) {
  const elevationImage = inputs[0];
  const width = elevationImage.width;
  const height = elevationImage.height;
  const elevationData = elevationImage.data;
  const pixel = [0, 0, 0, 0];
  const contourData = new Uint8ClampedArray(elevationData.length);
  const interval = data.interval;

  let offset, pixelY, pixelX;
  for (pixelY = 0; pixelY < height; ++pixelY) {
    for (pixelX = 0; pixelX < width; ++pixelX) {
      offset = (pixelY * width + Math.max(pixelX - 1, 0)) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      const leftElevation = calculateElevation(pixel);

      offset = (pixelY * width + Math.min(pixelX + 1, width - 1)) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      const rightElevation = calculateElevation(pixel);

      offset = (Math.max(pixelY - 1, 0) * width + pixelX) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      const topElevation = calculateElevation(pixel);

      offset = (Math.min(pixelY + 1, height - 1) * width + pixelX) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      const bottomElevation = calculateElevation(pixel);

      offset = (pixelY * width + pixelX) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      const centerElevation = calculateElevation(pixel);

      if (
        leftElevation !== undefined &&
        rightElevation !== undefined &&
        topElevation !== undefined &&
        bottomElevation !== undefined &&
        centerElevation !== undefined
      ) {
        const contour = Math.floor(centerElevation / interval);
        if (
          contour !== Math.floor(leftElevation / interval) ||
          contour !== Math.floor(rightElevation / interval) ||
          contour !== Math.floor(topElevation / interval) ||
          contour !== Math.floor(bottomElevation / interval)
        ) {
          if (
            centerElevation > 0 &&
            leftElevation > 0 &&
            rightElevation > 0 &&
            topElevation > 0 &&
            bottomElevation > 0
          ) {
            contourData[offset] = 0xe0;
            contourData[offset + 1] = 0x94;
            contourData[offset + 2] = 0x5e;
            contourData[offset + 3] = 255;
          } else {
            contourData[offset] = 0x00;
            contourData[offset + 1] = 0xa9;
            contourData[offset + 2] = 0xca;
            contourData[offset + 3] = 255;
          }
        }
      }
    }
  }
  return {data: contourData, width: width, height: height};
};

const contourSource = new RasterSource({
  sources: [elevation1],
  operationType: 'image',
  operation: contours,
  lib: {
    calculateElevation: calculateElevation,
  },
  resolutions: null,
});

contourSource.on('beforeoperations', function (event) {
  const data = event.data;
  if (event.resolution < 5) {
    data.interval = 10;
  } else if (event.resolution < 25) {
    data.interval = 50;
  } else if (event.resolution < 50) {
    data.interval = 100;
  } else if (event.resolution < 250) {
    data.interval = 500;
  } else {
    data.interval = 1000;
  }
});

const contourLayer = new ImageLayer({
  source: contourSource,
  opacity: 0.5,
});

const elevationLayer = new ImageLayer({
  source: new RasterSource({
    sources: [elevation2],
    operation: function (pixels) {
      return pixels[0];
    },
    resolutions: null,
  }),
  opacity: 0,
});

const dataLayer = new TileLayer({
  source: source,
  style: {color: ['array', 0, 0, 0, 0]},
});

dataLayer.once('postrender', function (event) {
  const gl = event.context;
  if (!gl.getSupportedExtensions().includes('OES_texture_float_linear')) {
    alert('Device does not support float interpolation');
  }
});

const map = new Map({
  layers: [
    new TileLayer({
      source: new OSM(),
    }),
    dataLayer,
    contourLayer,
    elevationLayer,
  ],
  target: 'map',
  view: new View({
    center: fromLonLat([-78.8175, -1.469167]),
    zoom: 17,
    maxZoom: 21,
  }),
});

const mousePositionControl = new MousePosition({
  className: 'custom-mouse-position',
  target: 'info',
  wrapX: false,
  coordinateFormat: function (coordinate) {
    let position = 'Position ' + toStringHDMS(toLonLat(coordinate)) + '<br>';
    let pixel = elevationLayer.getData(map.getPixelFromCoordinate(coordinate));
    if (pixel) {
      const elevation = calculateElevation(pixel);
      if (elevation !== undefined) {
        position += 'Elevation ' + elevation.toFixed(1) + ' meters';
      }
    }
    position += '<br><br>';

    pixel = dataLayer.getData(map.getPixelFromCoordinate(coordinate));
    if (pixel) {
      const elevation = calculateElevationFromData(pixel);
      if (elevation !== undefined) {
        position += 'Data value ' + elevation.toFixed(1) + ' meters';
      }
    }
    position += '<br>';
    return position;
  },
  placeholder: '<br><br><br><br>',
});
map.addControl(mousePositionControl);
index.html
<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="UTF-8">
    <title>Interpolating Contours from a DEM</title>
    <link rel="stylesheet" href="node_modules/ol/ol.css">
    <style>
      .map {
        width: 100%;
        height: 400px;
      }
    </style>
  </head>
  <body>
    <div id="map" class="map"></div>
    <div id="info"></div>

    <script type="module" src="main.js"></script>
  </body>
</html>
package.json
{
  "name": "contour-interpolation",
  "dependencies": {
    "ol": "10.4.1-dev"
  },
  "devDependencies": {
    "vite": "^3.2.3"
  },
  "scripts": {
    "start": "vite",
    "build": "vite build"
  }
}