Repository navigation
Geodesy and Tiling
This page covers vgl's coordinate systems and the functions that convert between them: geodetic latitude, longitude and height, Earth-centred Earth-fixed (ECEF) coordinates, and local East-North-Up (ENU) frames, together with the ellipsoid constants, angle helpers, the Earth rotation angle and heights above mean sea level from a geoid pack. It then summarizes the tiling scheme that packs use to address tiles, and ends with worked examples: placing a host object on the globe and building a camera that looks at a latitude and longitude.
The API is in include/vgl/geodesy.h and the C++ wrapper in include/vgl/geodesy.hpp. The geodesy functions take no context, allocate nothing and are free-threaded, so you can call them from any thread at any time. They are tested against GeographicLib reference values.
| System | Type | Definition |
|---|---|---|
| Geodetic | vgl_geodetic |
latitude (geodetic, in radians, in [-pi/2, pi/2]), longitude (radians, any finite value on input, [-pi, pi] on output) and height (metres above the ellipsoid). |
| ECEF | vgl_dvec3 |
Metres. Right-handed, origin at the centre of the ellipsoid, +Z through the north pole, +X through latitude 0, longitude 0, so +Y through latitude 0, longitude 90 degrees east. |
| ENU |
vgl_dvec3 in a frame |
Metres. Origin at a chosen geodetic position, +X east, +Y north, +Z up along the ellipsoid normal. |
All angles in the API are radians and all lengths metres (C API conventions, section 4.3). Everything vgl hands you in world space, from view matrices to ray hits and mesh centres, is ECEF in double precision. Matrices are vgl_dmat4: column-major, multiplying column vectors, element (row r, column c) at m[c * 4 + r], translation in m[12], m[13], m[14]. A Direct3D-style host with row-major matrices that multiply row vectors can pass its matrices unchanged, because the two differences cancel out (see Core concepts).
Heights are above the ellipsoid everywhere in the geodesy API. Heights above mean sea level (what maps, flight simulators and most people mean by altitude) differ from them by the geoid height, from about -107 m to +86 m; see Heights above mean sea level. GPS receivers report ellipsoidal heights natively.
typedef struct vgl_ellipsoid
{
double semi_major_axis; /* equatorial radius in metres */
double flattening; /* (a - b) / a */
} vgl_ellipsoid;
#define VGL_WGS84_SEMI_MAJOR_AXIS 6378137.0
#define VGL_WGS84_FLATTENING (1.0 / 298.257223563)
vgl_result vgl_ellipsoid_wgs84_get(vgl_ellipsoid* out_ellipsoid);A valid ellipsoid has a finite semi_major_axis greater than 0 and a flattening in [0, 0.5); a flattening of 0 is a sphere. Every conversion function takes the ellipsoid explicitly, so there is no hidden global. A context uses WGS84 unless vgl_context_create_info.ellipsoid names another, and its packs must have been baked on the same one. Unless you render another planet, use WGS84.
In C++, vgl::ellipsoid_wgs84_get() returns the ellipsoid directly.
#define VGL_PI 3.14159265358979323846
#define VGL_DEGREES_TO_RADIANS(degrees) ((degrees) * (VGL_PI / 180.0))
#define VGL_RADIANS_TO_DEGREES(radians) ((radians) * (180.0 / VGL_PI))
vgl_result vgl_degrees_to_radians(double degrees, double* out_radians);
vgl_result vgl_radians_to_degrees(double radians, double* out_degrees);The functions do the same as the macros, for language bindings that cannot use macros; NaN and infinities pass through. In C++, vgl::pi, vgl::degrees_to_radians and vgl::radians_to_degrees are constexpr.
vgl_result vgl_ellipsoid_geodetic_to_ecef(const vgl_ellipsoid* ellipsoid,
const vgl_geodetic* geodetic, vgl_dvec3* out_ecef);
vgl_result vgl_ellipsoid_ecef_to_geodetic(const vgl_ellipsoid* ellipsoid,
const vgl_dvec3* ecef, vgl_geodetic* out_geodetic);-
vgl_ellipsoid_geodetic_to_ecefreturnsVGL_RESULT_INVALID_ARGUMENTfor a null pointer, an invalid ellipsoid, a value that is not finite or a latitude out of range. A latitude up to 1e-12 radians beyond a pole is accepted and treated as the pole, so thatVGL_DEGREES_TO_RADIANS(90.0)is always valid despite rounding. -
vgl_ellipsoid_ecef_to_geodeticis closed-form and accurate to a few nanometres for any position, including positions deep inside the ellipsoid. On the polar axis the longitude is 0; at the very centre, the result is the north pole at a height of minus the semi-minor axis.
vgl_ellipsoid wgs84;
vgl_geodetic utrecht;
vgl_dvec3 ecef;
vgl_geodetic back;
vgl_ellipsoid_wgs84_get(&wgs84);
utrecht.latitude = VGL_DEGREES_TO_RADIANS(52.1);
utrecht.longitude = VGL_DEGREES_TO_RADIANS(5.1);
utrecht.height = 10.0;
if (vgl_ellipsoid_geodetic_to_ecef(&wgs84, &utrecht, &ecef) == VGL_RESULT_SUCCESS)
{
/* ecef.x is about 3910649.03 m */
vgl_ellipsoid_ecef_to_geodetic(&wgs84, &ecef, &back); /* back.height is 10.0 */
}const vgl_ellipsoid wgs84 = vgl::ellipsoid_wgs84_get();
const vgl_geodetic utrecht{vgl::degrees_to_radians(52.1), vgl::degrees_to_radians(5.1), 10.0};
const vgl::expected<vgl_dvec3> ecef = vgl::ellipsoid_geodetic_to_ecef(wgs84, utrecht);
const vgl::expected<vgl_geodetic> back = vgl::ellipsoid_ecef_to_geodetic(wgs84, *ecef);vgl_result vgl_ellipsoid_enu_frame_get(const vgl_ellipsoid* ellipsoid, const vgl_geodetic* origin,
vgl_dmat4* out_enu_to_ecef, vgl_dmat4* out_ecef_to_enu);An ENU frame is the local tangent frame at a position: +X east, +Y north, +Z up along the ellipsoid normal, with its origin at the position (including its height). Both matrices are rigid transforms:
- In
out_enu_to_ecef, columns 0, 1 and 2 are the east, north and up axes in ECEF, and column 3 (m[12..14]) is the origin in ECEF. It takes a point in local metres to ECEF: it is the model matrix of an object standing at the origin. -
out_ecef_to_enuis its inverse: it takes an ECEF point to local east, north and up offsets, for example to tell how far north of a reference point something is.
Either output may be null, but not both. In C++, vgl::ellipsoid_enu_frame_get(ellipsoid, origin) returns an expected<vgl::enu_frame> holding enu_to_ecef and ecef_to_enu.
vgl_result vgl_earth_rotation_angle_get(double julian_date, double* out_angle);Returns the angle, in [0, 2 pi), by which the Earth has turned from the celestial frame (the equator and equinox of J2000) at a time given as a Julian date in UT1 (UTC is within a second of it; 2451545.0 is 1 January 2000 at 12:00). A direction d in the celestial frame is, in ECEF:
(d.x cos(a) + d.y sin(a), -d.x sin(a) + d.y cos(a), d.z)
This is the IAU 2000 definition, which leaves out precession and nutation (a fraction of a degree over decades): enough to place stars behind the globe. The backends take it in their draw parameters (earth_rotation_angle) to orient the star field; see Atmosphere, stars and clouds. From a Unix time in seconds, julian_date = unix_seconds / 86400.0 + 2440587.5.
const double julian_date = static_cast<double>(std::time(nullptr)) / 86400.0 + 2440587.5;
const double angle = vgl::earth_rotation_angle_get(julian_date).value();vgl reads the geoid from a geoid pack (an EGM2008 pack is in the sample set, samples/data/packs/egm2008), added to a context with vgl_pack_create like any other pack; a context takes at most one. The geoid height N is the height of the geoid above the ellipsoid, and:
height above mean sea level = ellipsoidal height - N
There is no function that returns N alone. Instead:
-
vgl_context_height_querywithVGL_HEIGHT_REFERENCE_MEAN_SEA_LEVELreturns the terrain height above mean sea level, and withVGL_HEIGHT_REFERENCE_ELLIPSOIDabove the ellipsoid. The difference of the two at the same position is N there (the headless sample checks exactly that). Both need an elevation pack with data at the position; without a geoid pack the mean-sea-level query returnsVGL_RESULT_INVALID_OPERATION, and while the geoid is loading it returnsVGL_RESULT_NOT_READYwith a height above the ellipsoid. - For rendering, each draw item of the render interface carries the geoid height at its tile's four corners (
geoid_heights,has_geoid), which consumers interpolate for a sea surface or for heights above mean sea level in shaders.
The geoid is interpolated bilinearly from the pack's deepest level, as section 6.4 of the pack format specification specifies. To convert an altitude above mean sea level, such as an aircraft's, to the ellipsoidal height the geodesy functions need, add N: h_ellipsoid = h_msl + N.
// N at a position, from two height queries (radians).
std::optional<double> geoid_height(vgl::context& ctx, double lat, double lon)
{
const auto e = ctx.height_query(lat, lon, vgl::height_reference::ellipsoid);
const auto m = ctx.height_query(lat, lon, vgl::height_reference::mean_sea_level);
if (!e || !m || e->status != vgl::result::success || m->status != vgl::result::success)
{
return std::nullopt; // not loaded yet (update and ask again), or no data here
}
return e->height - m->height;
}Two functions in vgl/context.h complete the set, for labels and markers: vgl_view_project takes an ECEF position and vgl_view_geodetic_project a geodetic one (with an optional ellipsoid, null for WGS84), and both return the position in the viewport of a vgl_view: x and y in pixels from the top left, y downwards, and z the clip-space w (the distance along the view axis for a perspective projection). They return VGL_RESULT_NOT_FOUND for a position behind the camera, and do not tell whether the position is hidden by terrain; vgl_context_ray_query from the camera does that. See Context API.
Packs divide the globe into tiles with a fixed scheme, specified in section 2 of the pack format specification. You do not need it to use vgl: the core selects, reads and decodes tiles by itself, and there is no public tiling API. The scheme is implemented once in src/tiling and shared by the runtime and the Baker. It is summarized here because it explains pack layouts, file names, what vgl-bake options such as the levels mean, and what a host sees in its I/O reads.
The globe is covered by a geographic (equirectangular) quadtree on longitude and geodetic latitude. A tile is addressed by its level L, column x and row y:
- Level 0 has two tiles of 180 by 180 degrees: (0, 0, 0) is the western hemisphere and (0, 1, 0) the eastern.
- Level L has 2L+1 columns and 2L rows, and each tile spans s = 180 / 2L degrees in both longitude and latitude.
- Columns count eastwards from the antimeridian, rows southwards from the north pole, so tile rows run the same way as image rows.
- Tile (L, x, y) covers longitudes [-180 + x s, -180 + (x + 1) s] and latitudes [90 - (y + 1) s, 90 - y s].
- The children of (L, x, y) are (L + 1, 2x + i, 2y + j), numbered i + 2j: 0 north-west, 1 north-east, 2 south-west, 3 south-east.
- Valid levels are 0 to 28; a pack uses any contiguous range of them.
A position on a shared edge belongs to the tile east of a meridian and south of a parallel (tiles own their west and north edges), except at the east and south limits of the grid.
Within a pack, tiles are sorted by a 64-bit key, level first, then along a Morton (Z-order) curve:
key = (L << 58) | morton(x, y) // bit i of x at bit 2i, bit i of y at bit 2i + 1
Siblings are adjacent in key order and every subtree is a contiguous key range. For example, tile (7, 181, 37) has key 0x1C00000000004D33.
-
Elevation tiles are grid-registered: a tile of
tile_sizeN (2k + 1, by default 129) has N by N samples on and between its edges, and neighbouring tiles share their border samples exactly. Sample (i, j) lies at longitude west + i s / (N - 1) and latitude north - j s / (N - 1). -
Imagery tiles are cell-registered: a tile of
tile_sizeP (a power of two, by default 256) has P by P pixels covering it exactly, plus aborderof 0 or 1 pixel (by default 1) repeated from its neighbours, so that textures filter across tile edges without seams.
A pack groups its tiles into files at its shard_level S (by default max_level - 5):
<pack>/
pack.vglm manifest
base-<digest>.vgls every tile above level S
shards/<y>/<x>-<digest>.vgls the subtree of tile (S, x, y)
A tile (L, x, y) with L >= S lives in the shard (x >> (L - S), y >> (L - S)). The digest in each name changes whenever the shard's content does, which is what lets the disk cache and web caches keep shards forever (see Pack format overview and Hosting packs).
For 129-sample elevation tiles at the equator on WGS84 (pack format, section 2.5):
| Level | Tile span | Sample spacing | Used by default for |
|---|---|---|---|
| 0 | 180 degrees | 156.5 km | |
| 5 | 5.625 degrees | about 4.9 km | The EGM2008 geoid pack (max_level 5) |
| 7 | 1.40625 degrees | about 1.2 km | Blue Marble imagery (max_level 7, 256-pixel tiles, about 610 m per pixel) |
| 8 | 0.703125 degrees | 611.5 m | The GEBCO global base (max_level 8) |
| 10 | 0.17578125 degrees | 152.9 m | |
| 12 | 0.0439453125 degrees | 38.2 m | Copernicus GLO-30 (max_level 12) |
| 13 | 0.02197265625 degrees | 19.1 m |
This host-side code mirrors the scheme, for example to predict which shard files a region needs before you copy or prefetch them. It is not part of vgl's API.
#include <cmath>
#include <cstdint>
struct tile { std::uint32_t level, x, y; };
// The tile at a level that holds a position in degrees.
tile tile_at(double latitude, double longitude, std::uint32_t level)
{
const double s = 180.0 / static_cast<double>(1u << level);
double lon = std::fmod(longitude + 180.0, 360.0);
if (lon < 0.0)
{
lon += 360.0;
}
const double columns = static_cast<double>(2u << level);
const double rows = static_cast<double>(1u << level);
const double lat = std::fmin(std::fmax(latitude, -90.0), 90.0);
return {level,
static_cast<std::uint32_t>(std::fmin(std::floor(lon / s), columns - 1.0)),
static_cast<std::uint32_t>(std::fmin(std::floor((90.0 - lat) / s), rows - 1.0))};
}
std::uint64_t tile_key(const tile& t)
{
std::uint64_t morton = 0;
for (std::uint32_t i = 0; i < 32; ++i)
{
morton |= static_cast<std::uint64_t>((t.x >> i) & 1u) << (2 * i);
morton |= static_cast<std::uint64_t>((t.y >> i) & 1u) << (2 * i + 1);
}
return (static_cast<std::uint64_t>(t.level) << 58) | morton;
}
// Galdhopiggen (61.6363 N, 8.3125 E) in a GLO-30 pack with max_level 12 and shard_level 7:
// tile_at(61.6363, 8.3125, 12) is (12, 4285, 645), and its shard is the subtree of
// (7, 4285 >> 5, 645 >> 5) = (7, 133, 20): the file shards/20/133-<digest>.vgls.To put an object (a building, a marker, an aircraft) on the terrain, build its model matrix from an ENU frame at its position, with the height taken from a height query:
#include <vgl/vgl.hpp>
// c = a * b for column-major matrices that multiply column vectors.
vgl_dmat4 multiply(const vgl_dmat4& a, const vgl_dmat4& b)
{
vgl_dmat4 c{};
for (int col = 0; col < 4; ++col)
{
for (int row = 0; row < 4; ++row)
{
double sum = 0.0;
for (int k = 0; k < 4; ++k)
{
sum += a.m[k * 4 + row] * b.m[col * 4 + k];
}
c.m[col * 4 + row] = sum;
}
}
return c;
}
// The model matrix of an object standing on the terrain at a latitude and longitude in
// degrees, turned to a heading (degrees clockwise from north). The object's local axes are
// x east, y north and z up before the heading is applied. Returns false while the terrain
// height is not known yet; call ctx.update() and try again.
bool model_on_terrain(vgl::context& ctx, double lat_deg, double lon_deg, double heading_deg,
vgl_dmat4& out_model)
{
const double lat = vgl::degrees_to_radians(lat_deg);
const double lon = vgl::degrees_to_radians(lon_deg);
const vgl::expected<vgl::height_sample> ground = ctx.height_query(lat, lon);
if (!ground || std::isnan(ground->height))
{
return false; // not_found: no elevation data here; NaN: nothing loaded yet
}
// ground->status may be not_ready: a coarser height for now, finer after later updates.
const vgl_ellipsoid wgs84 = vgl::ellipsoid_wgs84_get();
const vgl::expected<vgl::enu_frame> frame =
vgl::ellipsoid_enu_frame_get(wgs84, {lat, lon, ground->height});
if (!frame)
{
return false;
}
// Turn about the local up axis: heading h clockwise from north.
const double h = vgl::degrees_to_radians(heading_deg);
vgl_dmat4 turn{};
turn.m[0] = std::cos(h); turn.m[1] = -std::sin(h); // local x (right of the object)
turn.m[4] = std::sin(h); turn.m[5] = std::cos(h); // local y (the object's forward)
turn.m[10] = 1.0;
turn.m[15] = 1.0;
out_model = multiply(frame->enu_to_ecef, turn);
return true;
}Precision. ECEF coordinates are millions of metres, so single-precision floats are about half a metre coarse at the Earth's surface, which shows as jitter. Compose the model-view matrix in double precision, multiply(view_matrix, model), and convert it to float only at the end, as vgl does for its own tiles (render interface, "no jitter"). Keep the object's vertices in local metres.
Depth. For the object to depth-test against the globe, draw it with the same depth convention as the backend (logarithmic by default). The OpenGL backend and Direct3D 11 backend pages give the matching shader snippet; see also Core concepts.
Altitude above sea level. For an aircraft whose altitude is above mean sea level, use h_ellipsoid = h_msl + N (see Heights above mean sea level) as the frame's height instead of the terrain height.
An "orbit" camera looks at a target on the globe from a distance, a heading (clockwise from north) and a pitch (down from the horizon, 90 degrees straight down). The ENU frame at the target gives the directions, and a look-at matrix turns them into a view matrix for vgl_view. The samples use the same construction in samples/common/sample_camera.hpp.
#include <vgl/vgl.hpp>
#include <cmath>
vgl_dvec3 add(vgl_dvec3 a, vgl_dvec3 b) { return {a.x + b.x, a.y + b.y, a.z + b.z}; }
vgl_dvec3 scale(vgl_dvec3 a, double s) { return {a.x * s, a.y * s, a.z * s}; }
double dot(vgl_dvec3 a, vgl_dvec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
vgl_dvec3 cross(vgl_dvec3 a, vgl_dvec3 b)
{
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
vgl_dvec3 normalize(vgl_dvec3 v) { return scale(v, 1.0 / std::sqrt(dot(v, v))); }
// A right-handed view matrix (the camera looks down its -Z axis), as OpenGL's lookAt.
vgl_dmat4 look_at(vgl_dvec3 eye, vgl_dvec3 target, vgl_dvec3 up)
{
const vgl_dvec3 f = normalize(add(target, scale(eye, -1.0)));
const vgl_dvec3 s = normalize(cross(f, up));
const vgl_dvec3 u = cross(s, f);
vgl_dmat4 m{};
m.m[0] = s.x; m.m[4] = s.y; m.m[8] = s.z;
m.m[1] = u.x; m.m[5] = u.y; m.m[9] = u.z;
m.m[2] = -f.x; m.m[6] = -f.y; m.m[10] = -f.z;
m.m[12] = -dot(s, eye);
m.m[13] = -dot(u, eye);
m.m[14] = dot(f, eye);
m.m[15] = 1.0;
return m;
}
// Looks at (lat, lon, height) in degrees and metres above the ellipsoid, from range metres
// away, turned to heading and tilted down by pitch, both in degrees.
vgl_dmat4 orbit_view(double lat_deg, double lon_deg, double height, double range,
double heading_deg, double pitch_deg)
{
const vgl_ellipsoid wgs84 = vgl::ellipsoid_wgs84_get();
const vgl_geodetic target{vgl::degrees_to_radians(lat_deg), vgl::degrees_to_radians(lon_deg),
height};
const vgl_dmat4 enu = vgl::ellipsoid_enu_frame_get(wgs84, target)->enu_to_ecef;
const vgl_dvec3 east{enu.m[0], enu.m[1], enu.m[2]};
const vgl_dvec3 north{enu.m[4], enu.m[5], enu.m[6]};
const vgl_dvec3 up{enu.m[8], enu.m[9], enu.m[10]};
const vgl_dvec3 centre{enu.m[12], enu.m[13], enu.m[14]};
const double h = vgl::degrees_to_radians(heading_deg);
const double p = vgl::degrees_to_radians(pitch_deg);
const vgl_dvec3 ahead = add(scale(north, std::cos(h)), scale(east, std::sin(h)));
const vgl_dvec3 forward = add(scale(ahead, std::cos(p)), scale(up, -std::sin(p)));
const vgl_dvec3 eye = add(centre, scale(forward, -range));
// Looking straight down, up is parallel to the view direction: use ahead instead.
return look_at(eye, centre, pitch_deg > 89.0 ? ahead : up);
}
// Usage: Galdhopiggen from 15 km, looking south-east, 30 degrees down.
auto view = vgl::init<vgl_view>();
view.view_matrix = orbit_view(61.6363, 8.3125, 2469.0, 15000.0, 135.0, 30.0);
view.projection_matrix = projection; // e.g. an OpenGL perspective; see Core concepts
view.viewport_width = width;
view.viewport_height = height;
static_cast<void>(ctx->view_set(&view));For the camera's own position as a latitude, longitude and altitude (a "fly-over" camera, as the samples' sample::camera does), use vgl::ellipsoid_geodetic_to_ecef for the eye and the ENU axes at the eye for the directions. To keep the camera above the terrain, query the height under it, or cast vgl_context_ray_query downwards, and clamp. For the projection, near and far planes, and depth conventions, see Core concepts and the backend pages.
-
Degrees passed as radians. Every angle in the API is in radians; convert with
VGL_DEGREES_TO_RADIANSorvgl::degrees_to_radians. A latitude of 52 "radians" is rejected as out of range, but a longitude is accepted silently. - Ellipsoidal versus mean-sea-level heights. Mixing them puts objects tens of metres too high or low. Elevation sources referenced to EGM2008 are converted to ellipsoidal heights by the baker.
- Single-precision world positions. Keep ECEF in doubles until the final model-view matrix.
- A different ellipsoid in the context than in the packs. Packs must be baked on the context's ellipsoid; for Earth, keep both at WGS84.
- Expecting a tiling API. The tile helpers are internal; the sketch above is all a host normally needs.
- Core concepts: coordinates, matrices and depth conventions in the bigger picture.
- Context API: height and ray queries, views and projection helpers.
- Pack format overview and the pack format specification.
- Baker: choosing levels and shard levels when baking.
Getting vgl
Using the library
Rendering
Data
Reference
Specifications