D3js raster tools docs : Shaded relief

The shaded relief technique is a method for representing the topography which is prettier and intuitive.

It’s possible to get a simple shaded relief representation using only DEM data. I based all the calculations on this example made with python.

You can find the whole code here.

The way to color the different height values is the same as in the drawing raster data section, but clipping the image:

context.clip();
context.drawImage(canvasRaster.node(), 0, 0);
context.drawImage(canvasShaded.node(), 0, 0);
context.restore();

Note that the canvasShaded is pasted twice. The first time without clipping and adding a white rectangle to make the background smoother. The second one, with the clipping, to show a stronger effect on the colored zone

Calculating the shaded relief

To calculate the shaded relief, a hidden canvas is created with this data:

var azimuth = 315;
var angleAltitude = 45;
var azimuthrad = azimuth*Math.PI / 180;
var altituderad = angleAltitude*Math.PI / 180;

Changing these parameters, the relief will change its aspect
- azimuth is the angle where the sun is, calculated from the north. Note that in the northern hemisphere, this angle is impossible, but it’s more natural.
- angleAltitude is the position of the sun from the horizon. Increasing the angle will make the image to have less contrast

var shadedData = new Array(image.getHeight());
for (var j = 0; j<image.getHeight(); j++){
  shadedData[j] = new Array(image.getWidth());
  for (var i = 0; i<image.getWidth(); i++){
    var gradX, gradY;
    if(i==0) gradX = altData[j][i+1] - altData[j][i];
    else if(i==image.getWidth()-1) gradX = altData[j][i] - altData[j][i-1];
    else gradX = (altData[j][i+1] - altData[j][i])/2 + (altData[j][i] - altData[j][i-1])/2;

    if(j==0) gradY = altData[j+1][i] - altData[j][i];
    else if(j==image.getHeight()-1) gradY = altData[j][i] - altData[j-1][i];
    else gradY = (altData[j+1][i] - altData[j][i])/2 + (altData[j][i] - altData[j-1][i])/2;

    var slope = Math.PI/2 - Math.atan(Math.sqrt(gradX*gradX + gradY*gradY));
    var aspect = Math.atan2(-gradY, gradX);

    shadedData[j][i] = Math.sin(altituderad) * Math.sin(slope)
      + Math.cos(altituderad) * Math.cos(slope)
      * Math.cos(azimuthrad - aspect);
  }
}

Here is where the data is really calculated
The loop iterates all the pixels, as always when reading rasters, but the stored value is calculated
The slope and aspect need the gradients. Those are calculated with the central difference method, since is the one that gave me the same results as the Numpy gradient function
- Since the borders don’t have both side pixels, a check has to be done and apply the correct case
- The central difference, in this case, is a simple average of the left and right differences
- The gradient is calculated in the x and y directions
The slope is simply the arctan of the module of the gradient in both directions
The aspect is the direction of the slope. It’s important to note which component and sign goes first.
Once the slope and aspect are calculated, the shaded relief data is directly calculated

var canvasShaded = d3.select("body").append("canvas")
    .attr("width", width)
    .attr("height", height)
    .style("display","none");

var contextShaded = canvasShaded.node().getContext("2d");

var idShaded = contextShaded.createImageData(width,height);
var dataShaded = idShaded.data;
var posShaded = 0;
for(var j = 0; j<height; j++){
  for(var i = 0; i<width; i++){
    var pointCoords = projection.invert([i,j]);
    var px = invGeoTransform[0] + pointCoords[0]* invGeoTransform[1];
    var py = invGeoTransform[3] + pointCoords[1] * invGeoTransform[5];

    var shadedValue;
    if(Math.floor(px) >= 0 && Math.ceil(px) < image.getWidth() && Math.floor(py) >= 0 && Math.ceil(py) < image.getHeight()){
      shadedValue = 255*(1+shadedData[Math.floor(py)][Math.floor(px)])/2;

    } else {
      shadedValue = 255;
    }
    dataShaded[posShaded]   = shadedValue;
    dataShaded[posShaded+1]   = shadedValue;
    dataShaded[posShaded+2]   = shadedValue;
    dataShaded[posShaded+3]   = 255 - shadedValue;

    posShaded=posShaded+4;
  }
}

contextShaded.putImageData( idShaded, 0, 0);

Once the shaded relief data is known, the image of the shaded relief is calculated taking the closest pixel for each image pixel
Note that the alpha value is the inverse of the calculated hillshade value, as explained here. Using this, the shades will be much sharper, and the flat zones won’t add noise to the colored image

D3js raster tools docs Draw raster data using D3js in the browser or nodejs

Shaded relief

Calculating the shaded relief