Public Documentation

conformal_cubed_sphere_coordinates(Nx, Ny;
                                   spacing = UniformSpacing(),
                                   spacing_parameters = (; ratio_raised_to_Nx_minus_one = 10.5,
                                                           k₀ByNx = 0.45))

Generate computational-space coordinates x and y on [-1, 1] × [-1, 1] and map them to Cartesian coordinates (X, Y, Z) on the sphere using conformal_cubed_sphere_mapping. The arrays X, Y, and Z are of size (Nx, Ny) and correspond to Nx × Ny grid vertices, defining (Nx-1) × (Ny-1) spherical quadrilateral cells of a conformal cubed sphere panel.

If spacing = UniformSpacing() then x and y have uniform spacing (equal increments), so their tensor product defines a uniform grid on [-1, 1] × [-1, 1]. Otherwise, x and y both have a symmetric graded spacing. Specifically:

• spacing = GeometricSpacing(): uses geometric_spacing(N, ratio_raised_to_Nx_minus_one) for each axis.
• spacing = ExponentialSpacing(): uses exponential_spacing(N, k₀ByNx) for each axis.

Arguments

• Nx, Ny: Number of grid vertices along the x and y directions (≥ 2).

Keyword arguments

• spacing: Either UniformSpacing() (default), GeometricSpacing(), or ExponentialSpacing().
• spacing_parameters: Parameters required for various spacings. Default: (; ratio_raised_to_Nx_minus_one = 10.5, k₀ByNx = 0.45). For GeometricSpacing(), ratio_raised_to_Nx_minus_one is interpreted as r^(Nx-1); for ExponentialSpacing(), k₀ByNx is used as k₀ = k₀ByNx * Nx.

Returns

• x, y: Computational-space vertex coordinates of lengths Nx and Ny.
• X, Y, Z: (Nx, Ny) arrays of Cartesian coordinates on the sphere, with X[i, j], Y[i, j], Z[i, j] = conformal_cubed_sphere_mapping(x[i], y[j]).

conformal_cubed_sphere_inverse_mapping(X, Y, Z; Z_map=Z_Rancic)

Inverse mapping for conformal cube sphere for quadrant of North-pole face in which X and Y are both positive. All other mappings to other cube face coordinates can be recovered from rotations of this map. There a 3 other quadrants for the north-pole face and five other faces for a total of twenty-four quadrants. Because of symmetry only the reverse for a single quadrant is needed. Because of branch cuts and the complex transform the inverse mappings are multi-valued in general, using a single quadrant case allows a simple set of rules to be applied.

The mapping is valid for the cube face quadrant defined by $0 < x < 1$ and $0 < y < 1$, where a full cube face has extent $-1 < x < 1$ and $-1 < y < 1$. The quadrant for the mapping is from a cube face that has "north-pole" at its center $(x=0, y=0)$. i.e., has X, Y, Z = (0, 0, 1) at its center. The valid ranges of X and Y for this mapping and convention are a quadrant defined be geodesics that connect the points A, B, C and D, on the shell of a sphere of radius $R$ with X, Y coordinates as follows

A = (0, 0)
B = (√2, 0)
C = (√3/3, √3/3)
D = (0, √2)

conformal_cubed_sphere_mapping(x, y; W_map=W_Rancic)

Conformal mapping from a face of a cube onto the equivalent sector of a sphere with unit radius.

Map the north-pole face of a cube with coordinates $(x, y)$ onto the equivalent sector of the sphere with coordinates $(X, Y, Z)$.

The cube's face oriented normal to $z$-axis and its coordinates must lie within the range $-1 ≤ x ≤ 1$, $-1 ≤ y ≤ 1$ with its center at $(x, y) = (0, 0)$. The coordinates $X, Y$ increase in the same direction as $x, y$.

The numerical conformal mapping used here is described by Rančić et al. [1].

This is a Julia translation of MATLAB code from MITgcm that is based on Fortran 77 code from Jim Purser & Misha Rančić.

Examples

The center of the cube's face $(x, y) = (0, 0)$ is mapped onto $(X, Y, Z) = (0, 0, 1)$

julia> using CubedSphere

julia> conformal_cubed_sphere_mapping(0, 0)
(0, 0, 1.0)

and the edge of the cube's face at $(x, y) = (1, 1)$ is mapped onto $(X, Y, Z) = (√3/3, √3/3, √3/3)$

julia> using CubedSphere

julia> conformal_cubed_sphere_mapping(1, 1)
(0.5773502691896256, 0.5773502691896256, 0.5773502691896257)

References

• [1] Rančić et al., Q. J. R. Meteorol., (1996).

optimized_non_uniform_conformal_cubed_sphere_coordinates(Nx, Ny, spacing)

High-level driver that uses EKI to optimize the non-uniform spacing parameter for a conformal cubed sphere panel, then builds and returns the corresponding grid.

Procedure:

1. Create a random ensemble of parameters within limits (N_ensemble = 40, reproducible seed).
2. Run optimize! to fit the weighted diagnostics to their ideal targets.
3. Build the optimized grid with conformal_cubed_sphere_coordinates(Nx, Ny; spacing, ...) using the mean optimized parameter.

For GeometricSpacing(), the parameter is ratio^(Nx-1); for ExponentialSpacing(), the parameter is k₀/Nx.

Arguments

• Nx, Ny: Number of grid vertices along panel coordinates.
• spacing: GeometricSpacing() or ExponentialSpacing().

Returns

• x, y: Computational-space coordinates of lengths Nx and Ny.
• X, Y, Z: (Nx, Ny) Cartesian coordinates of the vertices of the optimized non-uniform conformal cubed sphere panel.

cartesian_to_lat_lon(x, y, z)
cartesian_to_lat_lon(X)

Convert 3D Cartesian coordinates (x, y, z) or a 3-element Cartesian vector X = (x, y, z) on the sphere to latitude–longitude. Returns a tuple (latitude, longitude) in degrees.

• Latitude is the angle measured from the equatorial plane (z = 0).
• Longitude is measured anti-clockwise (eastward) from the x-axis (y = 0) about the z-axis.

Arguments

• x, y, z: Cartesian coordinates (numbers), or
• X: 3-element Cartesian vector.

Returns

• (latitude, longitude): Latitude and longitude angles in degrees.

Examples

Find latitude–longitude of the North Pole:

julia> using CubedSphere.SphericalGeometry

julia> x, y, z = (0, 0, 6.4e6);  # Cartesian coordinates of the North Pole [in meters]

julia> cartesian_to_lat_lon(x, y, z)
(90.0, 0.0)

Let's confirm that for few points on the unit sphere we get the answers we expect.

julia> cartesian_to_lat_lon(√2/4, -√2/4, √3/2)
(59.99999999999999, -45.0)

julia> cartesian_to_lat_lon(-√6/4, √2/4, -√2/2)
(-45.00000000000001, 150.0)

cartesian_to_latitude(x, y, z)

Convert Cartesian coordinates (x, y, z) to latitude (in degrees) on the sphere.

cartesian_to_longitude(x, y, z)

Convert Cartesian coordinates (x, y, z) to longitude (in degrees) on the sphere.

compute_cell_areas(X, Y, Z; radius = 1)

Compute the spherical surface areas of all quadrilateral cells in a spherical grid (e.g., a conformal cubed-sphere panel), defined by the Cartesian coordinate arrays X, Y, and Z. Each of X, Y, and Z is a 2D array of size (Nx, Ny) holding the Cartesian coordinates of the grid vertices on the sphere, such that the point at (i, j) is (X[i, j], Y[i, j], Z[i, j]). The grid therefore contains (Nx−1) × (Ny−1) spherical quadrilateral cells.

For each quadrilateral cell with vertices (i, j), (i+1, j), (i+1, j+1), and (i, j+1), the function computes the spherical area using spherical_area_quadrilateral and stores the result in a 2D array. If radius = 1, the returned areas correspond to the unit sphere. For radius ≠ 1, the physical areas are returned.

Arguments

• X, Y, Z: (Nx, Ny) arrays of Cartesian coordinates of grid vertices on the sphere.
• radius: Sphere radius (optional). Default is 1.

Returns

• cell_areas: An (Nx−1, Ny−1) array of quadrilateral cell areas.
  • For radius = 1, these are unit-sphere areas.
  • For radius ≠ 1, these are physical areas (e.g., in m²).

Notes

• The function assumes that X, Y, and Z represent vertices lying on the same sphere.

Examples

julia> using CubedSphere.SphericalGeometry

julia> Nx, Ny = 3, 3;

julia> lons = range(-π/6, π/6, length = Nx);

julia> lats = range(-π/6,  π/6, length = Ny);

julia> X = [cos(φ)*cos(λ) for λ in lons, φ in lats];

julia> Y = [cos(φ)*sin(λ) for λ in lons, φ in lats];

julia> Z = [sin(φ)        for λ in lons, φ in lats];

julia> A = compute_cell_areas(X, Y, Z);

julia> size(A)
(2, 2)

julia> all(isapprox.(A, fill(A[1, 1], 2, 2); rtol=1e-12))
true

julia> isapprox(A[1, 1], 0.26636308214247195; rtol=1e-12)
true

julia> R = 6.371e6;

julia> A_R = compute_cell_areas(X .* R, Y .* R, Z .* R; radius = R);

julia> isapprox.(A_R, A .* R^2; rtol=1e-9) |> all
true

lat_lon_to_cartesian(φ, λ; radius = 1)

Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinates (x, y, z) on the sphere.

Arguments

• φ: Latitude in degrees.
• λ: Longitude in degrees.
• radius: Sphere radius (optional). Default is 1.

Returns

• (x, y, z): Cartesian coordinates on the sphere.

Examples

Find the Cartesian coordinates of the North Pole on a unit sphere:

julia> using CubedSphere.SphericalGeometry

julia> lat_lon_to_cartesian(90, 0)
(0.0, 0.0, 1.0)

Find the Cartesian coordinates of a point on the equator with longitude 90°E:

julia> lat_lon_to_cartesian(0, 90)
(0.0, 1.0, 0.0)

latlonto_x(φ, λ; radius = 1)

Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate x on the sphere.

latlonto_y(φ, λ; radius = 1)

Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate y on the sphere.

latlonto_z(φ; radius = 1)

Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate z on the sphere.

spherical_area_quadrilateral(a₁, a₂, a₃, a₄; radius = 1)

Compute the area of a spherical quadrilateral on a sphere of radius, given the position of its four vertices as vectors a₁, a₂, a₃, and a₄ in 3D Cartesian coordinates. The origin is assumed to be at the center of the sphere.

The quadrilateral area is evaluated as half the sum of the areas of all four spherical triangles formed by the vertices. This approach avoids the need to explicitly choose a diagonal that splits the quadrilateral into two non-overlapping triangles.

Arguments

• a₁, a₂, a₃, a₄: 3-element Cartesian vectors representing the four vertices of the spherical quadrilateral. All four must lie on the same sphere.
• radius: Sphere radius (optional). Default is 1.

Returns

• The area of the spherical quadrilateral on a sphere of radius.

Note

• This method is numerically robust and works for convex spherical quadrilaterals without requiring explicit diagonal selection.

Examples

julia> using CubedSphere.SphericalGeometry

julia> a₁ = [1.0, 0.0, 0.0];
       a₂ = [0.0, 1.0, 0.0];
       a₃ = [0.0, 0.0, 1.0];
       a₄ = [1.0, 1.0, 0.0] ./ √2;  # mid-edge on unit sphere

julia> spherical_area_quadrilateral(a₁, a₂, a₃, a₄)
1.5707963267948966

julia> R = 6.371e6;

julia> spherical_area_quadrilateral(a₁ .* R, a₂ .* R, a₃ .* R, a₄ .* R; radius = R)
6.375805898872353e13

spherical_area_triangle(a₁, a₂, a₃; radius = 1)

Compute the area of a spherical triangle on a sphere of radius, given its three vertex position vectors a₁, a₂, and a₃ in 3D Cartesian coordinates. The origin is assumed to be at the center of the sphere.

For a unit sphere (radius = 1), the area equals the spherical excess E. For a sphere of radius R, the area is R² * E.

Arguments

• a₁, a₂, a₃: 3-element Cartesian vectors representing the vertices of the spherical triangle. All three must lie on the same sphere.
• radius: Sphere radius (optional). Default is 1.

Returns

• The area of the spherical triangle on a sphere of radius.

Notes

• This function generalizes the classical Euler–Lagrange formula for spherical excess by expressing the quantity

  \[P = \sqrt{1 - \cos^2 a - \cos^2 b - \cos^2 c + 2 \cos a \cos b \cos c}\]

  in terms of the scalar triple product

  \[P = |𝐚₁ ⋅ (𝐚₂ × 𝐚₃)|\]

  where a, b, and c are the side lengths of the spherical triangle formed by the vertices a₁, a₂, and a₃. The above formula was first derived by Eriksson (1990).

• The inputs a₁, a₂, and a₃ must have the same norm.

References

• Euler, L. (1778) De mensura angulorum solidorum, Opera omnia, 26, 204-233 (Orig. in Acta adac. sc. Petrop. 1778)
• Lagrange, J.-L. (1798) Solutions de quelques problèmes relatifs au triangles sphériques, Oeuvres, 7, 331-359.
• Eriksson, F. (1990) On the measure of solid angles, Mathematics Magazine, 63 (3), 184-187,

doi:10.1080/0025570X.1990.11977515

Examples

julia> using CubedSphere.SphericalGeometry

julia> a₁ = [1.0, 0.0, 0.0];
       a₂ = [0.0, 1.0, 0.0];
       a₃ = [0.0, 0.0, 1.0];

julia> spherical_area_triangle(a₁, a₂, a₃)
1.5707963267948966

julia> spherical_area_triangle(a₁ .* 6.371e6, a₂ .* 6.371e6, a₃ .* 6.371e6; radius = 6.371e6)
6.375805898872353e13

spherical_area_triangle(a::Number, b::Number, c::Number; radius = 1)

Compute the area of a spherical triangle on a sphere of radius, given its three side lengths a, b, and c (in radians).

For a unit sphere (radius = 1), the area equals the spherical excess E = A + B + C − π, where A, B, and C are the triangle’s interior angles. For a sphere of radius R, the area is R² * E.

Arguments

• a, b, c: Side lengths of the spherical triangle, measured as central angles (in radians).
• radius: Sphere radius (optional). Default is 1.

Returns

• The area of the spherical triangle on a sphere of radius.

Notes

• Euler (1778) and Lagrange (1798) showed that the spherical excess E on the unit sphere is computed as

  \[\tan\frac{E}{2} = \frac{\sqrt{1 - \cos^2 a - \cos^2 b - \cos^2 c + 2 \cos a \cos b \cos c}}{1 + \cos a + \cos b + \cos c}.\]

References

• Euler, L. (1778) De mensura angulorum solidorum, Opera omnia, 26, 204-233 (Orig. in Acta adac. sc. Petrop. 1778)
• Lagrange, J.-L. (1798) Solutions de quelques problèmes relatifs au triangles sphériques, Oeuvres, 7, 331-359.

Examples

julia> using CubedSphere.SphericalGeometry

julia> a = b = c = π/2;  # Right spherical triangle with 90° sides on unit sphere

julia> spherical_area_triangle(a, b, c)
1.5707963267948966

julia> spherical_area_triangle(a, b, c; radius = 6371e3)
6.375805898872353e13

spherical_distance(a₁, a₂; radius = 1)

Compute the great-circle distance between two Cartesian points a₁ and a₂ on a sphere of radius.

For radius = 1, this is equivalent to the central angle (in radians) between the two points.

Arguments

• a₁, a₂: 3-element Cartesian vectors on the sphere.
• radius: Sphere radius (optional). Default is 1.

Returns

• Great-circle distance between a₁ and a₂

Notes

• Both a₁ and a₂ must lie on the surface of the same sphere (i.e., have the same norm).

Examples

julia> using CubedSphere.SphericalGeometry

julia> a₁ = (1.0, 0.0, 0.0);  # point on unit sphere
       a₂ = (0.0, 1.0, 0.0);  # 90° away

julia> spherical_distance(a₁, a₂)
1.5707963267948968

spherical_quadrilateral_vertices(X, Y, Z, i, j)

Returns the four Cartesian vertex vectors of the spherical grid cell whose corners are indexed by (i, j), (i+1, j), (i+1, j+1), and (i, j+1) in the arrays X, Y, and Z. Each of X, Y, and Z is a 2D array of size (Nx, Ny) holding the Cartesian coordinates of grid vertices on the sphere, such that the point at (i, j) is (X[i, j], Y[i, j], Z[i, j]).

Arguments

• X, Y, Z: (Nx, Ny) arrays of Cartesian coordinates on the sphere.
• i, j: Indices of the lower-left corner of the cell (in array order).

Returns

• (a₁, a₂, a₃, a₄): The four 3-element Cartesian vertex vectors at (i, j), (i+1, j), (i+1, j+1), and (i, j+1).