C++ · Geospatial Metrics

Geospatial Metrics

#include <numkong/geospatial.hpp> 8.0.0-rc.1 · b1f866ab

The native SDK exposes both the fast spherical approximation and the more accurate ellipsoidal solver, each between every row and every column of two packed point sets. Inputs are in radians. Outputs are in meters.

// Statue of Liberty (40.6892°N, 74.0445°W) and Big Ben (51.5007°N, 0.1246°W)
nk_f64_t latitudes[] = {0.7101605100, 0.8988567821}, longitudes[] = {-1.2923203180, -0.0021746802};

// A pack is built once per point set, on an `nk_default_alignment_k` boundary
nk_size_t bytes = 0;
nk_vincentys_pack_size_f64_best(2, capabilities, &bytes);
void *vincentys_pack = aligned_alloc(nk_default_alignment_k, bytes);
nk_vincentys_pack_f64_best(latitudes, longitudes, 2, vincentys_pack, capabilities, NULL);

// Every row against every column, so the diagonal is zero
nk_f64_t distances[2][2];
nk_vincentys_packed_f64_best(vincentys_pack, vincentys_pack, &distances[0][0], 2 * sizeof(nk_f64_t), capabilities,
                             NULL); // distances[0][1] ≈ 5,589,857 m on WGS-84

nk_haversines_pack_size_f64_best(2, capabilities, &bytes);
void *haversines_pack = aligned_alloc(nk_default_alignment_k, bytes);
nk_haversines_pack_f64_best(latitudes, longitudes, 2, haversines_pack, capabilities, NULL);
nk_haversines_packed_f64_best(haversines_pack, haversines_pack, &distances[0][0], 2 * sizeof(nk_f64_t),
                              capabilities, NULL); // distances[0][1] ≈ 5,574,848 m on the mean sphere, 15 km less

API

Functions

haversinesGreat-circle metres between every row and every column point on the mean sphere, through the haversine 2R · atan2(√h, √(1 − h)) evaluated in precision_type_, where h = sin²(Δφ / 2) + cos φ₁ · cos φ₂ · sin²(Δλ / 2).
geospatial.hpp:39 ↗
template <numeric_dtype precision_type_>
void haversines(
    f64_t const *a_latitudes,
    f64_t const *a_longitudes,
    std::size_t rows,
    f64_t const *b_latitudes,
    f64_t const *b_longitudes,
    std::size_t columns,
    f64_t *results,
    std::size_t results_stride) noexcept;

Parameters

a_latitudesinThe rows row latitudes, in radians.
a_longitudesinThe rows row longitudes, in radians.
rowsinNumber of row points.
b_latitudesinThe columns column latitudes, in radians.
b_longitudesinThe columns column longitudes, in radians.
columnsinNumber of column points.
resultsoutThe [rows,columns] distances, results_stride bytes apart.
results_strideinBytes between the rows of results.
vincentysGeodesic metres on the WGS-84 ellipsoid between every row and every column point, through Vincenty's inverse formula in precision_type_, with the reduced latitudes taken through atan2 to stay finite at the poles; parameters as haversines.
geospatial.hpp:65 ↗
template <numeric_dtype precision_type_>
void vincentys(
    f64_t const *a_latitudes,
    f64_t const *a_longitudes,
    std::size_t rows,
    f64_t const *b_latitudes,
    f64_t const *b_longitudes,
    std::size_t columns,
    f64_t *results,
    std::size_t results_stride) noexcept;
sqeuclideans3dSquared Euclidean distances between every row point and every column point in 3D.
geospatial.hpp:173 ↗
template <numeric_dtype in_type_, numeric_dtype result_type_ = std::conditional_t<std::is_same_v<in_type_, i32_t>, u64_t, f64_t>>
requires (std::is_same_v<in_type_, i32_t> || std::is_same_v<in_type_, f64_t>)
status_t sqeuclideans3d(
    in_type_ const *a,
    std::size_t rows,
    in_type_ const *b,
    std::size_t columns,
    result_type_ *results,
    std::size_t results_stride,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;

Parameters

ainThe rows points, 3 interleaved coordinates each.
rowsinNumber of row points.
binThe columns points, 3 interleaved coordinates each.
columnsinNumber of column points.
resultsoutThe [rows,columns] squares, results_stride bytes apart.
results_strideinBytes between the rows of results.
capabilitiesinCapabilities to pick from, or zero for the reference in result_type_.
streaminNull on the CPU, or the stream of the device capabilities describes.
haversines_pack_sizeReturns the bytes a Haversine pack of points takes.
geospatial.hpp:213 ↗
expected<std::size_t> haversines_pack_size(
    std::size_t points,
    nk_capability_t capabilities = default_capabilities()) noexcept;

Parameters

pointsinNumber of points in the pack.
capabilitiesinCapabilities to pick from, or zero for the serial layout.
haversines_packPacks points as unit vectors on the mean sphere, for haversines_packed.
geospatial.hpp:230 ↗
status_t haversines_pack(
    f64_t const *latitudes,
    f64_t const *longitudes,
    std::size_t count,
    void *packed,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;

Parameters

latitudesinThe count latitudes, in radians.
longitudesinThe count longitudes, in radians.
countinNumber of points.
packedoutAt least haversines_pack_size bytes, aligned to nk_default_alignment_k.
capabilitiesinCapabilities to pick from, or zero for the serial kernel.
streaminNull on the CPU, or the stream of the device capabilities describes.
haversines_packedGreat-circle metres between every row and every column point of two Haversine packs.
geospatial.hpp:249 ↗
status_t haversines_packed(
    void const *a_packed,
    void const *b_packed,
    f64_t *results,
    std::size_t results_stride,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;

Parameters

a_packedinThe row points' Haversine pack.
b_packedinThe column points' Haversine pack.
resultsoutThe [rows,columns] distances, results_stride bytes apart.
results_strideinBytes between the rows of results.
capabilitiesinCapabilities to pick from, or zero for the serial kernel.
streaminNull on the CPU, or the stream of the device capabilities describes.
vincentys_pack_sizeReturns the bytes a Vincenty pack of points takes, as haversines_pack_size does.
geospatial.hpp:260 ↗
expected<std::size_t> vincentys_pack_size(
    std::size_t points,
    nk_capability_t capabilities = default_capabilities()) noexcept;
vincentys_packPacks points as their longitudes and the sines and cosines of their reduced latitudes, for vincentys_packed, taking what haversines_pack takes.
geospatial.hpp:270 ↗
status_t vincentys_pack(
    f64_t const *latitudes,
    f64_t const *longitudes,
    std::size_t count,
    void *packed,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;
vincentys_packedGeodesic metres on the WGS-84 ellipsoid between every row and every column point of two Vincenty packs, taking what haversines_packed takes.
geospatial.hpp:282 ↗
status_t vincentys_packed(
    void const *a_packed,
    void const *b_packed,
    f64_t *results,
    std::size_t results_stride,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;
terrain_pack_sizeReturns the bytes the terrain pack of a rows × columns grid takes.
geospatial.hpp:298 ↗
expected<std::size_t> terrain_pack_size(
    std::size_t rows,
    std::size_t columns,
    nk_capability_t capabilities = default_capabilities()) noexcept;

Parameters

rowsinNumber of grid rows.
columnsinNumber of grid columns.
capabilitiesinCapabilities to pick from, or zero for the serial layout.
terrain_packPacks the maxima of every 2ᵏ × 2ᵏ block of a north-up grid, its shape and scale.
geospatial.hpp:318 ↗
status_t terrain_pack(
    f32_t const *heights,
    std::size_t rows,
    std::size_t columns,
    std::size_t heights_stride,
    f64_t row_metres,
    f64_t column_metres,
    void *packed,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;

Parameters

heightsinThe [rows,columns] heights in metres, heights_stride bytes apart.
rowsinNumber of grid rows.
columnsinNumber of grid columns.
heights_strideinBytes between the rows of heights.
row_metresinMetres between adjacent rows.
column_metresinMetres between adjacent columns.
packedoutAt least terrain_pack_size bytes, aligned to nk_default_alignment_k.
capabilitiesinCapabilities to pick from, or zero for the serial kernel.
streaminNull on the CPU, or the stream of the device capabilities describes.
terrain_heights_packedBilinear terrain heights under count interleaved [column, row] points, NaN off the grid, dispatched as terrain_pack is.
geospatial.hpp:332 ↗
status_t terrain_heights_packed(
    f32_t const *heights,
    void const *packed,
    f32_t const *points,
    std::size_t count,
    f32_t *results,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;
clearances_packedClearances in metres from every row viewpoint to every column target over one terrain, as nk_clearances_packed_f32_best defines them, dispatched as terrain_pack is.
geospatial.hpp:344 ↗
status_t clearances_packed(
    f32_t const *heights,
    void const *packed,
    f32_t const *viewpoints,
    f32_t const *views,
    std::size_t rows,
    f32_t const *targets,
    std::size_t columns,
    f64_t refraction,
    f32_t *results,
    std::size_t results_stride,
    nk_capability_t capabilities = default_capabilities(),
    nk_stream_t stream = nullptr) noexcept;