3d Tile to Raster - jgoffeney/Cesium4Unreal GitHub Wiki

Back

Elevation Data

To convert a 3d tile to an elevation raster there are a few steps. The general idea is to render the model to a framebuffer using something like OpenGL. The model vertices are transformed within a vertex shader and the interpolated elevation values are written via a fragment shader. However the 3d tile GLTF data is stored in a transformed Earth Centered Earth Fixed (ECEF) format which is a cartesian format so directly rendering a model of a section of the earth would contain the Earth curvature. To create an actual digital elevation map (DEM) a planar format like cartographic of latitude, longitude and altitude should be used.

Restore GLTF ECEF

Geospatial vertex data stored as a GTLF model has been transformed and has to be restored to the original format. At a minimum GLTF expect the Y axis to be up while elevation data typically uses the Z axis. Additionally the model root node may contain a transform matrix (or a rotation and translation) which is used to restore the data. Finally there may be in the scene's "extensionsUsed" section the CESIUM_RTC tag which indicates an additional field containing three values that were used to translate the original data.

I have worked with data that contains a node matrix so the complete transform matrix is created by:

fixECEFMat4x4 = yUpToZUpMat4x4 * nodeTransformMat4x4

To get the original ECEF vertices then multiply the transform matrix by the GLTF vertex expanded to 4 elements as (x, y, z, 1).

originECEFVec4 = fixECEFMat4x4 * gltfEVEFVec4

This is described within the Transforms section of the 3d tiles specification.

Node Matrix

The GLTF Node JSON object may contain a matrix field containing either 1 or 16 values. The 16 values are a 4x4 transformation matrix stored in column major format (every 4 values define a matrix column). For example the values below generates the matrix.

-0.441601 0.824768 0.353195 0 -0.624602 0 -0.780943 0 -0.644097 -0.565472 0.515152 0 4.11305e+06 3.5868e+06 -3.28964e+06 1

|-0.441601|-0.624602|-0.644097| 4.11305e+06|
| 0.824768|        0|-0.565472|  3.5868e+06|
| 0.353195|-0.780943| 0.515152|-3.28964e+06|
|        0|        0|        0|           1|

Y Up to Z Up

The matrix below is used to convert the data back to having the Z axis as "up" (or toward the camera in OpenGL).

| 1| 0| 0| 0|
| 0| 0|-1| 0|
| 0| 1| 0| 0|
| 0| 0| 0| 1|

ECEF to LLA

Once the ECEF data has been restored then straight forward equations can be applied to convert to cartographic coordinates as latitude, longitude and altitude (LLA). The outputs are in radians and meters.

This is the closed form solution from Datum Transformations of GPS Positions.

Constants

  • WGS 84 ellipsoid semi-major axis
    • a = 6378137
  • WGS 84 ellipsoid semi-major axis
    • b = 6356752.314
  • WGS 84 ellipsoid first eccentricity squared
    • ecc1Sqrd = 0.00669438
  • WGS 84 ellipsoid second eccentricity squared
    • ecc2Sqrd = 0.006739497

Inputs

  • ECEF Vertex
    • X, Y, Z

Formulas

// Compute Longitude
longitude = atan2(Y, X)

// Compute Latitude
p = sqrt(X*X + Y*Y)
theta = atan2(Z*a, p*b)
sinTheta = sin(theta)
cosTheta = cos(theta)
latY = Z + ecc2Sqrd * b * sinTheta * sinTheta * sinTheta
latX = p - ecc1Sqrd * a * cosTheta * cosTheta * cosTheta

latitude = atan(latY, latX)       

// Compute Altitude
sinLat = sin(latitude)
N = a / sqrt(1.0 - ecc1Sqrd * sinLat * sinLat)

altitude = (p / cos(latitude)) - N

LLA Vertices to Raster

When using GLSL the general flow is vertex data is supplied to a vertex shader to be transformed and before the shader returns the vertex positions have to be converted into clip space in which the X, Y and Z axes have limits of -1.0 to 1.0 values. Any vertex that lies outside of these ranges are discarded so for a vertex to pass it has to be transformed by a projection matrix. An orthographic projection is used in this case since a perspective projection would warp the data. Within the fragment shader the clip coordinates are automatically mapped to the target array of a window or framebuffer (for an array of 1000 pixels then -1.0 maps to pixel 0 and 1.0 maps to pixel 999).

The desired geospatial bounds for a 3d tile are defined within its JSON file as radians and meters. Therefore the orthographic projection needs to map geospatial vertices to the clip space limits. In the vertex shader once the cartographic coordinates are computed are are transformed via a matrix which combines the projection matrix and a translation matrix.

I generate the orthographic matrix by first converting the 3d tile bounds latitude and longitude to degrees. Then I translate the bounds values by computing the average and subtracting it from the minimum value so for each translated range the minimum value = -(maximum value) which makes them centered at 0 (like with the clip space).

The translated bounds values are used to generate the orthographic projection and the negative average bounds values are used for a translation bounds. The clip transformation matrix for the vertex shader's LLA values is:

clipMat4x4 = orthoMat4x4 * translateMat4x4

Orthographic Projection

The formulas for creating an orthographic projection matrix is as shown below.

|A|0|0|D|
|0|B|0|E|
|0|0|C|F|
|0|0|0|G|
  • A = 2.0 / (right - left)
  • B = 2.0 / (top - bottom)
  • C = -2.0 / (far - near)
  • D = -(right + left) / (right - left)
  • E = -(top + bottom) / (top - bottom)
  • F = -(far + near) / (far - near)
  • G = 1.0

Vertex Shader

/// <summary>
/// A GLSL vertex shader to convert 3d Tile GLTF geospatial vertices to cartographic coordinates 
///
/// GLTF Transforms
///     https://github.com/CesiumGS/3d-tiles/blob/main/specification/README.md#transforms
/// ECEF to LLA 
///     https://microem.ru/files/2012/08/GPS.G1-X-00006.pdf
/// </summary>
#version 330
// Inputs
in vec3    in_vertex;
in vec3    in_normal;
in vec2    in_texcoord;
// Uniform Inputs
uniform mat4 gltfMatrix; // Converts transformed ECEF coords 
uniform mat4 mvpMatrix;  // Converts LLA coords in degrees / meters to clip space coordinates
// Constants
const float a = 6378137;              // WGS 84 ellipsoid semi-major axis
const float b = 6356752.314;          // WGS 84 ellipsoid semi-minor axis
const float ecc1Sqrd = 0.00669438;    // WGS 84 ellipsoid first eccentricity squared 
const float ecc2Sqrd = 0.006739497;   // WGS 84 ellipsoid second eccentricity squared 
const float RAD2DEG = 57.2957795132;  // Radians to degrees conversion
// Outputs
out vec3      normal;
out vec2      texcoord;
out vec3      lonLatAlt;
// Source
void main(void){
   normal = normalize(in_normal);
   texcoord = in_texcoord;
   
   // Apply GLTF transform on GLTF XYZ to ECEF XYZ 
   vec4 transformedVertex = gltfMatrix * vec4(in_vertex, 1.0); 
   
   float X = transformedVertex.x;
   float Y = transformedVertex.y;
   float Z = transformedVertex.z;

   // Compute Longitude in radians
   lonLatAlt.x = atan(Y, X);    

   float p = sqrt(X*X + Y*Y);
   float theta = atan(Z*a,p*b);
   float sinTheta = sin(theta);
   float cosTheta = cos(theta);
   float latY = Z + ecc2Sqrd * b * sinTheta * sinTheta * sinTheta;
   float latX = p - ecc1Sqrd * a * cosTheta * cosTheta * cosTheta;

   //Compute latitude in radians
   lonLatAlt.y = atan(latY,latX);       
   
   float sinLat = sin(lonLatAlt.y);
   float N = a / sqrt(1.0 - ecc1Sqrd * sinLat * sinLat);

   // Compute elevation
   lonLatAlt.z = (p / cos(lonLatAlt.y)) - N;
   
   // Convert to degrees
   lonLatAlt *= vec3(RAD2DEG, RAD2DEG, 1.0);

   //Convert to clip space 
   gl_Position = mvpMatrix * vec4(lonLatAlt, 1.0);  
   gl_Position.z = 0;  // Prevents any altitude clipping due to numerical issues

Fragment Shader

The fragment shader mainly exists to just write the values to the output array bound to the framebuffer. But when the gl_FragCoord.z value is 0 it is part of the background and should be set for the NODATA value.

/// <summary>
/// A GLSL fragment shader to write elevation data 
/// </summary>
#version 330
// Inputs for vertex shader
uniform float noDataValue; // The NODATA value to assign to background pixels
in vec4         lonLatAlt; // The input LLA value from the vertex shader
void main(void){
   if(gl_FragCoord.z == 0){gl_FragColor.r = -noDataValue;}
   else {
        gl_FragColor.r = lonLatAlt.z;
   }
}

Raster Parameters

Raster dimensions are determined by the size of the framebuffer. The pixel sizes for raster file generation are computed by dividing the geospatial bounds by the raster dimensions. The ESRI Ascii Format is useful as a simple method to generate a raster file although replace the CELLSIZE with DX and DY to allow for non-square pixel sizes (as shown in the example header below).

NCOLS 2000
NROWS 2000
XLLCORNER 37.29720442
YLLCORNER 32.72223247
DX 0.00134201512
DY 0.00170822024
NODATA_VALUE -9999
⚠️ **GitHub.com Fallback** ⚠️