-
Notifications
You must be signed in to change notification settings - Fork 24
/
octahedral_gaussian.jl
131 lines (109 loc) · 7.12 KB
/
octahedral_gaussian.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
"""An `OctahedralGaussianArray` is an array of octahedral grids, subtyping `AbstractReducedGridArray`,
that use Gaussian latitudes for each ring. First dimension of the underlying `N`-dimensional `data`
represents the horizontal dimension, in ring order (0 to 360˚E, then north to south),
other dimensions are used for the vertical and/or time or other dimensions.
The resolution parameter of the horizontal grid is `nlat_half` (number of latitude rings
on one hemisphere, Equator included) and the ring indices are precomputed in `rings`.
These grids are called octahedral because after starting with 20 points on the first ring around the
north pole (default) they increase the number of longitude points for each ring by 4, such that
they can be conceptually thought of as lying on the 4 faces of an octahedron on each hemisphere.
Hence, these grids have 20, 24, 28, ... longitude points for ring 1, 2, 3, ... There is no
ring on the Equator and the two rings around it have 16 + 4nlat_half longitude points before
reducing the number of longitude points per ring by 4 towards the southern-most ring
j = nlat. `rings` are the precomputed ring indices, the the example above
`rings = [1:20, 21:44, 45:72, ...]`. For efficient looping see `eachring` and `eachgrid`.
Fields are
$(TYPEDFIELDS)"""
struct OctahedralGaussianArray{T, N, ArrayType <: AbstractArray{T, N}} <: AbstractReducedGridArray{T, N, ArrayType}
data::ArrayType # data array, ring by ring, north to south
nlat_half::Int # number of latitudes on one hemisphere
rings::Vector{UnitRange{Int}} # TODO make same array type as data?
OctahedralGaussianArray(data::A, nlat_half, rings) where {A <: AbstractArray{T, N}} where {T, N} =
check_inputs(data, nlat_half, rings, OctahedralGaussianArray) ?
new{T, N, A}(data, nlat_half, rings) :
error_message(data, nlat_half, rings, OctahedralGaussianArray, T, N, A)
end
# TYPES
const OctahedralGaussianGrid{T} = OctahedralGaussianArray{T, 1, Vector{T}}
nonparametric_type(::Type{<:OctahedralGaussianArray}) = OctahedralGaussianArray
horizontal_grid_type(::Type{<:OctahedralGaussianArray}) = OctahedralGaussianGrid
full_array_type(::Type{<:OctahedralGaussianArray}) = FullGaussianArray
"""An `OctahedralGaussianArray` but constrained to `N=1` dimensions (horizontal only) and data is a `Vector{T}`."""
OctahedralGaussianGrid
# SIZE
nlat_odd(::Type{<:OctahedralGaussianArray}) = false
"""$(TYPEDSIGNATURES) [EXPERIMENTAL] additional number of longitude points on the first and last ring.
Change to 0 to start with 4 points on the first ring."""
npoints_pole(::Type{<:OctahedralGaussianArray}) = 16
"""$(TYPEDSIGNATURES) [EVEN MORE EXPERIMENTAL] number of longitude points added (removed) for every ring
towards the Equator (on the southern hemisphere towards the south pole)."""
npoints_added_per_ring(::Type{<:OctahedralGaussianArray}) = 4
function get_npoints2D(::Type{<:OctahedralGaussianArray}, nlat_half::Integer)
m, o = npoints_added_per_ring(OctahedralGaussianArray), npoints_pole(OctahedralGaussianArray)
return m*nlat_half^2 + (2o+m)*nlat_half
end
function get_nlat_half(::Type{<:OctahedralGaussianArray}, npoints2D::Integer)
m, o = npoints_added_per_ring(OctahedralGaussianArray), npoints_pole(OctahedralGaussianArray)
return round(Int, -(2o + m)/2m + sqrt(((2o+m)/2m)^2 + npoints2D/m))
end
function get_nlon_per_ring(Grid::Type{<:OctahedralGaussianArray}, nlat_half::Integer, j::Integer)
nlat = get_nlat(Grid, nlat_half)
@assert 0 < j <= nlat "Ring $j is outside O$nlat_half grid."
m, o = npoints_added_per_ring(OctahedralGaussianArray), npoints_pole(OctahedralGaussianArray)
j = j > nlat_half ? nlat - j + 1 : j # flip north south due to symmetry
return o + m*j
end
# maybe define at some point for Matrix(::OctahedralGaussianGrid)
# matrix_size(G::OctahedralGaussianGrid) = (2*(4+G.nlat_half), 2*(4+G.nlat_half+1))
## COORDINATES
get_colat(::Type{<:OctahedralGaussianArray}, nlat_half::Integer) = get_colat(FullGaussianArray, nlat_half)
function get_lon_per_ring(Grid::Type{<:OctahedralGaussianArray}, nlat_half::Integer, j::Integer)
nlon = get_nlon_per_ring(Grid, nlat_half, j)
return collect(0:2π/nlon:2π-π/nlon)
end
## QUADRATURE
get_quadrature_weights(::Type{<:OctahedralGaussianArray}, nlat_half::Integer) = gaussian_weights(nlat_half)
## INDEXING
"""$(TYPEDSIGNATURES) `UnitRange{Int}` to index grid points on ring `j` of `Grid`
at resolution `nlat_half`."""
function each_index_in_ring(Grid::Type{<:OctahedralGaussianArray},
j::Integer, # ring index north to south
nlat_half::Integer) # resolution param
nlat = get_nlat(Grid, nlat_half)
# TODO make m, o dependent
m, o = npoints_added_per_ring(OctahedralGaussianArray), npoints_pole(OctahedralGaussianArray)
m != 4 || o != 16 && @warn "This algorithm has not been generalised for m!=4, o!=16."
@boundscheck 0 < j <= nlat || throw(BoundsError) # ring index valid?
if j <= nlat_half # northern hemisphere incl Equator
index_1st = 2j*(j+7) - 15 # first in-ring index i
index_end = 2j*(j+9) # last in-ring index i
else # southern hemisphere excl Equator
j = nlat - j + 1 # mirror ring index around Equator
n = get_npoints2D(Grid, nlat_half) + 1 # number of grid points + 1
index_1st = n - 2j*(j+9) # count backwards
index_end = n - (2j*(j+7) - 15)
end
return index_1st:index_end # range of i's in ring
end
"""$(TYPEDSIGNATURES) precompute a `Vector{UnitRange{Int}} to index grid points on
every ring `j` (elements of the vector) of `Grid` at resolution `nlat_half`.
See `eachring` and `eachgrid` for efficient looping over grid points."""
function each_index_in_ring!( rings::Vector{<:UnitRange{<:Integer}},
Grid::Type{<:OctahedralGaussianArray},
nlat_half::Integer) # resolution param
nlat = length(rings)
@boundscheck nlat == get_nlat(Grid, nlat_half) || throw(BoundsError)
m, o = npoints_added_per_ring(OctahedralGaussianArray), npoints_pole(OctahedralGaussianArray)
index_end = 0
@inbounds for j in 1:nlat_half # North incl Eq only
index_1st = index_end + 1 # 1st index is +1 from prev ring's last index
index_end += o + m*j # add number of grid points per ring
rings[j] = index_1st:index_end # turn into UnitRange
end
@inbounds for (j, j_mirrored) in zip( nlat_half+1:nlat, # South only
nlat-nlat_half:-1:1) # reverse index
index_1st = index_end + 1 # 1st index is +1 from prev ring's last index
index_end += o + m*j_mirrored # add number of grid points per ring
rings[j] = index_1st:index_end # turn into UnitRange
end
end