# Geospatial Metrics

> NumKong for Python

Geospatial kernels score every row point against every column point into a `[rows, columns]` matrix.
Inputs are in radians, and outputs are in metres.
Haversine and Vincenty take both sides as packs, built once per point set, so the many-to-many kernels evaluate no trigonometry before their per-pair work.

```python
import numpy as np
import numkong as nk

# Statue of Liberty, Big Ben, Eiffel Tower
latitudes = np.radians([40.6892, 51.5007, 48.8584])
longitudes = np.radians([-74.0445, -0.1246, 2.2945])

ellipsoid = nk.vincentys_pack(latitudes, longitudes)  # any input dtype is cast to float64
sphere = nk.haversines_pack(latitudes, longitudes)
assert len(sphere) == sphere.points == 3 and sphere.shape == (3, 2)

vincentys = np.asarray(nk.vincentys_packed(ellipsoid, ellipsoid))  # Liberty → Big Ben ≈ 5,589,857 m on WGS-84
haversines = np.asarray(nk.haversines_packed(sphere, sphere))  # Liberty → Big Ben ≈ 5,574,848 m on the mean sphere
assert vincentys.shape == haversines.shape == (3, 3)

# Earth-fixed centimetres square exactly in 64 bits
stations = np.array([[631_000_000, 0, 0], [0, 637_800_000, 0]], dtype=np.int32)
squares = np.asarray(nk.sqeuclideans3d(stations, stations))
assert squares.dtype == np.uint64 and squares[0, 1] == 631_000_000**2 + 637_800_000**2
```

`nk.haversines_packed` measures 2R · atan2(‖a − b‖, ‖a + b‖) with the IUGG mean radius of 6,371,008.8 m.
That form stays conditioned from a millimetre to the antipode, and reads within 0.5% of the ellipsoidal geodesic.
`nk.vincentys_packed` iterates Vincenty's inverse formula on the WGS-84 ellipsoid until the longitude settles within 1e-12 radians.
A pair near the antipode that never settles reports its last iterate after 100 iterations.
`nk.sqeuclideans3d` takes two `[n, 3]` matrices, both `int32` or both `float64`, and returns `uint64` or `float64` squares.
The `int32` squares are exact whenever every axis difference fits 32 bits, which Earth-fixed centimetres below 4,359 km of altitude always satisfy.
Its `threads=` argument splits the row points across workers.
`nk.HaversinesPackedMatrix.pack_size(points)` and `nk.VincentysPackedMatrix.pack_size(points)` return the bytes a pack takes.
A pack keeps its point count in a header ahead of its planar arrays, so a view of fewer points would need a header of its own.
Indexing therefore returns the pack itself for a key selecting every point in order, and raises `ValueError` for any other.

`import numkong`

`from numkong import HaversinesPackedMatrix, VincentysPackedMatrix`

## HaversinesPackedMatrix

Points packed by haversines_pack() as unit vectors on the mean sphere.

```python
class HaversinesPackedMatrix
```

The kernels read the point count from a header ahead of the planar arrays, so the only view is
the whole pack, and a key selecting fewer points raises ValueError.

## pack_*_size

Return the bytes a pack of the given number of points takes.

```python
@classmethod
def pack__size(points: int, /) -> int
```

## __len__*__len__

Return the number of packed points.

```python
def __len____len__() -> int
```

## __repr__*__repr__

Return a string representation.

```python
def __repr____repr__() -> str
```

## points

Number of packed points.

```python
points: int
```

## shape

Dimensions (points, 2) of the packed latitude and longitude pairs.

```python
shape: tuple[int, int]
```

## nbytes

Size of the pack in bytes.

```python
nbytes: int
```

## VincentysPackedMatrix

Points packed by vincentys_pack() as longitudes and reduced latitudes on the WGS-84 ellipsoid.

```python
class VincentysPackedMatrix
```

The kernels read the point count from a header ahead of the planar arrays, so the only view is
the whole pack, and a key selecting fewer points raises ValueError.

## points

Number of packed points.

```python
points: int
```

## shape

Dimensions (points, 2) of the packed latitude and longitude pairs.

```python
shape: tuple[int, int]
```

## nbytes

Size of the pack in bytes.

```python
nbytes: int
```

## sqeuclideans3d

Compute [rows,columns] squared distances of [rows,3] and [columns,3] points, both int32 or both float64.

```python
def sqeuclideans3d(
    a: _BufferType,
    b: _BufferType,
    /,
    *,
    threads: int = 1,
    capabilities: Capability | int | None = ...,
    stream: int | None = ...,
) -> Tensor
```

- `capabilities`: Capabilities to run within: every CPU one by default, or for packed operands the mask that packed them.
- `stream`: A GPU stream pointer, left None on the CPU.

## haversines_pack

Pack every point as a unit vector on the mean sphere, for haversines_packed().

```python
def haversines_pack(
    latitudes: _BufferType,
    longitudes: _BufferType,
    /,
    *,
    capabilities: Capability | int | None = ...,
    stream: int | None = ...,
) -> HaversinesPackedMatrix
```

- `capabilities`: Capabilities to run within: every CPU one by default, or for packed operands the mask that packed them.
- `stream`: A GPU stream pointer, left None on the CPU.

## vincentys_pack

Pack every point as its longitude and reduced latitude on the WGS-84 ellipsoid, for vincentys_packed().

```python
def vincentys_pack(
    latitudes: _BufferType,
    longitudes: _BufferType,
    /,
    *,
    capabilities: Capability | int | None = ...,
    stream: int | None = ...,
) -> VincentysPackedMatrix
```

- `capabilities`: Capabilities to run within: every CPU one by default, or for packed operands the mask that packed them.
- `stream`: A GPU stream pointer, left None on the CPU.

## haversines_packed

Compute great-circle distances on the sphere of the IUGG mean Earth radius, 6,371,008.8 m.

```python
def haversines_packed(
    a: HaversinesPackedMatrix,
    b: HaversinesPackedMatrix,
    /,
    *,
    capabilities: Capability | int | None = ...,
    stream: int | None = ...,
) -> Tensor
```

- `capabilities`: Capabilities to run within: every CPU one by default, or for packed operands the mask that packed them.
- `stream`: A GPU stream pointer, left None on the CPU.

## vincentys_packed

Compute geodesic distances on the WGS-84 ellipsoid with Vincenty's inverse formula.

```python
def vincentys_packed(
    a: VincentysPackedMatrix,
    b: VincentysPackedMatrix,
    /,
    *,
    capabilities: Capability | int | None = ...,
    stream: int | None = ...,
) -> Tensor
```

- `capabilities`: Capabilities to run within: every CPU one by default, or for packed operands the mask that packed them.
- `stream`: A GPU stream pointer, left None on the CPU.