Public Documentation
Documentation for CubedSphere.jl
's public interface.
See the Internals section of the manual for internal package docs covering all submodules.
CubedSphere
CubedSphere.conformal_cubed_sphere_coordinates
— Methodconformal_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()
: usesgeometric_spacing(N, ratio_raised_to_Nx_minus_one)
for each axis.spacing = ExponentialSpacing()
: usesexponential_spacing(N, k₀ByNx)
for each axis.
Arguments
Nx, Ny
: Number of grid vertices along thex
andy
directions (≥ 2).
Keyword arguments
spacing
: EitherUniformSpacing()
(default),GeometricSpacing()
, orExponentialSpacing()
.spacing_parameters
: Parameters required for various spacings. Default:(; ratio_raised_to_Nx_minus_one = 10.5, k₀ByNx = 0.45)
. ForGeometricSpacing()
,ratio_raised_to_Nx_minus_one
is interpreted asr^(Nx-1)
; forExponentialSpacing()
,k₀ByNx
is used ask₀ = k₀ByNx * Nx
.
Returns
x
,y
: Computational-space vertex coordinates of lengthsNx
andNy
.X
,Y
,Z
:(Nx, Ny)
arrays of Cartesian coordinates on the sphere, withX[i, j], Y[i, j], Z[i, j] = conformal_cubed_sphere_mapping(x[i], y[j])
.
CubedSphere.conformal_cubed_sphere_inverse_mapping
— Methodconformal_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)
CubedSphere.conformal_cubed_sphere_mapping
— Methodconformal_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).
CubedSphere.optimized_non_uniform_conformal_cubed_sphere_coordinates
— Methodoptimized_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:
- Create a random ensemble of parameters within limits (
N_ensemble = 40
, reproducible seed). - Run
optimize!
to fit the weighted diagnostics to their ideal targets. - 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()
orExponentialSpacing()
.
Returns
x, y
: Computational-space coordinates of lengthsNx
andNy
.X, Y, Z
:(Nx, Ny)
Cartesian coordinates of the vertices of the optimized non-uniform conformal cubed sphere panel.
CubedSphere.SphericalGeometry
CubedSphere.SphericalGeometry.cartesian_to_lat_lon
— Methodcartesian_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 thez
-axis.
Arguments
x, y, z
: Cartesian coordinates (numbers), orX
: 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)
CubedSphere.SphericalGeometry.cartesian_to_latitude
— Methodcartesian_to_latitude(x, y, z)
Convert Cartesian coordinates (x, y, z)
to latitude (in degrees) on the sphere.
CubedSphere.SphericalGeometry.cartesian_to_longitude
— Methodcartesian_to_longitude(x, y, z)
Convert Cartesian coordinates (x, y, z)
to longitude (in degrees) on the sphere.
CubedSphere.SphericalGeometry.compute_cell_areas
— Methodcompute_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 is1
.
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²).
- For
Notes
- The function assumes that
X
,Y
, andZ
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
CubedSphere.SphericalGeometry.lat_lon_to_cartesian
— Methodlat_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 is1
.
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)
CubedSphere.SphericalGeometry.lat_lon_to_x
— Methodlatlonto_x(φ, λ; radius = 1)
Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate x on the sphere.
CubedSphere.SphericalGeometry.lat_lon_to_y
— Methodlatlonto_y(φ, λ; radius = 1)
Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate y on the sphere.
CubedSphere.SphericalGeometry.lat_lon_to_z
— Methodlatlonto_z(φ; radius = 1)
Convert (latitude, longitude) coordinates (in degrees) to Cartesian coordinate z on the sphere.
CubedSphere.SphericalGeometry.spherical_area_quadrilateral
— Methodspherical_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 is1
.
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
CubedSphere.SphericalGeometry.spherical_area_triangle
— Methodspherical_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 is1
.
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
, andc
are the side lengths of the spherical triangle formed by the verticesa₁
,a₂
, anda₃
. The above formula was first derived by Eriksson (1990).The inputs
a₁
,a₂
, anda₃
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
CubedSphere.SphericalGeometry.spherical_area_triangle
— Methodspherical_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 is1
.
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
CubedSphere.SphericalGeometry.spherical_distance
— Methodspherical_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 is1
.
Returns
- Great-circle distance between
a₁
anda₂
Notes
- Both
a₁
anda₂
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
CubedSphere.SphericalGeometry.spherical_quadrilateral_vertices
— Methodspherical_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)
.