Geometry
A geometry answers one question: what is the metric? It carries no points and no topology — only how to measure distance, area and volume, and what the coordinate directions are called.
using FlowGeometries: FlowGeometries as FG
cart = FG.Geometry.CartesianGeometry() # flat metric; carries no spacing
sph = FG.Geometry.SphericalGeometry() # default Earth radius
unit = FG.Geometry.SphericalGeometry(1.0) # unit sphere
wgs = FG.Geometry.SpheroidGeometry() # WGS 84 oblate spheroidFlowGeometries.Geometry.SpheroidGeometry{Float64}(6.378137e6, 0.0033528106647474805)Element type
A geometry is parameterized by the float type it computes in, and that type fixes the width of everything measured against it — coordinates, distances, metric factors. float_type reads it and similar_geometry carries a geometry to another width, keeping its shape:
FG.Geometry.float_type(sph), FG.Geometry.similar_geometry(Float32, wgs)(Float64, FlowGeometries.Geometry.SpheroidGeometry{Float32}(6.378137f6, 0.0033528106f0))Anything that builds at a chosen width takes the type FIRST, positionally, as zeros and rand do — a type given as a keyword takes no part in dispatch, so the result's element type would be unknown to the caller. Where a constructor takes both a width and a geometry, the width wins and the geometry is carried to it; where only a geometry is given, its own width is the one meant:
g32 = FG.Connectivity.structured_grid(Float32, FG.SphericalSampling.GaussLegendreSampling(), 16)
eltype(FG.Grids.axis(g32, 1)), FG.Geometry.float_type(FG.Grids.grid_geometry(g32))(Float32, Float32)Coordinate names
The geometry decides what a point's components are called, and the grid types inherit that. Cartesian is (x, y[, z]); spherical is (λ, φ[, r]) — longitude, geographic latitude, radius.
This is enforced rather than conventional: grid.x on a spherical grid is a FieldError, not a silent alias for longitude. Reading x and getting longitude is the kind of bug that survives code review and shows up as a wrong answer months later.
FG.Grids.coords(grid, 2, 3) # (λ = …, φ = …) on a spherical grid
FG.Grids.coordinate_names(grid) # (:λ, :φ)(:λ, :φ)Distance
distance is great-circle on a sphere and Euclidean on a plane, and takes points in whatever representation you have — tuples, NamedTuples, or (with the StaticArrays extension) static vectors.
FG.Geometry.distance(sph, (0.0, 0.0), (π/2, 0.0)) # quarter of the equator1.0007543398010286e7The spherical version uses the haversine form, which stays accurate for nearby points where the law-of-cosines form loses precision to cancellation.
Cell measures
FG.Geometry.area_element(sph, φ, Δλ, Δφ) # R²·cosφ·Δλ·Δφ
FG.Geometry.volume_element(sph, r, φ, Δλ, Δφ, Δr) # r²·cosφ·Δλ·Δφ·Δr
FG.Geometry.area_element(cart, 2.0, 3.0) # dx·dy, from the cell's own extents6.0These are the pointwise elements. A grid's per-cell measure is built from them once at construction — see Grids.
Cartesian ↔ spherical
p = FG.Geometry.spherical_to_cartesian(sph, (λ, φ)) # (x = …, y = …, z = …)
FG.Geometry.cartesian_to_spherical(sph, (p.x, p.y, p.z)) # (λ = …, φ = …, r = …)(λ = 0.29999999999999993, φ = 0.4, r = 6.371e6)Vectors transform differently from points, because a vector at a point is expressed in that point's local frame. Those have their own pair:
# components first, then the position they are expressed at
v = FG.Geometry.vector_to_cartesian(sph, 1.0, -0.5, 0.0, λ, φ) # (u_λ, u_φ, u_r) → (x, y, z)
FG.Geometry.vector_from_cartesian(sph, v.x, v.y, v.z, λ, φ) # and back(λ = 1.0, φ = -0.5, r = 0.0)Tangent-plane geometry
ê = FG.Geometry.local_tangent_basis(sph, (λ, φ)) # (; λ = ê_λ, φ = ê_φ)
FG.Geometry.project_to_tangent_plane(sph, (λ, φ), (λ + 1e-6, φ))(λ = 5.868079592699626, φ = 1.1424685396610457e-6)project_to_tangent_plane gives a neighbour's offset in the tangent plane at centre, which is what finite-difference and structure-function calculations on a sphere actually need.
Nonuniform derivatives
FG.Geometry.nonuniform_first_derivative(1.0, 0.0, 4.0, 1.0, 2.0)0.0Second-order accurate on an unequally spaced stencil — the usual three-point formula degrades to first order when h₋ ≠ h₊, which is the common case on a stretched grid.
A spheroid drives the whole stack
SpheroidGeometry overrides distance (Vincenty), area_element, volume_element and scale_factors; the grid, sampling and connectivity layers are inherited.
FG.Geometry.distance(wgs, (0.0, 0.0), (0.0, π/2)) # quarter meridian, WGS 841.0001965729311792e7FG.Geometry.prime_vertical_radius(wgs, 0.0), FG.Geometry.meridional_radius(wgs, π/2)(6.378137e6, 6.399593625758492e6)Grid directions are (λ, φ, h), with h the height above the ellipsoid — not the absolute radius that a spherical grid's third direction carries.
λ = range(0, 2π; length = 9)[1:8]
φ = range(-π/2, π/2; length = 7)
gs = FG.Grids.StructuredGrid(wgs, λ, φ)
FG.Grids.coordinate_names(gs), FG.Grids.isperiodic(gs, 1), size(gs)((:λ, :φ), true, (8, 7))The surface element M(φ)·N(φ)cosφ·Δλ·Δφ factors, so the measure stays separable and an interior cell matches the geometry's own element exactly:
Δλ, Δφ = FG.Grids.spacing(gs, 1), FG.Grids.spacing(gs, 2)
FG.Grids.measure(gs, 3, 4) ≈ FG.Geometry.area_element(wgs, φ[4], Δλ, Δφ)trueAdding a height direction changes that. The geodetic volume element offsets both curvature radii by h, coupling φ and h, so no product of per-axis factors reproduces it and the measure is stored dense — which measure_factors reports by returning nothing:
g3 = FG.Grids.StructuredGrid(wgs, λ, φ, range(0.0, 2000.0; length = 3))
FG.Grids.measure_factors(gs) !== nothing, FG.Grids.measure_factors(g3) === nothing(true, true)Δh = FG.Grids.spacing(g3, 3)
FG.Grids.measure(g3, 3, 4, 2) ≈ FG.Geometry.volume_element(wgs, φ[4], 1000.0, Δλ, Δφ, Δh)trueAdding a geometry
Subtype AbstractCartesianGeometry, AbstractSphericalGeometry or AbstractEllipsoidalGeometry and you inherit the whole grid, sampling and connectivity stack. Each hierarchy asks for its shape parameters through accessors, so that is all a new geometry has to supply — a sphere defines radius, an ellipsoid semimajor_axis and flattening, and no method reads a field, so store them however you like:
struct UnitSphere{T} <: FG.Geometry.AbstractSphericalGeometry{T} end
FG.Geometry.radius(::UnitSphere{T}) where {T} = one(T)
u = UnitSphere{Float64}()
gu = FG.Grids.StructuredGrid(u, λ, φ)
sum(FG.Grids.measure(gu)), FG.Geometry.distance(u, (0.0, 0.0), (0.0, π/2))(12.277955025156961, 1.5707963267948966)That total is the unit sphere's area to the discretization's accuracy, and the distance is a quarter great circle — both from the one method above.
Rotated frames
rot = FG.Geometry.PoleRotation(0.7, 0.3) # the frame whose north pole is at (0.7, 0.3)
FG.Geometry.rotate(rot, 0.7, 0.3) # that pole maps to φ = π/2(0.0, 1.5707963267948966)FG.Geometry.unrotate(rot, FG.Geometry.rotate(rot, 1.2, -0.4)...) # round-trips(1.2, -0.39999999999999986)A whole point set rotates in place — the shape a sampling's spherical_points output has — and the array forms allocate nothing:
λs = [0.1, 1.2, 3.0, 5.5]
φs = [0.0, -0.4, 0.9, 0.2]
FG.Geometry.rotate!(λs, φs, rot)
λs4-element Vector{Float64}:
5.120140589607384
0.6258905466991032
2.6523991954146995
4.545510895963549Rotating a rectilinear spherical grid gives a curvilinear one, because that is what it is: only its own frame's axes are separable. unrotate is the usual direction, taking a rotated-pole grid's (λ′, φ′) axes to the geographic coordinates of each cell.
sph = FG.Geometry.SphericalGeometry()
λr = range(0, 2π; length = 25)[1:24]
φr = range(-1.2, 1.2; length = 13)
grot = FG.Grids.unrotate(FG.Grids.StructuredGrid(sph, λr, φr), rot)
typeof(grot).name.name, size(grot), FG.Grids.coordinate_names(grot)(:CurvilinearGrid, (24, 13), (:λ, :φ))The cell measure carries over exactly rather than being recomputed, since a rotation is an isometry of the sphere — recomputing from the rotated corners would only add roundoff. The index topology carries over too: same mesh, same neighbours, so a direction that wrapped still wraps.
gplain = FG.Grids.StructuredGrid(sph, λr, φr)
FG.Grids.measure(grot, 3, 4) == FG.Grids.measure(gplain, 3, 4),
FG.Grids.isperiodic(grot, 1)(true, true)