Edit

Shaded Relief

raster4 shadedrelief1

x
°
°

Calculate shaded relief from elevation data

This example uses a ol/source/Raster to generate data based on another source. The raster source accepts any number of input sources (tile or image based) and runs a pipeline of operations on the input data. The return from the final operation is used as the data for the output source.

In this case, a single tiled source of elevation data is used as input. The shaded relief is calculated in a single "image" operation. By setting operationType: 'image' on the raster source, operations are called with an ImageData object for each of the input sources. Operations are also called with a general purpose data object. In this example, the sun elevation and azimuth data from the inputs above are assigned to this data object and accessed in the shading operation. The shading operation returns an array of ImageData objects. When the raster source is used by an image layer, the first ImageData object returned by the last operation in the pipeline is used for rendering.

main.js
import 'ol/ol.css';
import Map from 'ol/Map';
import View from 'ol/View';
import {Image as ImageLayer, Tile as TileLayer} from 'ol/layer';
import {OSM, Raster, XYZ} from 'ol/source';

/**
 * Generates a shaded relief image given elevation data.  Uses a 3x3
 * neighborhood for determining slope and aspect.
 * @param {Array<ImageData>} inputs Array of input images.
 * @param {Object} data Data added in the "beforeoperations" event.
 * @return {ImageData} Output image.
 */
function shade(inputs, data) {
  var elevationImage = inputs[0];
  var width = elevationImage.width;
  var height = elevationImage.height;
  var elevationData = elevationImage.data;
  var shadeData = new Uint8ClampedArray(elevationData.length);
  var dp = data.resolution * 2;
  var maxX = width - 1;
  var maxY = height - 1;
  var pixel = [0, 0, 0, 0];
  var twoPi = 2 * Math.PI;
  var halfPi = Math.PI / 2;
  var sunEl = (Math.PI * data.sunEl) / 180;
  var sunAz = (Math.PI * data.sunAz) / 180;
  var cosSunEl = Math.cos(sunEl);
  var sinSunEl = Math.sin(sunEl);
  var pixelX,
    pixelY,
    x0,
    x1,
    y0,
    y1,
    offset,
    z0,
    z1,
    dzdx,
    dzdy,
    slope,
    aspect,
    cosIncidence,
    scaled;
  for (pixelY = 0; pixelY <= maxY; ++pixelY) {
    y0 = pixelY === 0 ? 0 : pixelY - 1;
    y1 = pixelY === maxY ? maxY : pixelY + 1;
    for (pixelX = 0; pixelX <= maxX; ++pixelX) {
      x0 = pixelX === 0 ? 0 : pixelX - 1;
      x1 = pixelX === maxX ? maxX : pixelX + 1;

      // determine elevation for (x0, pixelY)
      offset = (pixelY * width + x0) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      z0 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);

      // determine elevation for (x1, pixelY)
      offset = (pixelY * width + x1) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      z1 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);

      dzdx = (z1 - z0) / dp;

      // determine elevation for (pixelX, y0)
      offset = (y0 * width + pixelX) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      z0 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);

      // determine elevation for (pixelX, y1)
      offset = (y1 * width + pixelX) * 4;
      pixel[0] = elevationData[offset];
      pixel[1] = elevationData[offset + 1];
      pixel[2] = elevationData[offset + 2];
      pixel[3] = elevationData[offset + 3];
      z1 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);

      dzdy = (z1 - z0) / dp;

      slope = Math.atan(Math.sqrt(dzdx * dzdx + dzdy * dzdy));

      aspect = Math.atan2(dzdy, -dzdx);
      if (aspect < 0) {
        aspect = halfPi - aspect;
      } else if (aspect > halfPi) {
        aspect = twoPi - aspect + halfPi;
      } else {
        aspect = halfPi - aspect;
      }

      cosIncidence =
        sinSunEl * Math.cos(slope) +
        cosSunEl * Math.sin(slope) * Math.cos(sunAz - aspect);

      offset = (pixelY * width + pixelX) * 4;
      scaled = 255 * cosIncidence;
      shadeData[offset] = scaled;
      shadeData[offset + 1] = scaled;
      shadeData[offset + 2] = scaled;
      shadeData[offset + 3] = elevationData[offset + 3];
    }
  }

  return {data: shadeData, width: width, height: height};
}

var elevation = new XYZ({
  url: 'https://{a-d}.tiles.mapbox.com/v3/aj.sf-dem/{z}/{x}/{y}.png',
  crossOrigin: 'anonymous',
});

var raster = new Raster({
  sources: [elevation],
  operationType: 'image',
  operation: shade,
});

var map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new OSM(),
    }),
    new ImageLayer({
      opacity: 0.3,
      source: raster,
    }) ],
  view: new View({
    extent: [-13675026, 4439648, -13580856, 4580292],
    center: [-13615645, 4497969],
    minZoom: 10,
    maxZoom: 16,
    zoom: 13,
  }),
});

var controlIds = ['vert', 'sunEl', 'sunAz'];
var controls = {};
controlIds.forEach(function (id) {
  var control = document.getElementById(id);
  var output = document.getElementById(id + 'Out');
  control.addEventListener('input', function () {
    output.innerText = control.value;
    raster.changed();
  });
  output.innerText = control.value;
  controls[id] = control;
});

raster.on('beforeoperations', function (event) {
  // the event.data object will be passed to operations
  var data = event.data;
  data.resolution = event.resolution;
  for (var id in controls) {
    data[id] = Number(controls[id].value);
  }
});
index.html
<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="UTF-8">
    <title>Shaded Relief</title>
    <!-- Pointer events polyfill for old browsers, see https://caniuse.com/#feat=pointer -->
    <script src="https://unpkg.com/elm-pep"></script>
    <style>
      .map {
        width: 100%;
        height:400px;
      }
      table.controls td {
        padding: 2px 5px;
      }
      table.controls td:nth-child(3) {
        text-align: right;
        min-width: 3em;
      }
    </style>
  </head>
  <body>
    <div id="map" class="map"></div>
    <table class="controls">
      <tr>
        <td><label for="vert">vertical exaggeration:</label></td>
        <td><input id="vert" type="range" min="1" max="5" value="1"/></td>
        <td><span id="vertOut"></span> x</td>
      </tr>
      <tr>
        <td><label for="sunEl">sun elevation:</label></td>
        <td><input id="sunEl" type="range" min="0" max="90" value="45"/></td>
        <td><span id="sunElOut"></span> °</td>
      </tr>
      <tr>
        <td><label for="sunAz">sun azimuth:</label></td>
        <td><input id="sunAz" type="range" min="0" max="360" value="45"/></td>
        <td><span id="sunAzOut"></span> °</td>
      </tr>
    </table>
    <script src="main.js"></script>
  </body>
</html>
package.json
{
  "name": "shaded-relief",
  "dependencies": {
    "ol": "6.5.0"
  },
  "devDependencies": {
    "parcel": "^2.0.0-beta.1"
  },
  "scripts": {
    "start": "parcel index.html",
    "build": "parcel build --public-url . index.html"
  }
}