Allow sampling data from neighboring pixels

examples/webgl-shaded-relief.css

@@ -0,0 +1,7 @@
+table.controls td {
+  padding: 2px 5px;
+}
+table.controls td:nth-child(3) {
+  text-align: right;
+  min-width: 3em;
+}

examples/webgl-shaded-relief.html

@@ -0,0 +1,32 @@
+---
+layout: example.html
+title: Shaded Relief (with WebGL)
+shortdesc: Calculate shaded relief from elevation data
+docs: >
+  <p>
+    For the shaded relief, a single tiled source of elevation data is used as input.
+    The shaded relief is calculated by the layer's <code>style</code> with a <code>color</code>
+    expression. The style variables are updated when the user drags one of the sliders. The
+    <code>band</code> operator is used to sample data from neighboring pixels for calculating slope and
+    aspect, which is done with the <code>['band', bandIndex, xOffset, yOffset]</code> syntax.
+  </p>
+tags: "webgl, shaded relief"
+---
+<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>

examples/webgl-shaded-relief.js

@@ -0,0 +1,97 @@
+import Map from '../src/ol/Map.js';
+import View from '../src/ol/View.js';
+import {OSM, XYZ} from '../src/ol/source.js';
+import {WebGLTile as TileLayer} from '../src/ol/layer.js';
+
+const variables = {};
+
+// The method used to extract elevations from the DEM.
+// In this case the format used is
+// red + green * 2 + blue * 3
+//
+// Other frequently used methods include the Mapbox format
+// (red * 256 * 256 + green * 256 + blue) * 0.1 - 10000
+// and the Terrarium format
+// (red * 256 + green + blue / 256) - 32768
+function elevation(xOffset, yOffset) {
+  return [
+    '+',
+    ['*', 256, ['band', 1, xOffset, yOffset]],
+    [
+      '+',
+      ['*', 2 * 256, ['band', 2, xOffset, yOffset]],
+      ['*', 3 * 256, ['band', 3, xOffset, yOffset]],
+    ],
+  ];
+}
+
+// Generates a shaded relief image given elevation data.  Uses a 3x3
+// neighborhood for determining slope and aspect.
+const halfPi = Math.PI / 2;
+const dp = ['*', 2, ['resolution']];
+const z0x = ['*', ['var', 'vert'], elevation(-1, 0)];
+const z1x = ['*', ['var', 'vert'], elevation(1, 0)];
+const dzdx = ['/', ['-', z1x, z0x], dp];
+const z0y = ['*', ['var', 'vert'], elevation(0, -1)];
+const z1y = ['*', ['var', 'vert'], elevation(0, 1)];
+const dzdy = ['/', ['-', z1y, z0y], dp];
+const slope = ['atan', ['^', ['+', ['^', dzdx, 2], ['^', dzdy, 2]], 0.5]];
+const rawAspect = ['atan', dzdy, ['-', 0, dzdx]];
+const aspect = [
+  'case',
+  ['>', rawAspect, halfPi],
+  ['+', halfPi, ['-', Math.PI * 2, rawAspect]],
+  ['-', halfPi, rawAspect],
+];
+const sunEl = ['*', Math.PI / 180, ['var', 'sunEl']];
+const sunAz = ['*', Math.PI / 180, ['var', 'sunAz']];
+const cosIncidence = [
+  '+',
+  ['*', ['sin', sunEl], ['cos', slope]],
+  ['*', ['*', ['cos', sunEl], ['sin', slope]], ['cos', ['-', sunAz, aspect]]],
+];
+const scaled = ['*', 255, cosIncidence];
+
+const shadedRelief = new TileLayer({
+  opacity: 0.3,
+  source: new XYZ({
+    url: 'https://{a-d}.tiles.mapbox.com/v3/aj.sf-dem/{z}/{x}/{y}.png',
+    crossOrigin: 'anonymous',
+  }),
+  style: {
+    variables: variables,
+    color: ['color', scaled, scaled, scaled],
+  },
+});
+
+const controlIds = ['vert', 'sunEl', 'sunAz'];
+controlIds.forEach(function (id) {
+  const control = document.getElementById(id);
+  const output = document.getElementById(id + 'Out');
+  function updateValues() {
+    output.innerText = control.value;
+    variables[id] = Number(control.value);
+  }
+  updateValues();
+  control.addEventListener('input', () => {
+    updateValues();
+    shadedRelief.updateStyleVariables(variables);
+  });
+});
+
+const map = new Map({
+  target: 'map',
+  layers: [
+    new TileLayer({
+      source: new OSM(),
+    }),
+    shadedRelief,
+  ],
+  view: new View({
+    extent: [-13675026, 4439648, -13580856, 4580292],
+    center: [-13615645, 4497969],
+    minZoom: 10,
+    maxZoom: 16,
+    zoom: 13,
+  }),
+});

src/ol/layer.js

@@ -13,3 +13,4 @@ export {default as Vector} from './layer/Vector.js';
 export {default as VectorImage} from './layer/VectorImage.js';
 export {default as VectorTile} from './layer/VectorTile.js';
 export {default as WebGLPoints} from './layer/WebGLPoints.js';
+export {default as WebGLTile} from './layer/WebGLTile.js';

src/ol/layer/WebGLTile.js

@@ -215,6 +215,10 @@ function parseStyle(style, bandCount) {
 
     varying vec2 v_textureCoord;
     uniform float ${Uniforms.TRANSITION_ALPHA};
+    uniform float ${Uniforms.TEXTURE_PIXEL_WIDTH};
+    uniform float ${Uniforms.TEXTURE_PIXEL_HEIGHT};
+    uniform float ${Uniforms.RESOLUTION};
+    uniform float ${Uniforms.ZOOM};
 
     ${uniformDeclarations.join('\n')}
 

src/ol/renderer/webgl/TileLayer.js

@@ -35,6 +35,10 @@ export const Uniforms = {
   TILE_TRANSFORM: 'u_tileTransform',
   TRANSITION_ALPHA: 'u_transitionAlpha',
   DEPTH: 'u_depth',
+  TEXTURE_PIXEL_WIDTH: 'u_texturePixelWidth',
+  TEXTURE_PIXEL_HEIGHT: 'u_texturePixelHeight',
+  RESOLUTION: 'u_resolution',
+  ZOOM: 'u_zoom',
 };
 
 export const Attributes = {
@@ -75,6 +79,23 @@ function addTileTextureToLookup(tileTexturesByZ, tileTexture, z) {
   tileTexturesByZ[z].push(tileTexture);
 }
 
+/**
+ *
+ * @param {import("../../PluggableMap.js").FrameState} frameState Frame state.
+ * @return {import("../../extent.js").Extent} Extent.
+ */
+function getRenderExtent(frameState) {
+  const layerState = frameState.layerStatesArray[frameState.layerIndex];
+  let extent = frameState.extent;
+  if (layerState.extent) {
+    extent = getIntersection(
+      extent,
+      fromUserExtent(layerState.extent, frameState.viewState.projection)
+    );
+  }
+  return extent;
+}
+
 /**
  * @typedef {Object} Options
  * @property {string} vertexShader Vertex shader source.
@@ -183,6 +204,9 @@ class WebGLTileLayerRenderer extends WebGLLayerRenderer {
    * @return {boolean} Layer is ready to be rendered.
    */
   prepareFrame(frameState) {
+    if (isEmpty(getRenderExtent(frameState))) {
+      return false;
+    }
     const source = this.getLayer().getSource();
     if (!source) {
       return false;
@@ -198,20 +222,9 @@ class WebGLTileLayerRenderer extends WebGLLayerRenderer {
   renderFrame(frameState) {
     this.preRender(frameState);
 
-    const layerState = frameState.layerStatesArray[frameState.layerIndex];
     const viewState = frameState.viewState;
-
-    let extent = frameState.extent;
-    if (layerState.extent) {
-      extent = getIntersection(
-        extent,
-        fromUserExtent(layerState.extent, viewState.projection)
-      );
-    }
-    if (isEmpty(extent)) {
-      return;
-    }
-
+    const layerState = frameState.layerStatesArray[frameState.layerIndex];
+    const extent = getRenderExtent(frameState);
     const tileLayer = this.getLayer();
     const tileSource = tileLayer.getSource();
     const tileGrid = tileSource.getTileGridForProjection(viewState.projection);
@@ -421,6 +434,19 @@ class WebGLTileLayerRenderer extends WebGLLayerRenderer {
 
         this.helper.setUniformFloatValue(Uniforms.TRANSITION_ALPHA, alpha);
         this.helper.setUniformFloatValue(Uniforms.DEPTH, depth);
+        this.helper.setUniformFloatValue(
+          Uniforms.TEXTURE_PIXEL_WIDTH,
+          tileSize[0]
+        );
+        this.helper.setUniformFloatValue(
+          Uniforms.TEXTURE_PIXEL_HEIGHT,
+          tileSize[1]
+        );
+        this.helper.setUniformFloatValue(
+          Uniforms.RESOLUTION,
+          viewState.resolution
+        );
+        this.helper.setUniformFloatValue(Uniforms.ZOOM, viewState.zoom);
 
         this.helper.drawElements(0, this.indices_.getSize());
       }

src/ol/source.js

@@ -5,6 +5,8 @@
 export {default as BingMaps} from './source/BingMaps.js';
 export {default as CartoDB} from './source/CartoDB.js';
 export {default as Cluster} from './source/Cluster.js';
+export {default as DataTile} from './source/DataTile.js';
+export {default as GeoTIFF} from './source/GeoTIFF.js';
 export {default as IIIF} from './source/IIIF.js';
 export {default as Image} from './source/Image.js';
 export {default as ImageArcGISRest} from './source/ImageArcGISRest.js';

src/ol/style/expressions.js

@@ -3,6 +3,7 @@
  * @module ol/style/expressions
  */
 
+import {Uniforms} from '../renderer/webgl/TileLayer.js';
 import {asArray, isStringColor} from '../color.js';
 
 /**
@@ -12,12 +13,13 @@ import {asArray, isStringColor} from '../color.js';
  * The following operators can be used:
  *
  * * Reading operators:
- *   * `['band', bandIndex]` fetches a pixel value from band `bandIndex` of the source's data. The first
- *     `bandIndex` of the source data is `1`. Fetched values are in the 0..1 range.
- *     {@link import("../source/TileImage.js").default} sources have 4 bands: red, green, blue and alpha.
- *     {@link import("../source/DataTile.js").default} sources can have any number of bands, depending on
- *     the underlying data source and
- *     {@link import("../source/GeoTIFF.js").Options configuration}.
+ *   * `['band', bandIndex, xOffset, yOffset]` For tile layers only. Fetches pixel values from band
+ *     `bandIndex` of the source's data. The first `bandIndex` of the source data is `1`. Fetched values
+ *     are in the 0..1 range. {@link import("../source/TileImage.js").default} sources have 4 bands: red,
+ *     green, blue and alpha. {@link import("../source/DataTile.js").default} sources can have any number
+ *     of bands, depending on the underlying data source and
+ *     {@link import("../source/GeoTIFF.js").Options configuration}. `xOffset` and `yOffset` are optional
+ *     and allow specifying pixel offsets for x and y. This is used for sampling data from neighboring pixels.
  *   * `['get', 'attributeName']` fetches a feature attribute (it will be prefixed by `a_` in the shader)
  *     Note: those will be taken from the attributes provided to the renderer
  *   * `['resolution']` returns the current resolution
@@ -34,6 +36,9 @@ import {asArray, isStringColor} from '../color.js';
  *   * `['%', value1, value2]` returns the result of `value1 % value2` (modulo)
  *   * `['^', value1, value2]` returns the value of `value1` raised to the `value2` power
  *   * `['abs', value1]` returns the absolute value of `value1`
+ *   * `['sin', value1]` returns the sine of `value1`
+ *   * `['cos', value1]` returns the cosine of `value1`
+ *   * `['atan', value1, value2]` returns `atan2(value1, value2)`. If `value2` is not provided, returns `atan(value1)`
  *
  * * Transform operators:
  *   * `['case', condition1, output1, ...conditionN, outputN, fallback]` selects the first output whose corresponding
@@ -416,7 +421,8 @@ Operators['band'] = {
     return ValueTypes.NUMBER;
   },
   toGlsl: function (context, args) {
-    assertArgsCount(args, 1);
+    assertArgsMinCount(args, 1);
+    assertArgsMaxCount(args, 3);
     const band = args[0];
     if (typeof band !== 'number') {
       throw new Error('Band index must be a number');
@@ -428,7 +434,22 @@ Operators['band'] = {
       // LUMINANCE_ALPHA - band 1 assigned to rgb and band 2 assigned to alpha
       bandIndex = 3;
     }
-    return `color${colorIndex}[${bandIndex}]`;
+    if (args.length === 1) {
+      return `color${colorIndex}[${bandIndex}]`;
+    } else {
+      const xOffset = args[1];
+      const yOffset = args[2] || 0;
+      assertNumber(xOffset);
+      assertNumber(yOffset);
+      const uniformName = Uniforms.TILE_TEXTURE_PREFIX + colorIndex;
+      return `texture2D(${uniformName}, v_textureCoord + vec2(${expressionToGlsl(
+        context,
+        xOffset
+      )} / ${Uniforms.TEXTURE_PIXEL_WIDTH}, ${expressionToGlsl(
+        context,
+        yOffset
+      )} / ${Uniforms.TEXTURE_PIXEL_HEIGHT}))[${bandIndex}]`;
+    }
   },
 };
 
@@ -570,6 +591,45 @@ Operators['abs'] = {
   },
 };
 
+Operators['sin'] = {
+  getReturnType: function (args) {
+    return ValueTypes.NUMBER;
+  },
+  toGlsl: function (context, args) {
+    assertArgsCount(args, 1);
+    assertNumbers(args);
+    return `sin(${expressionToGlsl(context, args[0])})`;
+  },
+};
+
+Operators['cos'] = {
+  getReturnType: function (args) {
+    return ValueTypes.NUMBER;
+  },
+  toGlsl: function (context, args) {
+    assertArgsCount(args, 1);
+    assertNumbers(args);
+    return `cos(${expressionToGlsl(context, args[0])})`;
+  },
+};
+
+Operators['atan'] = {
+  getReturnType: function (args) {
+    return ValueTypes.NUMBER;
+  },
+  toGlsl: function (context, args) {
+    assertArgsMinCount(args, 1);
+    assertArgsMaxCount(args, 2);
+    assertNumbers(args);
+    return args.length === 2
+      ? `atan(${expressionToGlsl(context, args[0])}, ${expressionToGlsl(
+          context,
+          args[1]
+        )})`
+      : `atan(${expressionToGlsl(context, args[0])})`;
+  },
+};
+
 Operators['>'] = {
   getReturnType: function (args) {
     return ValueTypes.BOOLEAN;