[OpenLayers] XYZ 방식의 타일 DEM으로 음영기복도(Shade) 레이어 만들기

타일 방식으로 제공되는 DEM을 이용해 음영기복도를 OpenLayers에서 동적으로 생성해 표시하는 내용을 설명한다.

음영기복도가 무었인지 이 글에서 작성할 코드를 통해 생성되는 최종 결과물을 통해 살펴보자. 아래와 같다.

DEM 데이터를 이용한 음영기복도의 생성은 동적으로 이루어지는데, DEM 데이터는 XYZ 방식의 타일맵으로 제공된다. 이 XYZ 방식의 타일맵을 통해 높이 값을 얻을 수 있다. 이렇게 얻은 높이값을 이용해 음영기복도를 위한 값을 동적으로, 즉 실행중에 계산하여 이미지를 만들어 낸다. 만들어진 이미지는 ImageLayer 클래스의 객체로써 레이어로 처리될 수 있으며, 레이어는 지도의 구성 단위이다. 또한 음영기복도 생성을 위한 연산을 Raster 연산이라고 한다. Raster 연산을 위해서는 Raster라는 데이터소스 클래스가 필요하다. 이처럼 지도(Map), ImageLayer, Raster, XYZ 방식의 타일맵 소스에 대한 클래스 관계도는 아래와 같다.

자, 이제 앞서 언급한 객체를 생성하고 조합해 보자. 먼저 필요한 모듈을 선언한다.

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

그리고 타일맵 방식으로 서비스 되는 DEM 데이터 소스에 대한 객체는 아래와 같다.

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

위의 elevation 데이터소스에 대한 Raster 연산 수행이 가능한 Raster 데이터소스 객체는 아래와 같다.

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

Raster 연산에 대한 함수는 operation 속성으로 지정된 shade 함수이다. 이 shade 함수가 바로 음영기복도를 계산하여 그 결과를 이미지로 반환하는 함수인데, 아래와 같다. (음영기복도 생성을 위한 알고리즘에 대한 설명은 생략함)

function shade(inputs, data) {
    var elevationImage = inputs[0]; // 첫번째(0) 데이터소스로 elevation 객체
    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(data.sunEl);
    var sinSunEl = Math.sin(data.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};
}

위 코드의 shade 함수의 인자 중 data는 Raster 객체의 beforeoperations 이벤트에 의해 추가적으로 속성값이 추가될 수 있는데 아래와 같다.

raster.on('beforeoperations', function(event) {
    var data = event.data;

    data.resolution = event.resolution;
    data.vert = 1; // 음영기복도의 과고감 정도
    data.sunEl = 45; // 태양의 Elevation(단위: Degree)
    data.sunAz = 45; // 태양의 Azimuth(단위: Degree)
});

beforeoperations 이벤트는 Raster 연산이 실행되기 직전에 호출되는 함수이며, Raster 연산은 Raster의 데이터소스, 즉 XYZ를 통한 DEM 데이터를 받았을 때 호출된다.

이제 음영기복도에 대한 데이터소스에 대한 모든 정의가 완성되었으므로, 이를 화면에 표시하기 위한 레이어를 아래처럼 생성한다.

var imageLayer = new ImageLayer({ source: raster });

마지막으로 지도 객체를 생성하고 위의 imageLayer를 레이어로써 아래 코드처럼 지도에 추가한다.

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