Geometry

AbstractCartesianGeometry{T} <: AbstractGeometry{T}

Cartesian metrics. Default: CartesianGeometry.

AbstractEllipsoidalGeometry{T} <: AbstractGeometry{T}

Oblate-spheroid metrics. Default: SpheroidGeometry. Coordinate names match the spherical convention, (λ, φ[, h]), with φ the geodetic latitude.

A subtype implements semimajor_axis and flattening; the rest — semiminor_axis, eccentricity², the curvature radii, distance, area_element, scale_factors — follow from those two.

AbstractGeometry{T<:AbstractFloat}

Supertype for coordinate metrics (T is the float type).

AbstractSphericalGeometry{T} <: AbstractGeometry{T}

Spherical metrics. Default: SphericalGeometry.

A subtype implements radius; every other method here is written in terms of it.

CartesianGeometry()
CartesianGeometry(T)
CartesianGeometry{T}()

Flat metric in T, of any dimension.

It carries no grid spacing. A cell's extent belongs to the grid's axes, which state it per cell and per direction, and the dimension likewise lives on the grid.

PoleRotation(λp, φp)

The frame whose north pole sits at (λp, φp) of the original one — a rotated-pole grid's coordinate change. Apply it with rotate and undo it with unrotate.

The tilt's sine and cosine are part of the rotation: θ = φp − π/2 is a property of the frame, so it is resolved once here and read by every point carried through it.

rotate and unrotate work at the width of the points given them, carrying the rotation there; T is the width the frame itself is stored and carried from. similar_rotation moves a rotation between widths explicitly.

SphericalGeometry{T} <: AbstractSphericalGeometry{T}

Default spherical geometry with sphere radius R (meters; default Earth 6.371e6).

SpheroidGeometry(a, f)
SpheroidGeometry()

Oblate spheroid of equatorial radius a and flattening f = (a-b)/a. The no-argument form is WGS 84, a = 6378137.0, f = 1/298.257223563.

distance, area_element, volume_element and scale_factors differ from a sphere. Grid directions are (λ, φ, h), with h the height above the ellipsoid, where a sphere's third direction is an absolute radius.

Rectilinear grids and the index-space connectivity built on them work as they do for a sphere; the samplings are purely angular and so carry over unchanged. So does everything written against the two primitives embed and local_tangent_basis — the local frame, the vector and tensor rotations, the tangent-plane projection, and therefore the least-squares gradient and the scattered interpolation. See AbstractLonLatGeometry for why the frame is shared.

The spherical area routines do not: unstructured_grid's Voronoi areas and the cell areas behind spherical_grid are built on spherical excess and 4πR²/n, which are sphere identities. They stay restricted to AbstractSphericalGeometry, and raise on an ellipsoid.

AbstractLonLatGeometry{T}

The geometries whose coordinates are (λ, φ, …) about a polar axis: AbstractSphericalGeometry and AbstractEllipsoidalGeometry.

They share their local frame. On an oblate spheroid the geodetic latitude is defined by the surface normal, n̂ = (cosφ·cosλ, cosφ·sinλ, sinφ), and the parallel through a point is a circle in a plane of constant z, so ê_λ = (-sinλ, cosλ, 0) and ê_φ = n̂ × ê_λ — the sphere's three vectors at the same (λ, φ). Every frame rotation here is therefore one implementation for both hierarchies. What differs is the position a (λ, φ) sits at, which is embed's business, and the physical length of a unit coordinate step, which is scale_factors'.

The two hierarchies are siblings, so that an ellipsoid does not inherit the area identities — 4πR²/n, spherical excess — that hold on a sphere alone. This Union names the part they share.

Component names of an ambient-frame vector (as returned by local_tangent_basis).

Convert positional point values to the geometry's own float type T.

_enu_frame(geo, λ, φ) -> (ê_λ, ê_φ, ê_r)

The local east/north/up triad at (λ, φ) in ambient Cartesian components.

The one definition of the frame, so a vector rotation, a tensor rotation and a tangent-plane projection cannot drift apart. local_tangent_basis is its first two vectors under their coordinate names.

_lonlat(p) -> (λ, φ)
_xyz(v) -> (x, y, z)

A point's components for a consumer that needs a fixed number of them, with the arity checked.

λ, φ = as_ntuple(p) does not check it: on a one-component point the destructure reads p[2] and reports an index out of bounds, when what is wrong is the point. These say so instead, and return a pair or a triple whose type does not depend on how many components the point had.

The arity is the dispatch, so the accepted width is a method signature and the tuple is indexed only where its type is known wide enough.

area_element(geo::AbstractCartesianGeometry, dx, dy)
area_element(geo::AbstractSphericalGeometry, φ, dλ, dφ)

Local cell area from the cell's own extents: dx·dy on the plane, R²·cosφ·dλ·dφ on the sphere.

area_element(geo::AbstractEllipsoidalGeometry, φ, dλ, dφ)

M(φ)·N(φ)cosφ·dλ·dφ, the exact ellipsoidal surface element — equal to a²(1-e²)cosφ / (1-e²sin²φ)²·dλ·dφ.

as_ntuple(p) -> Tuple
as_ntuple(p, Val(N)) -> NTuple{N}

Normalize a point (Tuple, NamedTuple, or AbstractVector of length 1–3) to a plain Tuple.

A vector's length is a runtime value, so for one the first form returns a union of the three tuple widths. Every entry point here reduces a point to a value whose type does not depend on that width — distance to a scalar, spherical_to_cartesian to three components — so the union splits across the branch and each arm infers concretely and allocates nothing.

scale_factors is the exception, its result being one factor per direction. Name the width with Val(N) there, or hand it a tuple.

as_tensor6(τ) -> NTuple{6}

Normalize a symmetric rank-2 tensor to a plain 6-tuple of its independent components, from a Tuple, NamedTuple or AbstractVector — the rank-2 counterpart of as_ntuple.

build_point(S, names::NTuple{N,Symbol}, vals::NTuple{N}) -> S

Assemble a point/vector in representation S from positional vals and their names. Bare NamedTuple takes names; a parameterized NamedTuple type takes its own. Everything else is built from vals alone, covering Tuple, NTuple{N,T}, Vector{T}, SVector{N,T}, MVector{N,T}, and any user type with a positional constructor.

cartesian_to_geodetic(geo, xyz) -> NamedTuple{(:λ,:φ,:h)}
cartesian_to_geodetic(S, geo, xyz) -> S

Inverse of geodetic_to_cartesian: Earth-centred Cartesian to geodetic (λ, φ, h) with λ ∈ (-π, π], φ the geodetic latitude in [-π/2, π/2], and h the height above the ellipsoid.

Unlike the spherical inverse this is not closed-form — geodetic latitude satisfies a quartic — so it is solved by Newton's method on

tan φ = z(N+h) / (p(N(1-e²)+h)),   p = √(x²+y²)

with h = p·cosφ + z·sinφ − a√(1-e²sin²φ) substituted each step, which is exact for any φ and stays finite at the poles where p/cosφ − N does not. Convergence is quadratic and reaches eps(T) in three or four steps at terrestrial flattening. On the polar axis λ = 0 and the latitude is ±π/2 exactly.

The foot of the normal is unique for any point outside the ellipsoid's evolute — every position on or above the surface, and far below it. Deep inside, near the centre, a point has several feet and no preferred one; the centre itself is reported as (0, 0, -a), the same convention cartesian_to_spherical uses at the origin.

cartesian_to_spherical(geo, xyz) -> NamedTuple{(:λ,:φ,:r)}
cartesian_to_spherical(S, geo, xyz) -> S

Inverse of spherical_to_cartesian: Cartesian (x, y, z) to (λ, φ, r) with λ ∈ (-π, π], φ ∈ [-π/2, π/2], and r the absolute radius from the origin. At the origin λ = φ = 0.

distance(geo::AbstractEllipsoidalGeometry, p1, p2)

2-D (λ, φ): geodesic distance by Vincenty's inverse method, iterated to 1e-12 in the auxiliary longitude, which holds the result to well under a millimetre at Earth scale.

Vincenty's iteration converges slowly for very nearly antipodal point pairs. It is capped, and on reaching the cap the result is the great-circle distance on a sphere of the mean radius.

3-D (λ, φ, h): the chord through Cartesian (ECEF), mirroring the spherical 3-D convention. A lone (λ,) is the shorter arc of the equator.

distance(geo::AbstractGeometry, pt1, pt2)

Distance between two points. Spherical 2D uses great-circle (Haversine); spherical 3D uses the chord through Cartesian, with r the absolute radius from the origin; a lone (λ,) is the shorter arc of its circle.

First eccentricity squared, e² = f(2-f).

embed(geo, coords) -> NTuple{D,T}
embed(S, geo, coords) -> S

A point's position in the ambient Cartesian space this metric is realized in: the coordinates themselves on the plane, the point on the sphere of radius R, the Earth-centred position of a geodetic (λ, φ[, h]) on a spheroid.

One of the two primitives the metric-independent constructions here are written against; the other is local_tangent_basis. A straight-line displacement, a chord distance, a tangent-plane projection and a least-squares gradient need a position and a frame and nothing else, so each is written once and serves every geometry supplying those two.

This is the ambient position, always. What a spatial index searches is a separate question — an index may scale a sphere's surface coordinates to make an arc comparable to a chord, which is FlowGeometries.Grids.embed_point and its FlowGeometries.Grids.AbstractEmbedding.

flattening(geo) -> T

Flattening f = (a-b)/a. Defaults to the f field.

float_type(geo) -> Type{<:AbstractFloat}

The element type geo computes in. Anything that pairs a geometry with an optional element type defaults to this, so a geometry given without one keeps its own width.

geodetic_to_cartesian(geo, coords) -> (; x, y, z)

(λ, φ, h) — geodetic latitude and height above the ellipsoid — to Earth-centred Cartesian: x = (N(φ)+h)cosφ·cosλ, y = (N(φ)+h)cosφ·sinλ, z = (N(φ)(1-e²)+h)sinφ. The 2-D form takes h = 0, the surface.

jacobian(geo, point) -> T

∏ of scale_factors: the volume element per unit coordinate volume at point.

local_displacement(geo, center, neighbor) -> NamedTuple

neighbor's displacement from center, resolved on the local coordinate frame at center, one component per coordinate — the quantity a least-squares fit differences against.

Two coordinates is project_to_tangent_plane, which this calls: a surface has a tangent plane and the displacement lies in it. Three is the same ambient chord on the full local frame — (x, y, z) on the plane, and eastward/northward/outward on a sphere or spheroid, the third direction being the normal unit_vector gives.

local_tangent_basis(geo, coords) -> NamedTuple
local_tangent_basis(S, geo, coords) -> NamedTuple of `S`

Coordinate-aligned unit vectors at coords, in the ambient Cartesian frame:

• Cartesian: (; x = ê_x, y = ê_y), 2-component
• Spherical and spheroidal: (; λ = ê_λ, φ = ê_φ), 3-component — the eastward and northward tangents, which are the same triad on either (see AbstractLonLatGeometry)

The basis vectors themselves are Tuples by default; pass a leading S to get them as SVectors (or any other representation) so they can be used with vector arithmetic directly.

meridional_radius(geo, φ) -> T

M(φ) = a(1-e²) / (1 - e²sin²φ)^{3/2}: the radius of curvature along the meridian. A unit step in latitude covers M(φ).

metric_invariant_directions(geo) -> NTuple{K,Int}

The coordinate directions along which scale_factors does not vary.

A bulk operation that divides by a scale factor then solves it once per line: on a sphere no factor depends on longitude, so a whole row shares one value.

The default is (), no direction invariant, and each concrete geometry declares its own. A subtype of AbstractSphericalGeometry may write its own scale_factors, and an inherited claim about them hoists a value that varies, giving a wrong derivative. Declare the directions for your own geometry to opt in.

named_point(geo, vals::NTuple{N}) -> NamedTuple

Build the geometry-appropriate named point from positional values (axis order).

nonuniform_first_derivative(f_m, f_0, f_p, h_m, h_p)

Standard 3-point, 2nd-order-accurate centered finite-difference approximation of the first derivative at the middle node on a possibly nonuniform stencil. f_m, f_0, f_p are the function values at the minus/center/plus nodes, and h_m = x_0 - x_{-}, h_p = x_{+} - x_0 > 0 are the (physical) left/right spacings.

Reduces exactly to the uniform central difference (f_p - f_m)/(2h) when h_m == h_p == h, and is exact for linear and quadratic f for any h_m, h_p.

point_names(geo, Val(N)) -> NTuple{N,Symbol}

Field names for an N-coordinate point in this geometry. Cartesian: (:x,), (:x,:y), (:x,:y,:z), then (:x1, …, :xN). Spherical: (:λ,), (:λ,:φ), (:λ,:φ,:r), then (:λ,:φ,:r,:q4,…).

prime_vertical_radius(geo, φ) -> T

N(φ) = a / √(1 - e²sin²φ): the radius of curvature perpendicular to the meridian. A unit step in longitude covers N(φ)cosφ.

project_to_tangent_plane(geo, center, neighbor) -> NamedTuple
project_to_tangent_plane(S, geo, center, neighbor) -> S

Tangent-plane displacement of neighbor relative to center along the local coordinate basis from local_tangent_basis: the ambient chord between the two embedded positions, resolved on that basis.

A least-squares gradient and a scattered interpolation are built on this, so both work on any geometry supplying those two primitives — a spheroid included, where the chord runs between geodetic positions and the basis is the same triad.

radius(geo) -> T

The sphere radius. The one method a spherical geometry supplies; defaults to the R field.

rotate!(λ, φ, rot) -> (λ, φ)
unrotate!(λ, φ, rot) -> (λ, φ)

Rotate a whole point set in place — the form a sampling's spherical_points output takes. λ and φ are any arrays of matching shape, so this covers a scattered node set and a grid's 2-D coordinate fields alike. Allocates nothing.

rotate(rot, λ, φ) -> (λ′, φ′)

(λ, φ) expressed in the rotated frame. The rotation's own pole maps to φ′ = π/2.

scale_factors(geo, point) -> NTuple
scale_factors(geo, point, Val(N)) -> NTuple{N}

Physical length of a unit coordinate step in each direction at point — see FlowGeometries.Discretization.scale_factors.

Cartesian: 1 in every direction, the metric being the identity. Spherical: (R·cosφ, R) on the surface of radius R, and (r·cosφ, r, 1) where a radius direction is present, r being that point's own radius.

The result is one factor per direction, so from a point whose length is a runtime value — an AbstractVector — its width is one too. Val(N) names that width, checks it, and gives a concrete NTuple{N}.

semimajor_axis(geo) -> T

Equatorial radius a. Defaults to the a field.

Semi-minor axis b = a(1-f).

similar_geometry(T, geo) -> AbstractGeometry{T}

geo with its element type changed to T, keeping its shape parameters.

A geometry fixes the width of everything computed against it — coordinates, distances, metric factors — so building a grid at one element type around a geometry at another silently promotes the whole grid back. Anything that takes an element type and a geometry together should pass the geometry through here first.

A geometry already at T is returned unchanged, whatever its type, so a geometry defined outside this package inherits the whole stack without supplying anything. Only an actual change of width needs a method, and only a geometry that wants to support one has to define it.

similar_rotation(T, rot) -> PoleRotation{T}

rot at element width T, the counterpart of similar_geometry for a frame. The tilt is re-resolved at the new width from (λp, φp).

spherical_excess(a, b, c) -> T

Spherical excess of the triangle spanned by three unit vectors, via Van Oosterom & Strackee (1983):

tan(E/2) = |a · (b × c)| / (1 + a·b + b·c + c·a)

One atan and no other transcendental. Taking directions lets a mesh convert each vertex once — see unit_vector — and share it across every triangle. Multiply by R² for an area, as triangle_area does.

The numerator is evaluated as a · (u × v) for u = b - a and v = c - a, an identity because a · (a × x) vanishes for every x. u and v carry the triangle's own scale, so the result is never assembled by cancelling terms of order one. Cancellation there costs the excess all of its correct digits below a side of about R·√eps(T) — 2 km on Earth at Float32, which is an ordinary cell.

spherical_to_cartesian(geo, coords) -> NamedTuple{(:x,:y,:z)}
spherical_to_cartesian(S, geo, coords) -> S

(λ, φ) on the reference sphere of radius R, or (λ, φ, r) with r the absolute radius from the origin, mapped to Cartesian coordinates of the same space.

tensor_from_local(geo, τλλ, τφφ, τrr, τλφ, τλr, τφr, λ, φ) -> NamedTuple
tensor_from_local(geo, τ, λ, φ) -> NamedTuple
tensor_from_local(S, geo, args...) -> S

The inverse of tensor_to_local: τ = Rᵀ τ' R, returning the ambient Cartesian components (; xx, yy, zz, xy, xz, yz). Composing the two in either order is the identity.

tensor_to_local(geo, τxx, τyy, τzz, τxy, τxz, τyz, λ, φ) -> NamedTuple
tensor_to_local(geo, τ, λ, φ) -> NamedTuple
tensor_to_local(S, geo, args...) -> S

Rotate a symmetric rank-2 tensor from the ambient Cartesian frame into the local (ê_λ, ê_φ, ê_r) frame at (λ, φ): τ' = R τ Rᵀ, returning (; λλ, φφ, rr, λφ, λr, φr).

The rank-1 counterpart is vector_from_cartesian; between them they cover essentially every field anyone rotates — a stress, a strain rate, a covariance, a subfilter flux. τ may be given as the six independent components or as any 6-component representation, in the order above.

R is the same basis local_tangent_basis defines, so the two cannot drift apart — and being a rotation of the frame alone, it holds on a spheroid as it does on a sphere.

triangle_area(geo, p1, p2, p3) -> T

Exact area of the spherical triangle through three (λ, φ) points, each in any accepted representation. R² times spherical_excess.

This converts all three points on every call. A mesh that already holds vertex directions should call spherical_excess on those and scale once.

triangle_area_from_unit_vectors(geo, u1, u2, u3) -> T
triangle_area_from_unit_vectors(R², u1, u2, u3) -> T

triangle_area for vertices already held as unit vectors — R² times spherical_excess, with no coordinate conversion.

It carries its own name because dispatch cannot tell the two apart: a unit vector and a 3-D spherical point (λ, φ, r) are both 3-tuples.

Pass R² itself to walk many triangles: a mesh shares vertices between cells, so the squaring belongs outside the loop alongside the one-per-vertex unit_vector call. That is the form the cell-area and Voronoi paths here use.

unit_vector(T, p) -> NTuple{3,T}

The direction of the (λ, φ) point p on the unit sphere, in any accepted point representation. This is spherical_to_cartesian with the radius divided out, as a bare tuple: the form the spherical-triangle kernels want, computed once per vertex and reused across every triangle sharing it.

unrotate(rot, λ′, φ′) -> (λ, φ)

Inverse of rotate.

vector_from_cartesian(geo, ux, uy, uz, λ, φ) -> NamedTuple{(:λ,:φ,:r)}
vector_from_cartesian(geo, v, λ, φ) -> NamedTuple{(:λ,:φ,:r)}
vector_from_cartesian(S, geo, args...) -> S

Map an ambient Cartesian vector (ux, uy, uz) at (λ, φ) into local components (λ = u_λ, φ = u_φ, r = u_r) — eastward, northward, and along the outward normal. v may be any 3-component representation. The inverse of vector_to_cartesian, and like it the same on a sphere and on a spheroid.

vector_to_cartesian(geo, u_λ, u_φ, u_r, λ, φ) -> NamedTuple{(:x,:y,:z)}
vector_to_cartesian(geo, u_λ, u_φ, λ, φ) -> NamedTuple{(:x,:y,:z)}
vector_to_cartesian(S, geo, args...) -> S

Map the components of a vector given in the local (ê_λ, ê_φ, ê_r) basis at (λ, φ) into the ambient Cartesian basis. This transforms vector components at a point; see embed for positions. The 2-component form assumes u_r = 0.

A rotation of the frame and nothing else, so it is one and the same on a sphere and on a spheroid — see AbstractLonLatGeometry.

volume_element(geo::AbstractCartesianGeometry, dx, dy, dz)
volume_element(geo::AbstractSphericalGeometry, r, φ, dλ, dφ, dr)

Local cell volume from the cell's own extents. The spherical form takes the local radius r at this level, giving the shell element r²·cosφ·dλ·dφ·dr; the reference radius does not enter.

volume_element(geo::AbstractEllipsoidalGeometry, φ, h, dλ, dφ, dh)

(N(φ)+h)cosφ · (M(φ)+h) · dλ·dφ·dh, the geodetic volume element at ellipsoidal height h.

Unlike the spherical form this does not factor into a function of φ times a function of h: both curvature radii are offset by h, so the two directions are coupled.