OpenLayers
OpenLayers
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.
import ImageTile from 'ol/source/ImageTile.js';
import Map from 'ol/Map.js';
import View from 'ol/View.js';
import {Image as ImageLayer, Tile as TileLayer} from 'ol/layer.js';
import {OSM, Raster} from 'ol/source.js';
/**
* 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) {
const elevationImage = inputs[0];
const width = elevationImage.width;
const height = elevationImage.height;
const elevationData = elevationImage.data;
const shadeData = new Uint8ClampedArray(elevationData.length);
const dp = data.resolution * 2;
const maxX = width - 1;
const maxY = height - 1;
const pixel = [0, 0, 0, 0];
const twoPi = 2 * Math.PI;
const halfPi = Math.PI / 2;
const sunEl = (Math.PI * data.sunEl) / 180;
const sunAz = (Math.PI * data.sunAz) / 180;
const cosSunEl = Math.cos(sunEl);
const sinSunEl = Math.sin(sunEl);
let pixelX,
pixelY,
x0,
x1,
y0,
y1,
offset,
z0,
z1,
dzdx,
dzdy,
slope,
aspect,
cosIncidence,
scaled;
function calculateElevation(pixel) {
// The method used to extract elevations from the DEM.
// In this case the format used is Terrarium
// red * 256 + green + blue / 256 - 32768
//
// Other frequently used methods include the Mapbox format
// (red * 256 * 256 + green * 256 + blue) * 0.1 - 10000
//
return pixel[0] * 256 + pixel[1] + pixel[2] / 256 - 32768;
}
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 * calculateElevation(pixel);
// 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 * calculateElevation(pixel);
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 * calculateElevation(pixel);
// 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 * calculateElevation(pixel);
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};
}
const elevation = new ImageTile({
url: 'https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png',
crossOrigin: 'anonymous',
maxZoom: 15,
attributions:
'<a href="https://github.com/tilezen/joerd/blob/master/docs/attribution.md" target="_blank">Data sources and attribution</a>',
});
const raster = new Raster({
sources: [elevation],
operationType: 'image',
operation: shade,
});
const map = new Map({
target: 'map',
layers: [
new TileLayer({
source: new OSM(),
}),
new ImageLayer({
opacity: 0.3,
source: raster,
}),
],
view: new View({
center: [-13615645, 4497969],
zoom: 13,
}),
});
const controlIds = ['vert', 'sunEl', 'sunAz'];
const controls = {};
controlIds.forEach(function (id) {
const control = document.getElementById(id);
const 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
const data = event.data;
data.resolution = event.resolution;
for (const id in controls) {
data[id] = Number(controls[id].value);
}
});
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title>Shaded Relief</title>
<link rel="stylesheet" href="node_modules/ol/ol.css">
<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 type="module" src="main.js"></script>
</body>
</html>{
"name": "shaded-relief",
"dependencies": {
"ol": "10.3.1"
},
"devDependencies": {
"vite": "^3.2.3"
},
"scripts": {
"start": "vite",
"build": "vite build"
}
}